Merge pull request #485 from ungarj/shapely_1.8.4_fix

adapt code to shapely 1.8.4

mapchete/tile.py

@@ -397,12 +397,8 @@ def count_tiles(
     unbuffered_pyramid = BufferedTilePyramid(
         pyramid.grid, tile_size=pyramid.tile_size, metatiling=pyramid.metatiling
     )
-    # make sure no rounding errors occur
-    geometry = geometry.buffer(-0.000000001)
-
     height = pyramid.matrix_height(init_zoom)
     width = pyramid.matrix_width(init_zoom)
-
     # rasterize to array and count cells if too many tiles are expected
     if width > rasterize_threshold or height > rasterize_threshold:
         logger.debug("rasterize tiles to count geometry overlap")
@@ -428,18 +424,18 @@ def _count_tiles(tiles, geometry, minzoom, maxzoom):
     count = 0
     for tile in tiles:
         # determine data covered by tile
-        tile_intersection = tile.bbox.intersection(geometry)
+        tile_intersection = geometry.intersection(tile.bbox)
 
-        # skip if there is no data
-        if tile_intersection.is_empty:
+        # skip if there is no intersection
+        if not tile_intersection.area:
             continue
+
         # increase counter as tile contains data
-        elif tile.zoom >= minzoom:
+        if tile.zoom >= minzoom:
             count += 1
 
         # if there are further zoom levels, analyze descendants
         if tile.zoom < maxzoom:
-
             # if tile is half full, analyze each descendant
             # also do this if the tile children are not four in which case we cannot use
             # the count formula below
@@ -467,18 +463,27 @@ def _count_tiles(tiles, geometry, minzoom, maxzoom):
 
 
 def _count_cells(pyramid, geometry, minzoom, maxzoom):
-    # rasterize geometry on maxzoom
     if geometry.is_empty:
         return 0
 
+    # for the rasterization algorithm we need to keep the all_touched flag True
+    # but slightly reduce the geometry area in order to get the same results as
+    # with the tiles/vector algorithm.
+    left, bottom, right, top = geometry.bounds
+    width = right - left
+    height = top - bottom
+    buffer_distance = ((width + height) / 2) * -0.0000001
+    # geometry will be reduced by a tiny fraction of the average from bounds width & height
+    geometry_reduced = geometry.buffer(buffer_distance)
+
     logger.debug(
         "rasterize polygon on %s x %s cells",
         pyramid.matrix_height(maxzoom),
         pyramid.matrix_width(maxzoom),
     )
     transform = pyramid.matrix_affine(maxzoom)
     raster = rasterize(
-        [(geometry, 1)],
+        [(geometry_reduced, 1)],
         out_shape=(pyramid.matrix_height(maxzoom), pyramid.matrix_width(maxzoom)),
         fill=0,
         transform=transform,

requirements.txt

@@ -21,7 +21,7 @@ pyproj
 retry>=0.9.2
 rasterio>=1.0.28
 s3fs>=0.5.2
-Shapely>=1.7.1,<1.8.1
+Shapely>=1.7.1
 tilematrix>=2022.3.0
 tqdm>=4.46.0
 werkzeug>=0.15

test/test_mapchete.py

@@ -514,7 +514,8 @@ def test_process_template(dummy1_tif, mp_tmpdir):
         mp.execute(process_tile)
 
 
-def test_count_tiles(zoom_mapchete):
+@pytest.mark.parametrize("minzoom", range(0, 14))
+def test_count_tiles(zoom_mapchete, minzoom):
     """Count tiles function."""
     maxzoom = 13
     conf = zoom_mapchete.dict
@@ -524,37 +525,39 @@ def test_count_tiles(zoom_mapchete):
         input=None,
     )
     conf["pyramid"].update(metatiling=8, pixelbuffer=5)
-    for minzoom in range(0, 14):
-        conf["zoom_levels"].update(min=minzoom)
-        with mapchete.open(conf) as mp:
-            assert len(list(mp.get_process_tiles())) == mapchete.count_tiles(
-                mp.config.area_at_zoom(), mp.config.process_pyramid, minzoom, maxzoom
-            )
-
-
-def test_count_tiles_mercator():
-    for metatiling in [1, 2, 4, 8, 16]:
-        tp = BufferedTilePyramid("mercator", metatiling=metatiling)
-        for zoom in range(13):
-            count_by_geom = count_tiles(box(*tp.bounds), tp, zoom, zoom)
-            count_by_tp = tp.matrix_width(zoom) * tp.matrix_height(zoom)
-            assert count_by_geom == count_by_tp
+    conf["zoom_levels"].update(min=minzoom)
+    with mapchete.open(conf) as mp:
+        assert len(list(mp.get_process_tiles())) == mapchete.count_tiles(
+            mp.config.area_at_zoom(), mp.config.process_pyramid, minzoom, maxzoom
+        )
 
 
-def test_count_tiles_large_init_zoom(geometrycollection):
-    tp = BufferedTilePyramid(
-        grid=dict(
-            shape=[100, 200],
-            bounds=[-180, -90, 180, 90],
-            srs=dict(epsg=4326),
-            is_global=True,
-        )
+@pytest.mark.parametrize("metatiling", [1, 2, 4, 8, 16])
+@pytest.mark.parametrize("zoom", range(15))
+def test_count_tiles_mercator(metatiling, zoom):
+    tp = BufferedTilePyramid("mercator", metatiling=metatiling)
+    count_by_geom = count_tiles(box(*tp.bounds), tp, zoom, zoom)
+    count_by_tp = tp.matrix_width(zoom) * tp.matrix_height(zoom)
+    assert count_by_geom == count_by_tp
+
+
+# This test only works until zoom 14. After this the results between the count_tiles()
+# algorithms (rasterized & tile-based) start to differ from the actual TilePyramid.tiles_from_geom()
+# implementation. Please also note that TilePyramid.tiles_from_geom(exact=True) ast to be activated
+# in order to pass
+@pytest.mark.parametrize("zoom", range(14))
+def test_count_tiles_large_init_zoom(geometrycollection, zoom):
+    tp = BufferedTilePyramid(grid="geodetic")
+    raster_count = count_tiles(
+        geometrycollection, tp, zoom, zoom, rasterize_threshold=0
     )
-    minzoom = 0
-    maxzoom = 7
-    assert count_tiles(
-        geometrycollection, tp, minzoom, maxzoom, rasterize_threshold=0
-    ) == count_tiles(geometrycollection, tp, minzoom, maxzoom, rasterize_threshold=1000)
+    tile_count = count_tiles(
+        geometrycollection, tp, zoom, zoom, rasterize_threshold=100_000_000
+    )
+    tiles = len(
+        list(tp.tiles_from_geom(geometry=geometrycollection, zoom=zoom, exact=True))
+    )
+    assert raster_count == tile_count == tiles
 
 
 def test_batch_process(mp_tmpdir, cleantopo_tl):
@@ -605,21 +608,21 @@ def test_execute_kwargs(example_mapchete, execute_kwargs_py):
         mp.execute((7, 61, 129))
 
 
-def test_snap_bounds_to_zoom():
+@pytest.mark.parametrize("pixelbuffer", [0, 5, 10])
+@pytest.mark.parametrize("metatiling", [1, 2, 4])
+@pytest.mark.parametrize("zoom", range(3, 5))
+def test_snap_bounds_to_zoom(pixelbuffer, metatiling, zoom):
     bounds = (-180, -90, -60, -30)
-    for pixelbuffer in [0, 5, 10]:
-        for metatiling in [1, 2, 4]:
-            pyramid = BufferedTilePyramid(
-                "geodetic", pixelbuffer=pixelbuffer, metatiling=metatiling
-            )
-            for zoom in range(3, 5):
-                snapped_bounds = mapchete.config.snap_bounds(
-                    bounds=bounds, pyramid=pyramid, zoom=zoom
-                )
-                control_bounds = unary_union(
-                    [t.bbox for t in pyramid.tiles_from_bounds(bounds, zoom)]
-                ).bounds
-                assert snapped_bounds == control_bounds
+    pyramid = BufferedTilePyramid(
+        "geodetic", pixelbuffer=pixelbuffer, metatiling=metatiling
+    )
+    snapped_bounds = mapchete.config.snap_bounds(
+        bounds=bounds, pyramid=pyramid, zoom=zoom
+    )
+    control_bounds = unary_union(
+        [t.bbox for t in pyramid.tiles_from_bounds(bounds, zoom)]
+    ).bounds
+    assert snapped_bounds == control_bounds
 
 
 def test_snap_bounds_errors():