Fix point requests for smooth geoblock (#55)

.travis.yml

@@ -7,8 +7,8 @@ matrix:
   - os: linux
     env: PYTHON_VERSION="3.6" MAINDEPS="numpy=1.14 gdal=2.2.4 scipy=1.1 pytz" DEPS="dask-core=0.20 toolz=0.10 pandas=0.23 geopandas=0.5"
   - os: linux
-    env: PYTHON_VERSION="3.7" MAINDEPS="numpy scipy pytz" DEPS="dask-core gdal=2.* toolz geopandas pyproj>=2"
-  - os: osx
+    env: PYTHON_VERSION="3.7" MAINDEPS="numpy=1.16 scipy=1.3 pytz" DEPS="dask-core=1.* gdal=2.* toolz pandas=0.25 geopandas=0.6 pyproj>=2"
+  - os: linux
     env: PYTHON_VERSION="3.7" MAINDEPS="numpy scipy pytz" DEPS="dask-core gdal=2.* toolz geopandas pyproj>=2"
 
 install:

CHANGES.rst

@@ -4,7 +4,7 @@ Changelog of dask-geomodeling
 2.2.10 (unreleased)
 -------------------
 
-- Nothing changed yet.
+- Fix point requests in raster.Smooth.
 
 
 2.2.9 (2020-06-23)

dask_geomodeling/geometry/sources.py

@@ -3,6 +3,7 @@
 """
 import fiona
 import geopandas as gpd
+import warnings
 
 from dask import config
 from dask_geomodeling import utils
@@ -133,7 +134,9 @@ def process(url, request):
 
         # only return geometries that truly intersect the requested geometry
         if request["mode"] == "centroid":
-            f = f[f["geometry"].centroid.within(filt_geom.iloc[0])]
+            with warnings.catch_warnings():    # geopandas warns if in WGS84
+                warnings.simplefilter("ignore")
+                f = f[f["geometry"].centroid.within(filt_geom.iloc[0])]
         else:
             f = f[f["geometry"].intersects(filt_geom.iloc[0])]
 
@@ -237,7 +240,10 @@ def process(data, request):
             return {"features": f, "projection": request['projection']}
 
         elif mode == 'centroid':
-            if not geometry.centroid.intersects(request["geometry"]):
+            with warnings.catch_warnings():    # geopandas warns if in WGS84
+                warnings.simplefilter("ignore")
+                centroid = geometry.centroid
+            if not centroid.intersects(request["geometry"]):
                 return {
                     "projection": request["projection"],
                     "features": gpd.GeoDataFrame([]),

dask_geomodeling/raster/spatial.py

@@ -77,17 +77,20 @@ def expand_request_meters(request, radius_m=1):
     x1, y1, x2, y2 = bbox
     shape_m = y2 - y1, x2 - x1
 
-    # omit the +1 in zoom for efficiency: zoom=0.2 is a zoom factor of 1.2
-    zoom = [2 * radius_m / s for s in shape_m]
-
-    # compute the size in pixels, and compute the margins (rounded size)
-    shape_px = request["height"], request["width"]
-    radius_px = [z * s / 2 for (z, s) in zip(zoom, shape_px)]
-    margins_px = [int(round(sz)) for sz in radius_px]
-
-    # use these (int-valued) margins to compute the actual zoom and margins
-    zoom = [2 * m / s for (s, m) in zip(shape_px, margins_px)]
-    margins_m = [z * s / 2 for (z, s) in zip(zoom, shape_m)]
+    if shape_m[0] > 0 and shape_m[1] > 0:
+        # Resolution in pixels per meter:
+        resolution = request["height"] / shape_m[0], request["width"] / shape_m[1]
+        # How many pixels to add:
+        radius_px = [radius_m * res for res in resolution]
+        # How many pixels to add, rounded to integers:
+        margins_px = [int(round(r)) for r in radius_px]
+        # How many meters to add (based on rounded pixels):
+        margins_m = [m / res for m, res in zip(margins_px, resolution)]
+    else:
+        # There is no resolution. Add MARGIN_THRESHOLD pixels to the request.
+        radius_px = margins_px = [Smooth.MARGIN_THRESHOLD] * 2
+        # Expand the request with radius_m exactly.
+        margins_m = [radius_m] * 2
 
     # assemble the request
     new_request = request.copy()

dask_geomodeling/tests/test_raster.py

@@ -1888,6 +1888,13 @@ def test_smooth(self):
         self.assertEqual(view.get_data(**self.meta_request)["meta"], self.expected_meta)
         self.assertEqual(view.get_data(**self.time_request)["time"], self.expected_time)
 
+        # point request
+        request["bbox"] = (50, 50, 50, 50)
+        request["height"] = 1
+        request["width"] = 1
+        data = view.get_data(**request)
+        data["values"][0, 0, 0] == peak
+
     def test_hill_shade(self):
         view = raster.HillShade(store=self.raster)
         self.assertEqual(view.dtype, "u1")

dask_geomodeling/tests/test_raster_spatial.py

@@ -81,7 +81,7 @@ def test_place_attrs(source, center):
 
 def test_place_invalid_statistic(source, center):
     with pytest.raises(ValueError):
-        place = raster.Place(
+        raster.Place(
             source, "EPSG:28992", center, [(50, 50)], statistic="nonexisting"
         )
 

dask_geomodeling/utils.py

@@ -475,7 +475,9 @@ def transform_min_size(min_size, geometry, src_srs, dst_srs):
     Note that in order to guarantee the minimum size, the reverse operation may
     result in a larger value than the original minimum size.
     """
-    source = geometry.centroid.buffer(min_size / 2)
+    with warnings.catch_warnings():  # geopandas warns if in WGS84
+        warnings.simplefilter("ignore")
+        source = geometry.centroid.buffer(min_size / 2)
     target = shapely_transform(source, src_srs=src_srs, dst_srs=dst_srs)
     x1, y1, x2, y2 = target.bounds
     return max(x2 - x1, y2 - y1)