Custom TMS with decimation

morecantile/models.py

@@ -660,6 +660,7 @@ def custom(
         ordered_axes: Optional[List[str]] = None,
         geographic_crs: pyproj.CRS = WGS84_CRS,
         screen_pixel_size: float = 0.28e-3,
+        decimation_base: int = 2,
         **kwargs: Any,
     ):
         """
@@ -694,6 +695,8 @@ def custom(
             Geographic (lat,lon) coordinate reference system (default is EPSG:4326)
         screen_pixel_size: float, optional
             Rendering pixel size. 0.28 mm was the actual pixel size of a common display from 2005 and considered as standard by OGC.
+        decimation_base: int, optional
+            How tiles are divided at each zoom level (default is 2). Must be greater than 1.
         kwargs: Any
             Attributes to forward to the TileMatrixSet
 
@@ -713,6 +716,11 @@ def custom(
             transform = pyproj.Transformer.from_crs(extent_crs, crs, always_xy=True)
             extent = transform.transform_bounds(*extent, densify_pts=21)
 
+        if decimation_base <= 1:
+            raise ValueError(
+                "Custom TileMatrixSet requires a decimation base that is greater than 1."
+            )
+
         bbox = BoundingBox(*extent)
         x_origin = bbox.left if not is_inverted else bbox.top
         y_origin = bbox.top if not is_inverted else bbox.left
@@ -723,8 +731,10 @@ def custom(
         tile_matrices: List[TileMatrix] = []
         for zoom in range(minzoom, maxzoom + 1):
             res = max(
-                width / (tile_width * matrix_scale[0]) / 2.0**zoom,
-                height / (tile_height * matrix_scale[1]) / 2.0**zoom,
+                width / (tile_width * matrix_scale[0]) / float(decimation_base) ** zoom,
+                height
+                / (tile_height * matrix_scale[1])
+                / float(decimation_base) ** zoom,
             )
             tile_matrices.append(
                 TileMatrix(
@@ -735,8 +745,8 @@ def custom(
                         "pointOfOrigin": [x_origin, y_origin],
                         "tileWidth": tile_width,
                         "tileHeight": tile_height,
-                        "matrixWidth": matrix_scale[0] * 2**zoom,
-                        "matrixHeight": matrix_scale[1] * 2**zoom,
+                        "matrixWidth": matrix_scale[0] * decimation_base**zoom,
+                        "matrixHeight": matrix_scale[1] * decimation_base**zoom,
                     }
                 )
             )

tests/test_models.py

@@ -185,6 +185,31 @@ def test_custom_tms_bounds_user_crs():
     assert custom_tms.bounds(0, 0, 0) == (-120, 30, -110, 40)
 
 
+def test_custom_tms_decimation():
+    """Check bounds with epsg6342 and custom decimation base."""
+    extent = (238170, 4334121, 377264, 4473215)
+    left, bottom, right, top = extent
+    for decimation_base in [2, 3, 4, 5]:
+        custom_tms = TileMatrixSet.custom(
+            extent,
+            pyproj.CRS.from_epsg(6342),
+            decimation_base=decimation_base,
+        )
+
+        if decimation_base == 2:
+            assert custom_tms.is_quadtree
+        else:
+            assert not custom_tms.is_quadtree
+
+        for zoom in [0, 1, 2, 3]:
+            tile_width = (right - left) / decimation_base**zoom
+            tile_height = (top - bottom) / decimation_base**zoom
+            expected = (left, top - tile_height, left + tile_width, top)
+            tile_bounds = custom_tms.xy_bounds(0, 0, zoom)
+            for a, b in zip(expected, tile_bounds):
+                assert round(a - b, 4) == 0
+
+
 def test_nztm_quad_is_quad():
     """Test NZTM2000Quad."""
     tms = morecantile.tms.get("NZTM2000Quad")