Merge pull request #7 from nens/casper-clip

Fix propagation of 'extent' through raster.Clip

CHANGES.rst

@@ -5,6 +5,12 @@ Changelog of dask-geomodeling
 2.0.4 (unreleased)
 ------------------
 
+- Fixed propagation of the 'extent' and 'geometry' attributes through the
+  raster.Clip. Both now return the intersection of the store and mask rasters.
+
+- The MemorySource and elementwise blocks now return None for 'extent' and
+  'geometry' if they are empty.
+
 - Preserve functionality of the geometry.Difference block with geopandas 0.6.
   When taking the difference of a geometry with a missing geometry (A - None),
   geopandas < 0.6 returned A as result, while >= 0.6 returns None as result.

dask_geomodeling/raster/elemwise.py

@@ -159,6 +159,8 @@ def geometry(self):
                 geometry = geometry.Clone()
                 geometry.TransformTo(sr)
             result = result.Intersection(geometry)
+        if result.GetArea() == 0.0:
+            return
         return result
 
     @property

dask_geomodeling/raster/misc.py

@@ -68,9 +68,37 @@ def process(data, source_data):
         values[mask] = data["no_data_value"]
         return {"values": values, "no_data_value": data["no_data_value"]}
 
+    @property
+    def extent(self):
+        """Intersection of bounding boxes of 'store' and 'source'. """
+        result, mask = [s.extent for s in self.args]
+        if result is None or mask is None:
+            return
+
+        # return the overlapping box
+        x1 = max(result[0], mask[0])
+        y1 = max(result[1], mask[1])
+        x2 = min(result[2], mask[2])
+        y2 = min(result[3], mask[3])
+        if x2 <= x1 or y2 <= y1:
+            return None  # no overlap
+        else:
+            return x1, y1, x2, y2
+
     @property
     def geometry(self):
-        return self.source.geometry
+        """Intersection of geometries of 'store' and 'source'. """
+        result, mask = [x.geometry for x in self.args]
+        if result is None or mask is None:
+            return
+        sr = result.GetSpatialReference()
+        if not mask.GetSpatialReference().IsSame(sr):
+            mask = mask.Clone()
+            mask.TransformTo(sr)
+        result = result.Intersection(mask)
+        if result.GetArea() == 0.0:
+            return
+        return result
 
 
 class Mask(BaseSingle):

dask_geomodeling/raster/sources.py

@@ -140,17 +140,24 @@ def geo_transform(self):
         return utils.GeoTransform((p, a, 0, q, 0, -d))
 
     def _get_extent(self):
+        if not self.data.size:
+            return
         bbox = self.geo_transform.get_bbox((0, 0), self.data.shape[1:])
         return utils.Extent(bbox, utils.get_sr(self.projection))
 
     @property
     def extent(self):
-        extent_epsg4326 = self._get_extent().transformed(utils.EPSG4326)
-        return extent_epsg4326.bbox
+        extent = self._get_extent()
+        if extent is None:
+            return
+        return extent.transformed(utils.EPSG4326).bbox
 
     @property
     def geometry(self):
-        return self._get_extent().as_geometry()
+        extent = self._get_extent()
+        if extent is None:
+            return
+        return extent.as_geometry()
 
     def __len__(self):
         return self.data.shape[0]

dask_geomodeling/tests/test_raster.py

@@ -1506,38 +1506,6 @@ def test_shift(self):
         )
         self.assertEqual(view2.time, view.time)
 
-    def test_clip(self):
-        # store returns None
-        clip = raster.Clip(store=self.raster_none, source=self.raster_none)
-        self.assertIsNone(clip.get_data(**self.vals_request))
-
-        # store returns no data
-        clip = raster.Clip(store=self.raster_nodata, source=self.raster_none)
-        assert_equal(clip.get_data(**self.vals_request)["values"], 255)
-
-        # source returns None
-        clip = raster.Clip(store=self.raster, source=self.raster_none)
-        self.assertIsNone(clip.get_data(**self.vals_request))
-
-        # source has nodata everywhere (everything will be masked)
-        clip = raster.Clip(store=self.raster, source=self.raster_nodata)
-        assert_equal(clip.get_data(**self.vals_request)["values"], 255)
-
-        # source has data everywhere (nothing will be masked)
-        clip = raster.Clip(store=self.raster, source=self.raster)
-        assert_equal(clip.get_data(**self.vals_request)["values"], 7)
-
-        self.assertEqual(clip.get_data(**self.meta_request)["meta"], self.expected_meta)
-        self.assertEqual(clip.get_data(**self.time_request)["time"], self.expected_time)
-
-        # Use a boolean mask with True everywhere, should propagate data
-        clip = raster.Clip(self.raster, self.raster == 7)
-        assert_equal(clip.get_data(**self.vals_request)["values"], 7)
-
-        # Use a boolean mask with False everywhere, should propagate nodata
-        clip = raster.Clip(self.raster, self.raster != 7)
-        assert_equal(clip.get_data(**self.vals_request)["values"], 255)
-
     def test_mask(self):
         view = raster.Mask(store=self.raster, value=8)
         data = view.get_data(**self.vals_request)

dask_geomodeling/tests/test_raster_misc.py

@@ -29,6 +29,32 @@ def source():
     )
 
 
+@pytest.fixture
+def empty_source():
+    yield MemorySource(
+        data=np.empty((0, 0, 0), dtype=np.uint8),
+        no_data_value=255,
+        projection="EPSG:28992",
+        pixel_size=0.5,
+        pixel_origin=(135000, 456000),
+    )
+
+
+@pytest.fixture
+def nodata_source():
+    time_first = datetime(2000, 1, 1)
+    time_delta = timedelta(hours=1)
+    yield MemorySource(
+        data=np.full((3, 10, 10), 255, dtype=np.uint8),
+        no_data_value=255,
+        projection="EPSG:28992",
+        pixel_size=0.5,
+        pixel_origin=(135000, 456000),
+        time_first=time_first,
+        time_delta=time_delta,
+    )
+
+
 @pytest.fixture
 def vals_request():
     bands = 3
@@ -45,6 +71,22 @@ def vals_request():
     }
 
 
+@pytest.fixture
+def vals_request_none():
+    bands = 3
+    time_first = datetime(2001, 1, 1)
+    time_delta = timedelta(hours=1)
+    yield {
+        "mode": "vals",
+        "start": time_first,
+        "stop": time_first + bands * time_delta,
+        "width": 4,
+        "height": 6,
+        "bbox": (135000, 456000 - 3, 135000 + 2, 456000),
+        "projection": "EPSG:28992",
+    }
+
+
 @pytest.fixture
 def expected_meta():
     bands = 3
@@ -59,6 +101,114 @@ def expected_time():
     return [time_first + i * time_delta for i in range(bands)]
 
 
+def test_clip_attrs_store_empty(source, empty_source):
+    # clip should propagate the (empty) extent of the store
+    clip = raster.Clip(empty_source, source)
+    assert clip.extent is None
+    assert clip.geometry is None
+
+
+def test_clip_attrs_mask_empty(source, empty_source):
+    # clip should propagate the (empty) extent of the clipping mask
+    clip = raster.Clip(source, empty_source)
+    assert clip.extent is None
+    assert clip.geometry is None
+
+
+def test_clip_attrs_intersects(source, empty_source):
+    # create a raster in that only partially overlaps the store
+    clipping_mask = MemorySource(
+        data=source.data,
+        no_data_value=source.no_data_value,
+        projection="EPSG:28992",
+        pixel_size=source.pixel_size,
+        pixel_origin=[o + 3 for o in source.pixel_origin],
+        time_first=source.time_first,
+        time_delta=source.time_delta,
+    )
+    clip = raster.Clip(source, clipping_mask)
+    expected_extent = (
+        clipping_mask.extent[0],
+        clipping_mask.extent[1],
+        source.extent[2],
+        source.extent[3],
+    )
+    expected_geometry = source.geometry.Intersection(clipping_mask.geometry)
+    assert clip.extent == expected_extent
+    assert clip.geometry.ExportToWkt() == expected_geometry.ExportToWkt()
+
+
+def test_clip_attrs_with_reprojection(source, empty_source):
+    # create a raster in WGS84 that contains the store
+    clipping_mask = MemorySource(
+        data=source.data,
+        no_data_value=source.no_data_value,
+        projection="EPSG:4326",
+        pixel_size=1,
+        pixel_origin=(4, 54),
+        time_first=source.time_first,
+        time_delta=source.time_delta,
+    )
+    clip = raster.Clip(source, clipping_mask)
+    assert clip.extent == source.extent
+    assert clip.geometry.GetEnvelope() == source.geometry.GetEnvelope()
+
+
+def test_clip_attrs_no_intersection(source, empty_source):
+    # create a raster in that does not overlap the store
+    clipping_mask = MemorySource(
+        data=source.data,
+        no_data_value=source.no_data_value,
+        projection="EPSG:28992",
+        pixel_size=source.pixel_size,
+        pixel_origin=[o + 5 for o in source.pixel_origin],
+        time_first=source.time_first,
+        time_delta=source.time_delta,
+    )
+    clip = raster.Clip(source, clipping_mask)
+    assert clip.extent is None
+    assert clip.geometry is None
+
+
+def test_clip_empty_source(source, empty_source, vals_request):
+    clip = raster.Clip(empty_source, source)
+    assert clip.get_data(**vals_request) is None
+
+
+def test_clip_with_empty_mask(source, empty_source, vals_request):
+    clip = raster.Clip(source, empty_source)
+    assert clip.get_data(**vals_request) is None
+
+
+def test_clip_with_nodata(source, nodata_source, vals_request):
+    # the clipping mask has nodata everywhere (everything will be masked)
+    clip = raster.Clip(source, nodata_source)
+    assert_equal(clip.get_data(**vals_request)["values"], 255)
+
+
+def test_clip_with_data(source, nodata_source, vals_request):
+    # the clipping mask has data everywhere (nothing will be masked)
+    clip = raster.Clip(source, source)
+    assert_equal(clip.get_data(**vals_request)["values"][:, 0, 0], [1, 7, 255])
+
+
+def test_clip_with_bool(source, vals_request):
+    clip = raster.Clip(source, source == 7)
+    assert_equal(clip.get_data(**vals_request)["values"][:, 0, 0], [255, 7, 255])
+
+
+def test_clip_meta_request(source, vals_request, expected_meta):
+    clip = raster.Clip(source, source)
+    vals_request["mode"] = "meta"
+    assert clip.get_data(**vals_request)["meta"] == expected_meta
+
+
+def test_clip_time_request(source, vals_request, expected_time):
+    clip = raster.Clip(source, source)
+    vals_request["mode"] = "time"
+    assert clip.get_data(**vals_request)["time"] == expected_time
+
+
 def test_reclassify(source, vals_request):
     view = raster.Reclassify(store=source, data=[[7, 1000]])
     data = view.get_data(**vals_request)