Point requests for place geoblock (#51)

CHANGES.rst

@@ -15,6 +15,11 @@ Changelog of dask-geomodeling
 - Added the 'statistics' argument to raster.spatial.Place to deal with
   overlapping features. 
 
+- Allow point requests in raster.spatial.Place.
+
+- Clarifications about raster cell validity ranges in MemorySource and
+  RasterFileSource.
+
 
 2.2.7 (2020-04-30)
 ------------------

dask_geomodeling/raster/sources.py

@@ -20,6 +20,10 @@ class MemorySource(RasterBlock):
     Nodata values are supported, but when upsampling the data, these are
     assumed to be 0 biasing data edges towards 0.
 
+    The raster pixel with its topleft corner at [x, y] will define ranges
+    [x, x + dx) and (y - dy, y]. Here [dx, dy] denotes the (unsigned) pixel
+    size. The topleft corner and top and left edges belong to a pixel.
+
     :param data: the pixel values
         this value will be transformed in a 3D array (t, y, x)
     :param no_data_value: the pixel value that designates 'no data'
@@ -282,6 +286,11 @@ def process(process_kwargs):
 class RasterFileSource(RasterBlock):
     """A raster source that interfaces data from a file path.
 
+    The value at raster cell with its topleft corner at [x, y] is assumed to
+    define a value for ranges [x, x + dx) and (y - dy, y]. Here [dx, dy]
+    denotes the (unsigned) pixel size. The topleft corner and top and
+    left edges belong to a pixel.
+
     :param url: the path to the file. File paths have to be contained inside
       the current root setting. Relative paths are interpreted relative to this
       setting (but internally stored as absolute paths).

dask_geomodeling/raster/spatial.py

@@ -583,29 +583,31 @@ def get_sources_and_requests(self, **request):
         size_x = (x2 - x1) / request["width"]
         size_y = (y2 - y1) / request["height"]
 
-        # check what the full source extent would require
-        full_height = math.ceil((ymax - ymin) / size_y)
-        full_width = math.ceil((xmax - xmin) / size_x)
-        if full_height * full_width <= request["width"] * request["height"]:
-            _request = request.copy()
-            _request["width"] = full_width
-            _request["height"] = full_height
-            _request["bbox"] = (
-                xmin,
-                ymin,
-                xmin + full_width * size_x,
-                ymin + full_height * size_y,
-            )
-            process_kwargs = {
-                "mode": "warp",
-                "anchor": anchor,
-                "coordinates": coordinates,
-                "src_bbox": _request["bbox"],
-                "dst_bbox": request["bbox"],
-                "cellsize": (size_x, size_y),
-                "statistic": self.statistic,
-            }
-            return [(process_kwargs, None), (self.store, _request)]
+        # point requests: never request the full source extent
+        if size_x > 0 and size_y > 0:
+            # check what the full source extent would require
+            full_height = math.ceil((ymax - ymin) / size_y)
+            full_width = math.ceil((xmax - xmin) / size_x)
+            if full_height * full_width <= request["width"] * request["height"]:
+                _request = request.copy()
+                _request["width"] = full_width
+                _request["height"] = full_height
+                _request["bbox"] = (
+                    xmin,
+                    ymin,
+                    xmin + full_width * size_x,
+                    ymin + full_height * size_y,
+                )
+                process_kwargs = {
+                    "mode": "warp",
+                    "anchor": anchor,
+                    "coordinates": coordinates,
+                    "src_bbox": _request["bbox"],
+                    "dst_bbox": request["bbox"],
+                    "cellsize": (size_x, size_y),
+                    "statistic": self.statistic,
+                }
+                return [(process_kwargs, None), (self.store, _request)]
 
         # generate a new (backwards shifted) bbox for each coordinate
         sources_and_requests = []
@@ -618,7 +620,9 @@ def get_sources_and_requests(self, **request):
                 y2 + anchor[1] - _y,
             ]
             # check the overlap with the source's extent
-            if bbox[0] >= xmax or bbox[1] >= ymax or bbox[2] <= xmin or bbox[3] <= ymin:
+            # Note that raster cells are defined [xmin, xmax) and (ymin, ymax]
+            # so points precisely at xmax or ymin certainly do not have data.
+            if bbox[0] >= xmax or bbox[1] > ymax or bbox[2] < xmin or bbox[3] <= ymin:
                 continue
             filtered_coordinates.append((_x, _y))
             _request = request.copy()

dask_geomodeling/tests/test_raster_sources.py

@@ -52,6 +52,7 @@ def test_geometry(self):
         self.assertEqual(expected, self.source.geometry.ExportToWkt())
 
     def test_point_single_pixel(self):
+        # data is defined at [136700, 136705) and (455795, 455800]
         for dx, dy in ((0, 0), (0, -4.99), (4.99, 0), (4.99, -4.99)):
             data = self.source.get_data(
                 mode="vals",
@@ -64,7 +65,10 @@ def test_point_single_pixel(self):
             assert_equal(data["values"], 5)
 
     def test_point_single_pixel_nodata(self):
-        for dx, dy in ((-0.01, -1), (1, 0.01), (5.01, -1), (1, -5.01)):
+        # data is defined at [136700, 136705) and (455795, 455800]
+        for dx, dy in (
+            (0, -5.0), (5.0, 0), (-5.0, 5.0), (-0.01, 0), (0, 0.01),
+        ):
             data = self.source.get_data(
                 mode="vals",
                 projection="EPSG:28992",

dask_geomodeling/tests/test_raster_spatial.py

@@ -205,3 +205,44 @@ def test_place_time_request(source, center):
 def test_place_meta_request(source, center):
     place = raster.Place(source, "EPSG:28992", center, [(150, 50)])
     assert source.get_data(mode="meta") == place.get_data(mode="meta")
+
+
+@pytest.mark.parametrize("point,expected", [
+    ((5, 15), 7),  # zone 1
+    ((15, 15), 255),  # zone 2
+    ((5, 5), 255),  # zone 3
+    ((15, 5), 7),  # zone 4
+    ((10, 15), 255),  # line 1-2
+    ((5, 10), 255),  # line 1-3
+    ((15, 10), 7),  # line 2-4
+    ((10, 5), 7),  # line 3-4
+    ((10, 10), 7),  # center
+])
+def test_place_point_request(source, center, point, expected):
+    # For point requests, edges are important. Let's do a drawing:
+    # 20  _______
+    #    |       |
+    #    |   1   |   2
+    # 10 |_______|_______
+    #            |       |
+    #        3   |   4   |
+    #  0         |_______|
+    #    0       10     20
+    # - zone 1 and 4 are filled; zone 2 and 3 are empty (see below coordinates)
+    # A pixel includes its topleft corner and top and left edges
+    # - line between 2-4 and 3-4 are filled; line 1-2 and 1-3 are empty
+    # - center point at (10, 10) is filled
+    coordinates = [(60, -40), (-40, 60)]
+    place = raster.Place(
+        source, "EPSG:28992", anchor=center, coordinates=coordinates
+    )
+    point_request = dict(
+        mode="vals",
+        bbox=point * 2,
+        projection="EPSG:28992",
+        width=1,
+        height=1,
+    )
+    values = place.get_data(**point_request)["values"]
+    assert values.shape == (1, 1, 1)
+    assert values.item() == expected