Repository: dlr-eoc/ukis-pysat

Issue:
Target align option for Image.warp
In many applications it is extremely useful to be able to force align two raster grids to match each others crs, extent, dimensions and grid alignment. With a combination of `warp` and `mask` we can come close to align two rasters (this actually works in many cases). It may, however, happen that the raster grids are shifted with subpixel resolution, which prevents them from being exactly on the same grid and therefore they cannot be compared at array level with for example numpy. 

I propose the following solution to add a `target_align` option to `warp`, which aligns one Image instance to another one:

```python
    def warp(self, dst_crs, resampling_method=0, num_threads=4, resolution=None, nodata=None, target_align=None):
        """Reproject a source raster to a destination raster.
        :param dst_crs: CRS or dict, Target coordinate reference system.
        :param resampling_method: Resampling algorithm, int, defaults to 0 (Nearest)
            numbers: https://github.com/mapbox/rasterio/blob/master/rasterio/enums.py#L28
        :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 target_align: raster to which to align resolution, extent and gridspacing, optional (Image).
        :param nodata: nodata value of source, int or float, optional.
        """
        if target_align:
            transform = target_align.dataset.transform
            width = target_align.dataset.width
            height = target_align.dataset.height

        else:
            if resolution:
                transform, width, height = calculate_default_transform(
                    self.dataset.crs,
                    dst_crs,
                    self.dataset.width,
                    self.dataset.height,
                    *self.dataset.bounds,
                    resolution=resolution,
                )

            else:
                transform, width, height = calculate_default_transform(
                    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)

        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()
        self.dataset = self.__update_dataset(dst_crs, transform, nodata=nodata)
```

Discussion:

Patch:
diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index 9422790..e34dfcf 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -19,11 +19,13 @@ Added
 *****
 - ``raster.Image()``: With Statement Context Manager for Image #45
 - ``raster.Image()``: Alter image array #67
+- ``raster.Image()``: Target align option for `warp()` #60
 
 Fixed
 *****
 - ``data.Source()``: Fixed query metadata return for new EarthExplorer API #71
 - ``raster.Image()``: Consider all image bands in `pad()` #59
+- ``raster.Image()``: Memory leak caused by `__update_dataset()` #62
 
 Changed
 *******
diff --git a/requirements.txt b/requirements.txt
index 0939c24..5e42d44 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,7 +1,7 @@
 shapely~=1.7.0
 numpy~=1.18.1
 requests~=2.23.0
-rasterio~=1.1.0
+rasterio>=1.1.0,<1.1.5
 pillow~=7.1.2
 sentinelsat~=0.13
 landsatxplore==0.6
diff --git a/setup.py b/setup.py
index afb5073..a92da5f 100644
--- a/setup.py
+++ b/setup.py
@@ -39,7 +39,7 @@ extras_require = {
         "python-dateutil~=2.8.1",
         "shapely~=1.7.0",
     ],
-    "raster": ["numpy~=1.18.1", "rasterio~=1.1.0", "rio_toa==0.3.0", "shapely~=1.7.0",],
+    "raster": ["numpy~=1.18.1", "rasterio!=1.1.5", "rio_toa==0.3.0", "shapely~=1.7.0",],
 }
 extras_require["complete"] = list(set([v for req in extras_require.values() for v in req]))
 
diff --git a/ukis_pysat/raster.py b/ukis_pysat/raster.py
index 97d0dc7..c4b3176 100644
--- a/ukis_pysat/raster.py
+++ b/ukis_pysat/raster.py
@@ -57,8 +57,9 @@ class Image:
             if dimorder == "first":
                 self.__arr = data
             else:
-                self.__arr = reshape_as_raster(data)
-            self.dataset = self.__update_dataset(crs, transform)
+                self.__arr = reshape_as_raster(data)
+            self.dataset = None
+            self.__update_dataset(crs, transform)
         else:
             raise TypeError("dataset must be of type rasterio.io.DatasetReader, str or np.ndarray")
 
@@ -114,6 +115,7 @@ class Image:
         """
         if pad:
             from warnings import warn
+
             warn("`pad` was renamed to `fill` in `mask()` and will be removed with version 0.7.0", DeprecationWarning)
             fill = pad
 
@@ -122,7 +124,7 @@ class Image:
             pad_width = self._get_pad_width(bbox)
             if pad_width > 0:
                 # only pad raster if it is smaller than bbox
-                self.dataset = self.pad(pad_width, mode, constant_values)
+                self.pad(pad_width, mode, constant_values)
 
         if isinstance(bbox, polygon.Polygon):
             self.__arr, transform = rasterio.mask.mask(self.dataset, [bbox], crop=crop)
@@ -131,10 +133,7 @@ class Image:
         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.crs, transform, nodata=self.dataset.nodata)
+        self.__update_dataset(self.dataset.crs, transform, nodata=self.dataset.nodata)
 
     def _get_pad_width(self, bbox):
         """Calculates biggest difference from raster bounds to bbox. Can be used with pad() to fill array with
@@ -175,9 +174,7 @@ class Image:
             )
 
         self.__arr = destination
-
-        self.dataset.close()
-        return self.__update_dataset(self.dataset.crs, transform, nodata=self.dataset.nodata)
+        self.__update_dataset(self.dataset.crs, transform, nodata=self.dataset.nodata)
 
     def __update_dataset(self, crs, transform, nodata=None):
         """Update dataset without writing to file after it theoretically changed.
@@ -197,15 +194,15 @@ class Image:
             "count": self.__arr.shape[0],
             "crs": crs,
             "transform": transform,
-        }
-
-        memfile = MemoryFile()
-        with memfile.open(**meta) as ds:
-            ds.write(self.__arr)
-
-        return memfile.open()
-
-    def warp(self, dst_crs, resampling_method=0, num_threads=4, resolution=None, nodata=None):
+        }
+
+        memfile = MemoryFile()
+        with memfile.open(**meta) as ds:
+            ds.write(self.__arr)
+        self.dataset = memfile.open()
+        memfile.close()
+
+    def warp(self, dst_crs, resampling_method=0, num_threads=4, resolution=None, nodata=None, target_align=None):
         """Reproject a source raster to a destination raster.
 
         :param dst_crs: CRS or dict, Target coordinate reference system.
@@ -214,22 +211,28 @@ 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 target_align: raster to which to align resolution, extent and gridspacing, optional (Image).
         :param nodata: nodata value of source, int or float, optional.
-        """
-        # output dimensions and transform for reprojection.
-        if resolution:
-            transform, width, height = calculate_default_transform(
-                self.dataset.crs,
-                dst_crs,
-                self.dataset.width,
-                self.dataset.height,
-                *self.dataset.bounds,
-                resolution=resolution,
-            )
+        """
+        if target_align:
+            transform = target_align.dataset.transform
+            width = target_align.dataset.width
+            height = target_align.dataset.height
+
         else:
-            transform, width, height = calculate_default_transform(
-                self.dataset.crs, dst_crs, self.dataset.width, self.dataset.height, *self.dataset.bounds,
-            )
+            if resolution:
+                transform, width, height = calculate_default_transform(
+                    self.dataset.crs,
+                    dst_crs,
+                    self.dataset.width,
+                    self.dataset.height,
+                    *self.dataset.bounds,
+                    resolution=resolution,
+                )
+            else:
+                transform, width, height = calculate_default_transform(
+                    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)
 
@@ -246,8 +249,7 @@ class Image:
             num_threads=num_threads,
         )
 
-        self.dataset.close()
-        self.dataset = self.__update_dataset(dst_crs, transform, nodata=nodata)
+        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:
@@ -317,8 +319,7 @@ class Image:
                 f"Sentinel-2]. "
             )
 
-        self.dataset.close()
-        self.dataset = self.__update_dataset(self.dataset.crs, self.dataset.transform, nodata=self.dataset.nodata)
+        self.__update_dataset(self.dataset.crs, self.dataset.transform, nodata=self.dataset.nodata)
 
     @staticmethod
     def _lookup_bands(platform, wavelengths):


Test patch:
diff --git a/tests/test_raster.py b/tests/test_raster.py
index 0bf67ca..1f65f72 100644
--- a/tests/test_raster.py
+++ b/tests/test_raster.py
@@ -175,6 +175,14 @@ class DataTest(unittest.TestCase):
         self.img.warp("EPSG:4326", resolution=1.0)
         self.assertEqual(1.0, self.img.dataset.transform.to_gdal()[1])
 
+        source_img = Image(TEST_FILE)
+        source_img.warp("EPSG:3857", resolution=10)
+        target_img = Image(TEST_FILE)
+        target_img.warp("EPSG:3857", resolution=25)
+        self.assertNotEqual(source_img.dataset.transform, target_img.dataset.transform)
+        source_img.warp("EPSG:3857", target_align=target_img)
+        self.assertEqual(source_img.dataset.transform, target_img.dataset.transform)
+
     def test_dn2toa(self):
         target_dir = os.path.join(os.path.dirname(__file__), "testfiles", "satellite_data")
         tests = [