Repository: dlr-eoc/ukis-pysat

Issue:
Build array subset from window tiles
Provide a method to build an array subset from window tiles.
The current workflow is to call the subset window in a for-loop through the [get_tiles](https://github.com/dlr-eoc/ukis-pysat/blob/master/ukis_pysat/data.py#L280) method. Subsetting the array is currently performed by the user by indexing the array with the window properties:
```
for window in data.Image(img_path).get_tiles():
    (row_start, row_stop), (col_start, col_stop) = window_tile.toranges()
    array[0][row_start: row_stop, col_start: col_stop]
```

This can be outsourced to a new method in [ukis_pysat.data](https://github.com/dlr-eoc/ukis-pysat/blob/master/ukis_pysat/data.py):
```
def get_subset(self, window_tile, array, band=0):
    (row_start, row_stop), (col_start, col_stop) = window_tile.toranges()
    return array[band][row_start: row_stop, col_start: col_stop]
```

Performance decreases when requesting the image array for each call:
```
def get_subset(self, window_tile, band=0):
    (row_start, row_stop), (col_start, col_stop) = window_tile.toranges()
    return self.arr[band][row_start: row_stop, col_start: col_stop]
```

This can be solved by handing the array to the new method directly. However, this requires refactoring the [constructor ](https://github.com/dlr-eoc/ukis-pysat/blob/master/ukis_pysat/data.py#L23) of the Image class so that arrays are accepted without handing a dataset as it is currently done [here](https://github.com/dlr-eoc/ukis-pysat/blob/master/ukis_pysat/data.py#L31).

Possible solution:
```
if isinstance(dataset, rasterio.io.DatasetReader) and isinstance(arr, np.ndarray):
    self.dataset = dataset
    self.arr = arr
elif isinstance(arr, np.ndarray):
    self.arr = arr
```


Discussion:

Patch:
diff --git a/ukis_pysat/data.py b/ukis_pysat/data.py
index 5c9f8ca..68a1847 100644
--- a/ukis_pysat/data.py
+++ b/ukis_pysat/data.py
@@ -278,7 +278,7 @@ class Image:
         return [wave_bands[platform][wavelength.lower()] for wavelength in wavelengths]
 
     def get_tiles(self, width=256, height=256, overlap=0):
-        """calculates rasterio.windows.Window, idea from https://stackoverflow.com/a/54525931
+        """Calculates rasterio.windows.Window, idea from https://stackoverflow.com/a/54525931
 
         TODO boundless with 'reflect' padding
 
@@ -306,6 +306,17 @@ class Image:
         # TODO using tukey --> https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.windows.tukey.html
         pass
 
+    def get_subset(self, tile, band=0):
+        """Get slice of array.
+
+        :param tile: rasterio.windows.Window tile from get_tiles().
+        :param band: Band number (default: 0).
+        :return: Sliced numpy array, bounding box of array slice.
+        """
+        # access window bounds
+        bounds = windows.bounds(tile, self.dataset.transform)
+        return self.arr[(band,) + tile.toslices()], bounds  # Shape of array is announced with (bands, height, width)
+
     def get_dask_array(self, chunk_size=(1, 6000, 6000)):
         """ transforms numpy to dask array
 


Test patch:
diff --git a/tests/test_data.py b/tests/test_data.py
index a7f7f87..99a36c8 100644
--- a/tests/test_data.py
+++ b/tests/test_data.py
@@ -93,6 +93,13 @@ class DataTest(unittest.TestCase):
 
         self.assertEqual(idx, 20807)
 
+    def test_get_subset(self):
+        for idx, each in enumerate(img.get_tiles(5, 5, 1)):
+            if idx == 2578:
+                array, bounds = img.get_subset(each)
+                self.assertTrue(np.array_equal(array, np.zeros(shape=(7, 7), dtype=array.dtype)))
+                self.assertEqual(bounds, (11.903960582768779, 51.45624717410995, 11.904589403469808, 51.45687599481152))
+
     def test_get_dask_array(self):
         self.assertIsInstance(img.get_dask_array(chunk_size=(1, 10, 10)), dask.array.core.Array)