Repository: dlr-eoc/ukis-pysat

Issue:
__update_dataset() in raster.Image() does not update dtype
[__update_dataset()](https://github.com/dlr-eoc/ukis-pysat/blob/b249ee8210ce43034889ef0452062d117513cc92/ukis_pysat/raster.py#L153) does not update dtype. The following is a simple example that shows the problem:

```python
from ukis_pysat.raster import Image
from ukis_pysat.members import Platform


file = "D:/tmp/test/T32UPU_20200413T100549_B01.jp2"

# read image
img = Image(data=file, dimorder="last")
print(img.arr.dtype)
>> uint16

# dn2toa
# NOTE: this modifies img dtype
img.dn2toa(platform=Platform.Sentinel2)
print(img.arr.dtype)
>> float32

# subset
# BUG: this overwrite the img dtype modified by dn2toa
img.mask_image(bbox=(654700.76, 663334.23, 5324479.29, 5330987.60), crop=True, pad=False)
print(img.arr.dtype)
>> uint16
```
Suggest to make sure that [__update_dataset()](https://github.com/dlr-eoc/ukis-pysat/blob/b249ee8210ce43034889ef0452062d117513cc92/ukis_pysat/raster.py#L153) upon calling takes what it can get from the array and enforce passing crs and transform information separately (similar to what we do with creating an Image from an array, just in memory).


Discussion:

Patch:
diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index d15c395..57aa9e6 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -23,11 +23,13 @@ Added
 Fixed
 *****
 - ``file.get_ts_from_sentinel_filename()``: Return datetime.datetime objects instead of timestamp strings #42
-- ``raster.Image()``: in-memory dataset could not be updated if not GTiff #48
+- ``raster.Image()``: in-memory dataset could not be updated if not GTiff and other improvements #48 #52
 
 Changed
 *******
+- ``raster.Image()``: renamed `mask_image()` to `mask()`
 - ``raster.Image()``: update of init signature to be less confusing #41 #50
+- ``raster.Image()``: in-memory dataset now always with "GTiff" driver #53
 
 
 [0.4.0]  (2020-06-05)
diff --git a/ukis_pysat/raster.py b/ukis_pysat/raster.py
index e0ac79c..15ab56a 100644
--- a/ukis_pysat/raster.py
+++ b/ukis_pysat/raster.py
@@ -43,14 +43,10 @@ class Image:
 
         if isinstance(data, rasterio.io.DatasetReader):
             self.dataset = data
-            self.crs = self.dataset.crs
-            self.transform = self.dataset.transform
             self.__arr = self.dataset.read()
 
         elif isinstance(data, str):
             self.dataset = rasterio.open(data)
-            self.crs = self.dataset.crs
-            self.transform = self.dataset.transform
             self.__arr = self.dataset.read()
 
         elif isinstance(data, np.ndarray):
@@ -62,11 +58,7 @@ class Image:
                 self.__arr = data
             else:
                 self.__arr = reshape_as_raster(data)
-            meta = {"dtype": self.__arr.dtype, "count": self.__arr.shape[0], "crs": crs, "driver": "GTiff"}
-            self.transform = transform
-            self.dataset = self.__update_dataset(meta).open()
-            self.crs = self.dataset.crs
-            self.transform = self.dataset.transform
+            self.dataset = self.__update_dataset(crs, transform)
         else:
             raise TypeError("dataset must be of type rasterio.io.DatasetReader, str or np.ndarray")
 
@@ -85,9 +77,15 @@ class Image:
         :return: tuple with valid data bounds
         """
         valid_data_window = windows.get_data_window(self.__arr, nodata=nodata)
-        return windows.bounds(valid_data_window, windows.transform(valid_data_window, self.transform))
+        return windows.bounds(valid_data_window, windows.transform(valid_data_window, self.dataset.transform))
 
     def mask_image(self, bbox, crop=True, pad=False, mode="constant", constant_values=0):
+        from warnings import warn
+
+        warn("`mask_image()` was renamed to `mask()` and will be removed with version 0.6.0", DeprecationWarning)
+        return self.mask(bbox, crop=crop, pad=pad, mode=mode, constant_values=constant_values)
+
+    def mask(self, bbox, crop=True, pad=False, mode="constant", constant_values=0):
         """Mask the area outside of the input shapes with no data.
 
         :param bbox: bounding box of type tuple or Shapely Polygon
@@ -100,19 +98,19 @@ class Image:
         if pad:
             self.dataset.close()
             # here we __update_dataset()
-            self.dataset = self._pad_to_bbox(bbox, mode, constant_values).open(**self.dataset.meta)
+            self.dataset = self._pad_to_bbox(bbox, mode, constant_values)
 
         if isinstance(bbox, polygon.Polygon):
-            self.__arr, self.transform = rasterio.mask.mask(self.dataset, [bbox], crop=crop)
+            self.__arr, 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, transform = rasterio.mask.mask(self.dataset, [box(*bbox)], crop=crop)
         else:
             raise TypeError(f"bbox must be of type tuple or Shapely Polygon")
 
         # update for further processing
         self.dataset.close()  # meta still available when DataSetReader is closed
         # passing meta to open, because driver gets lost if it's not GTiff
-        self.dataset = self.__update_dataset(self.dataset.meta).open(**self.dataset.meta)
+        self.dataset = self.__update_dataset(self.dataset.crs, transform, nodata=self.dataset.nodata)
 
     def _pad_to_bbox(self, bbox, mode="constant", constant_values=0):
         """Buffers array with biggest difference to bbox and adjusts affine transform matrix. Can be used to fill
@@ -134,7 +132,7 @@ class Image:
         max_diff_ll = np.max(np.subtract(tuple(self.dataset.bounds[:2]), bbox[:2]))
         max_diff = max(max_diff_ll, max_diff_ur)  # buffer in units
 
-        pad_width = math.ceil(max_diff / self.transform.to_gdal()[1])  # units / pixel_size
+        pad_width = math.ceil(max_diff / self.dataset.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,),
@@ -142,32 +140,42 @@ class Image:
         )
 
         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,
+            destination[i], transform = rasterio.pad(
+                self.__arr[0], self.dataset.transform, pad_width, mode, constant_values=constant_values,
             )
 
         self.__arr = destination
 
-        return self.__update_dataset(self.dataset.meta)
+        self.dataset.close()
+        return self.__update_dataset(self.dataset.crs, transform, nodata=self.dataset.nodata)
 
-    def __update_dataset(self, meta):
+    def __update_dataset(self, crs, transform, nodata=None):
         """Update dataset without writing to file after it theoretically changed.
 
-        :param meta: The basic metadata of the dataset as returned from the meta property of rasterio Datasets
-        :return: file in memory
+        :param crs: crs of the dataset
+        :param transform: transform of the dataset
+        :param nodata: nodata value, optional
+        :return: file in memory, open as dataset
         """
 
-        meta.update(
-            {"height": self.__arr.shape[-2], "width": self.__arr.shape[-1], "transform": self.transform,}
-        )
+        meta = {
+            "driver": "GTiff",
+            "dtype": self.__arr.dtype,
+            "nodata": nodata,
+            "height": self.__arr.shape[-2],
+            "width": self.__arr.shape[-1],
+            "count": self.__arr.shape[0],
+            "crs": crs,
+            "transform": transform,
+        }
 
         memfile = MemoryFile()
         with memfile.open(**meta) as ds:
             ds.write(self.__arr)
 
-        return memfile
+        return memfile.open()
 
-    def warp(self, dst_crs, resampling_method=0, num_threads=4, resolution=None):
+    def warp(self, dst_crs, resampling_method=0, num_threads=4, resolution=None, nodata=None):
         """Reproject a source raster to a destination raster.
 
         :param dst_crs: CRS or dict, Target coordinate reference system.
@@ -176,6 +184,7 @@ class Image:
         :param num_threads: int, number of workers, optional (default: 4)
         :param resolution: tuple (x resolution, y resolution) or float, optional.
             Target resolution, in units of target coordinate reference system.
+        :param nodata: nodata value of source, int or float, optional.
         """
         # output dimensions and transform for reprojection.
         if resolution:
@@ -194,26 +203,21 @@ class Image:
 
         destination = np.zeros((self.dataset.count, height, width), self.__arr.dtype)
 
-        for i in range(0, self.dataset.count):
-            reproject(
-                source=rasterio.band(self.dataset, i + 1),  # index starting at 1
-                destination=destination[i],
-                src_transform=self.dataset.transform,
-                src_crs=self.dataset.crs,
-                dst_transform=transform,
-                dst_crs=dst_crs,
-                resampling=resampling_method,
-                num_threads=num_threads,
-            )
-
-        # update for further processing
-        self.__arr = destination
-        self.transform = transform
-        self.crs = dst_crs
+        self.__arr, transform = reproject(
+            source=self.__arr,
+            destination=destination,
+            src_transform=self.dataset.transform,
+            src_crs=self.dataset.crs,
+            src_nodata=nodata,
+            dst_transform=transform,
+            dst_crs=dst_crs,
+            dst_nodata=nodata,
+            resampling=resampling_method,
+            num_threads=num_threads,
+        )
 
         self.dataset.close()
-        # passing meta to open, because driver gets lost if it's not GTiff
-        self.dataset = self.__update_dataset(self.dataset.meta).open(**self.dataset.meta)
+        self.dataset = self.__update_dataset(dst_crs, transform, nodata=nodata)
 
     def dn2toa(self, platform, mtl_file=None, wavelengths=None):
         """This method converts digital numbers to top of atmosphere reflectance, like described here:
@@ -284,6 +288,9 @@ class Image:
                 f"Sentinel-2]. "
             )
 
+        self.dataset.close()
+        self.dataset = self.__update_dataset(self.dataset.crs, self.dataset.transform, nodata=self.dataset.nodata)
+
     @staticmethod
     def _lookup_bands(platform, wavelengths):
         """lookup for bands and their according wavelength for Landsat-5, -7 & -8.
@@ -403,8 +410,8 @@ class Image:
                 "height": self.__arr.shape[-2],
                 "width": self.__arr.shape[-1],
                 "dtype": dtype,
-                "transform": self.transform,
-                "crs": self.crs,
+                "transform": self.dataset.transform,
+                "crs": self.dataset.crs,
             }
         )
 


Test patch:
diff --git a/tests/test_raster.py b/tests/test_raster.py
index b09158a..dc4c233 100644
--- a/tests/test_raster.py
+++ b/tests/test_raster.py
@@ -40,16 +40,16 @@ class DataTest(unittest.TestCase):
 
     def test_init_with_arry_fail_missing_transform(self):
         with self.assertRaises(TypeError):
-            Image(data=self.img.arr, crs=self.img.crs)
+            Image(data=self.img.arr, crs=self.img.dataset.crs)
 
     def test_init_with_array(self):
-        img = Image(self.img.arr, crs=self.img.crs, transform=self.img.transform)
+        img = Image(self.img.arr, crs=self.img.dataset.crs, transform=self.img.dataset.transform)
         self.assertTrue(np.array_equal(self.img.arr, img.arr))
         img.close()
 
     def test_init_with_arry_fail_missing_crs(self):
         with self.assertRaises(TypeError):
-            Image(data=self.img.arr, crs=self.img.transform)
+            Image(data=self.img.arr, crs=self.img.dataset.transform)
 
     def test_dimorder_error(self):
         with self.assertRaises(TypeError):
@@ -57,10 +57,10 @@ class DataTest(unittest.TestCase):
 
     def test_arr(self):
         img_first = Image(TEST_FILE, dimorder="first")
-        img_first.mask_image(box(11.9027457562112939, 51.4664152338322580, 11.9477435281016131, 51.5009522690838750,))
+        img_first.mask(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(img_first.dataset.transform),
             str(
                 from_bounds(
                     11.9027457562112939,
@@ -75,10 +75,10 @@ class DataTest(unittest.TestCase):
         img_first.close()
 
         img_last = Image(TEST_FILE, dimorder="last")
-        img_last.mask_image(box(11.9027457562112939, 51.4664152338322580, 11.9477435281016131, 51.5009522690838750,))
+        img_last.mask(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(img_last.dataset.transform),
             str(
                 from_bounds(
                     11.9027457562112939,
@@ -92,11 +92,15 @@ class DataTest(unittest.TestCase):
         )
         img_last.close()
 
-        img_first = Image(np.ones((1, 385, 502)), dimorder="first", crs=self.img.crs, transform=self.img.transform)
+        img_first = Image(
+            np.ones((1, 385, 502)), dimorder="first", crs=self.img.dataset.crs, transform=self.img.dataset.transform
+        )
         self.assertEqual(img_first.arr.shape, (1, 385, 502))
         img_first.close()
 
-        img_last = Image(np.ones((385, 502, 1)), dimorder="last", crs=self.img.crs, transform=self.img.transform)
+        img_last = Image(
+            np.ones((385, 502, 1)), dimorder="last", crs=self.img.dataset.crs, transform=self.img.dataset.transform
+        )
         self.assertEqual(img_last.arr.shape, (385, 502, 1))
         img_last.close()
 
@@ -110,9 +114,9 @@ class DataTest(unittest.TestCase):
 
     def test_mask_image(self):
         with self.assertRaises(TypeError, msg="bbox must be of type tuple or Shapely Polygon"):
-            self.img.mask_image([1, 2, 3])
+            self.img.mask([1, 2, 3])
 
-        self.img.mask_image(box(11.9027457562112939, 51.4664152338322580, 11.9477435281016131, 51.5009522690838750,))
+        self.img.mask(box(11.9027457562112939, 51.4664152338322580, 11.9477435281016131, 51.5009522690838750, ))
         self.assertEqual(
             self.img.dataset.bounds,
             BoundingBox(
@@ -120,7 +124,7 @@ class DataTest(unittest.TestCase):
             ),
         )
 
-        self.img.mask_image((11.9027457562112939, 51.4664152338322580, 11.9477435281016131, 51.5009522690838750,))
+        self.img.mask((11.9027457562112939, 51.4664152338322580, 11.9477435281016131, 51.5009522690838750,))
         self.assertEqual(
             self.img.dataset.bounds,
             BoundingBox(
@@ -128,7 +132,7 @@ class DataTest(unittest.TestCase):
             ),
         )
 
-        self.img.mask_image(
+        self.img.mask(
             box(11.8919236802142620, 51.4664152338322580, 11.9477435281016131, 51.5009522690838750,), pad=True,
         )
         self.assertEqual(
@@ -139,13 +143,14 @@ class DataTest(unittest.TestCase):
         )
 
     def test_warp(self):
-        self.assertEqual(self.img.crs, "EPSG:4326")
+        self.assertEqual(self.img.dataset.crs, "EPSG:4326")
 
         self.img.warp("EPSG:3857")
-        self.assertEqual(self.img.crs, "EPSG:3857")
+        self.assertEqual(self.img.dataset.crs, "EPSG:3857")
+        self.assertEqual(self.img.dataset.meta["crs"], "EPSG:3857")
 
         self.img.warp("EPSG:4326", resolution=1.0)
-        self.assertEqual(1.0, self.img.transform.to_gdal()[1])
+        self.assertEqual(1.0, self.img.dataset.transform.to_gdal()[1])
 
     def test_dn2toa(self):
         target_dir = os.path.join(os.path.dirname(__file__), "testfiles", "satellite_data")
@@ -242,7 +247,7 @@ class DataTest(unittest.TestCase):
         img2.close()
         os.remove(r"result.tif")
 
-        img3 = Image(self.img.arr, crs=self.img.crs, transform=self.img.transform)
+        img3 = Image(self.img.arr, crs=self.img.dataset.crs, transform=self.img.dataset.transform)
         img3.write_to_file(r"result.tif", np.uint16)
         img3.close()
         img4 = Image("result.tif")