Repository: dlr-eoc/ukis-pysat

Issue:
Image coordinate order
We use the coordinate order of rasterio when loading images, which is channels first (channels, rows, cols). Most image processing, visualization and machine learning libraries, however, use channels last (rows, cols, channels) for image coordinate convention (e.g. scikit-image, opencv, pillow, matplotlib, etc.).  

In order to avoid confusion and to adjust to user requirements more easily it would be great if we could provide an option to change the coordinate order upon instantiation of an image object. 

I know that Rasterio provides conversion methods [here](https://rasterio.readthedocs.io/en/latest/topics/image_processing.html?highlight=read#image-processing-software), but it would be convenient if the order could be set once (e.g. as env var or parameter in Image init) to comply with a user's requirements.

Tensorflow for example provides such a switch between channels first and channels last to account for performance differences when training on GPU or CPU.

Discussion:

Patch:
diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index 4f8ab01..bd7d7fd 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -18,6 +18,8 @@ something was changed
 
 Added
 *****
+- ``raster.Image()``: expanded test_arr to test AttributeError #31
+- ``raster.Image()``: optional dimorder for arr and according test #31
 - dn2toa tests and testfiles #17
 
 Fixed
diff --git a/README.md b/README.md
index 87d2ed4..3eeb6fd 100644
--- a/README.md
+++ b/README.md
@@ -36,7 +36,7 @@ from ukis_pysat.raster import Image
 
 
 # connect to Scihub and query metadata (returns MetadataCollection)
-src = Source(Datahub.Scihub)
+src = Source(source=Datahub.Scihub)
 meta = src.query_metadata(
     platform=Platform.Sentinel2,
     date=("20200101", "NOW"),
diff --git a/ukis_pysat/raster.py b/ukis_pysat/raster.py
index bcb97ca..398e040 100644
--- a/ukis_pysat/raster.py
+++ b/ukis_pysat/raster.py
@@ -10,6 +10,7 @@ try:
     import rasterio.mask
     from rasterio import windows
     from rasterio.io import MemoryFile
+    from rasterio.plot import reshape_as_image
     from rasterio.warp import calculate_default_transform, reproject
     from rio_toa import reflectance, brightness_temp, toa_utils
     from shapely.geometry import box, polygon
@@ -27,32 +28,42 @@ logger = logging.getLogger(__name__)
 
 
 class Image:
-    def __init__(self, path=None, dataset=None, arr=None):
+    def __init__(self, path=None, dataset=None, arr=None, dimorder="first"):
         if path:
             if isinstance(path, str):
                 self.dataset = rasterio.open(path)
-                self.arr = self.dataset.read()
+                self.__arr = self.dataset.read()
             else:
                 raise TypeError(f"path must be of type str")
         else:
             if isinstance(dataset, rasterio.io.DatasetReader) and isinstance(arr, np.ndarray):
                 self.dataset = dataset
-                self.arr = arr
+                self.__arr = arr
             else:
                 raise TypeError(
                     f"dataset must be of type rasterio.io.DatasetReader and arr must be of type " f"numpy.ndarray"
                 )
+        self.dimorder = dimorder
         self.transform = self.dataset.transform
         self.crs = self.dataset.crs
         self.da_arr = None
 
+    @property
+    def arr(self):
+        if self.dimorder == "first":
+            return self.__arr
+        elif self.dimorder == "last":
+            return reshape_as_image(self.__arr)
+        else:
+            raise AttributeError("dimorder for bands or channels must be either 'first' or 'last'.")
+
     def get_valid_data_bbox(self, nodata=0):
         """bounding box covering the input array's valid data pixels.
 
         :param nodata: nodata value, optional (default: 0)
         :return: tuple with valid data bounds
         """
-        valid_data_window = windows.get_data_window(self.arr, nodata=nodata)
+        valid_data_window = windows.get_data_window(self.__arr, nodata=nodata)
         return windows.bounds(valid_data_window, windows.transform(valid_data_window, self.transform))
 
     def mask_image(self, bbox, crop=True, pad=False, **kwargs):
@@ -66,9 +77,9 @@ class Image:
             self.dataset = self._pad_to_bbox(bbox, **kwargs).open()
 
         if isinstance(bbox, polygon.Polygon):
-            self.arr, self.transform = rasterio.mask.mask(self.dataset, [bbox], crop=crop)
+            self.__arr, self.transform = rasterio.mask.mask(self.dataset, [bbox], crop=crop)
         elif isinstance(bbox, tuple):
-            self.arr, self.transform = rasterio.mask.mask(self.dataset, [box(*bbox)], crop=crop)
+            self.__arr, self.transform = rasterio.mask.mask(self.dataset, [box(*bbox)], crop=crop)
         else:
             raise TypeError(f"bbox must be of type tuple or Shapely Polygon")
 
@@ -95,24 +106,25 @@ class Image:
         pad_width = math.ceil(max_diff / self.transform.to_gdal()[1])  # units / pixel_size
 
         destination = np.zeros(
-            (self.dataset.count, self.arr.shape[1] + 2 * pad_width, self.arr.shape[2] + 2 * pad_width,), self.arr.dtype,
+            (self.dataset.count, self.__arr.shape[1] + 2 * pad_width, self.__arr.shape[2] + 2 * pad_width,),
+            self.__arr.dtype,
         )
 
         for i in range(0, self.dataset.count):
             destination[i], self.transform = rasterio.pad(
-                self.arr[0], self.transform, pad_width, mode, constant_values=constant_values,
+                self.__arr[0], self.transform, pad_width, mode, constant_values=constant_values,
             )
 
-        self.arr = destination
+        self.__arr = destination
 
         mem_profile = self.dataset.meta
         mem_profile.update(
-            {"height": self.arr.shape[-2], "width": self.arr.shape[-1], "transform": self.transform,}
+            {"height": self.__arr.shape[-2], "width": self.__arr.shape[-1], "transform": self.transform,}
         )
 
         memfile = MemoryFile()
         ds = memfile.open(**mem_profile)
-        ds.write(self.arr)
+        ds.write(self.__arr)
 
         return memfile
 
@@ -141,7 +153,7 @@ class Image:
                 self.dataset.crs, dst_crs, self.dataset.width, self.dataset.height, *self.dataset.bounds,
             )
 
-        destination = np.zeros((self.dataset.count, height, width), self.arr.dtype)
+        destination = np.zeros((self.dataset.count, height, width), self.__arr.dtype)
 
         for i in range(0, self.dataset.count):
             reproject(
@@ -156,7 +168,7 @@ class Image:
             )
 
         # update for further processing
-        self.arr = destination
+        self.__arr = destination
         self.transform = transform
         self.crs = dst_crs
 
@@ -199,7 +211,7 @@ class Image:
                         # rescale thermal bands
                         toa.append(
                             brightness_temp.brightness_temp(
-                                self.arr[idx, :, :],
+                                self.__arr[idx, :, :],
                                 ML=multiplicative_rescaling_factors,
                                 AL=additive_rescaling_factors,
                                 K1=thermal_conversion_constant1,
@@ -213,16 +225,16 @@ class Image:
                     additive_rescaling_factors = metadata["RADIOMETRIC_RESCALING"][f"REFLECTANCE_ADD_BAND_{b}"]
                     toa.append(
                         reflectance.reflectance(
-                            self.arr[idx, :, :],
+                            self.__arr[idx, :, :],
                             MR=multiplicative_rescaling_factors,
                             AR=additive_rescaling_factors,
                             E=sun_elevation,
                         )
                     )
 
-                self.arr = np.array(np.stack(toa, axis=0))
+                self.__arr = np.array(np.stack(toa, axis=0))
         elif platform == Platform.Sentinel2:
-            self.arr = self.arr.astype(np.float32) / 10000.0
+            self.__arr = self.__arr.astype(np.float32) / 10000.0
         else:
             logger.warning(
                 f"Cannot convert dn2toa. Platform {platform} not supported [Landsat-5, Landsat-7, Landsat-8, "
@@ -286,8 +298,8 @@ class Image:
         :yields: window of tile
         """
         logger.info(f"Get tiles with size {width}x{height}.")
-        rows = self.arr.shape[-2]
-        cols = self.arr.shape[-1]
+        rows = self.__arr.shape[-2]
+        cols = self.__arr.shape[-1]
         offsets = product(range(0, cols, width), range(0, rows, height))
         bounding_window = windows.Window(col_off=0, row_off=0, width=cols, height=rows)
         for col_off, row_off in offsets:
@@ -313,7 +325,7 @@ class Image:
         """
         # 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)
+        return self.__arr[(band,) + tile.toslices()], bounds  # Shape of array is announced with (bands, height, width)
 
     def to_dask_array(self, chunk_size=(1, 6000, 6000)):
         """ transforms numpy to dask array
@@ -326,7 +338,7 @@ class Image:
         except ImportError:
             raise ImportError("to_dask_array requires optional dependency dask[array].")
 
-        self.da_arr = da.from_array(self.arr, chunks=chunk_size)
+        self.da_arr = da.from_array(self.__arr, chunks=chunk_size)
         return self.da_arr
 
     def write_to_file(self, path_to_file, dtype=rasterio.uint16, driver="GTiff"):
@@ -340,8 +352,8 @@ class Image:
         profile.update(
             {
                 "driver": driver,
-                "height": self.arr.shape[-2],
-                "width": self.arr.shape[-1],
+                "height": self.__arr.shape[-2],
+                "width": self.__arr.shape[-1],
                 "dtype": dtype,
                 "transform": self.transform,
                 "crs": self.crs,
@@ -349,7 +361,7 @@ class Image:
         )
 
         with rasterio.open(path_to_file, "w", **profile) as dst:
-            dst.write(self.arr.astype(dtype))
+            dst.write(self.__arr.astype(dtype))
 
     def close(self):
         """closes Image"""


Test patch:
diff --git a/tests/test_raster.py b/tests/test_raster.py
index ef9c99b..0fbe9e1 100644
--- a/tests/test_raster.py
+++ b/tests/test_raster.py
@@ -6,6 +6,7 @@ import numpy as np
 from rasterio import windows
 from rasterio.coords import BoundingBox
 from shapely.geometry import box
+from rasterio.transform import from_bounds
 
 from ukis_pysat.members import Platform
 from ukis_pysat.raster import Image
@@ -30,6 +31,45 @@ class DataTest(unittest.TestCase):
 
         self.assertTrue(np.array_equal(img.arr, Image(dataset=img.dataset, arr=img.arr).arr))
 
+    def test_arr(self):
+        img_first = Image(path=os.path.join(os.path.dirname(__file__), "testfiles", "dummy.tif"), dimorder="first")
+        img_first.mask_image(box(11.9027457562112939, 51.4664152338322580, 11.9477435281016131, 51.5009522690838750,))
+        self.assertEqual(img_first.arr.shape, (1, 385, 502))
+        self.assertEqual(
+            str(img_first.transform),
+            str(
+                from_bounds(
+                    11.9027457562112939,
+                    51.4664152338322580,
+                    11.9477435281016131,
+                    51.5009522690838750,
+                    img_first.arr.shape[2],
+                    img_first.arr.shape[1],
+                )
+            ),
+        )
+
+        img_last = Image(path=os.path.join(os.path.dirname(__file__), "testfiles", "dummy.tif"), dimorder="last")
+        img_last.mask_image(box(11.9027457562112939, 51.4664152338322580, 11.9477435281016131, 51.5009522690838750,))
+        self.assertEqual(img_last.arr.shape, (385, 502, 1))
+        self.assertEqual(
+            str(img_last.transform),
+            str(
+                from_bounds(
+                    11.9027457562112939,
+                    51.4664152338322580,
+                    11.9477435281016131,
+                    51.5009522690838750,
+                    img_last.arr.shape[1],
+                    img_last.arr.shape[0],
+                )
+            ),
+        )
+
+        with self.assertRaises(AttributeError, msg="dimorder for bands or channels must be either 'first' or 'last'."):
+            img_first = Image(path=os.path.join(os.path.dirname(__file__), "testfiles", "dummy.tif"), dimorder="middle")
+            img_first.arr
+
     def test_get_valid_data_bbox(self):
         self.assertEqual(
             img.get_valid_data_bbox(), (11.896863892, 51.515176657, 11.896863892, 51.515176657),
@@ -100,7 +140,9 @@ class DataTest(unittest.TestCase):
             {
                 "platform": Platform.Sentinel2,
                 "dn_file": os.path.join(target_dir, "S2B_MSIL1C_20200406T101559_N0209_R065_T32UPC_20200406T130159.tif"),
-                "toa_file": os.path.join(target_dir, "S2B_MSIL1C_20200406T101559_N0209_R065_T32UPC_20200406T130159_toa.tif"),
+                "toa_file": os.path.join(
+                    target_dir, "S2B_MSIL1C_20200406T101559_N0209_R065_T32UPC_20200406T130159_toa.tif"
+                ),
                 "mtl_file": None,
                 "wavelengths": None,
             },
@@ -109,7 +151,9 @@ class DataTest(unittest.TestCase):
         for i in range(len(tests)):
             img_dn = Image(path=tests[i]["dn_file"])
             img_toa = Image(path=tests[i]["toa_file"])
-            img_dn.dn2toa(platform=tests[i]["platform"], mtl_file=tests[i]["mtl_file"], wavelengths=tests[i]["wavelengths"])
+            img_dn.dn2toa(
+                platform=tests[i]["platform"], mtl_file=tests[i]["mtl_file"], wavelengths=tests[i]["wavelengths"]
+            )
             self.assertTrue(np.array_equal(img_dn.arr, img_toa.arr))
 
     def test__lookup_bands(self):