Add decimation base as input to cog_translate

docs/docs/Advanced.md

@@ -41,6 +41,20 @@ print(overviews)
 [2, 4, 8]
 ```
 
+### Decimation Base
+
+As described above, a decimation base of 2 is used by default. However you can provide a custom base, N > 1, with *--decimation-base N*. Optimal overviews are computed assuming a base 2 is used, so using *--decimation-base* also requires that *--overview-level* is provided. Similar to the default example, here are the overviews for base 3:
+
+```python
+overview_level = 3
+decimation_base = 3
+overviews = [decimation_base ** j for j in range(1, overview_level + 1)]
+print(overviews)
+[3, 9, 27]
+```
+
+This is primarily useful when working with [custom TileMatrixSets](https://developmentseed.org/morecantile/usage/#define-custom-grid) that also use a non-default decimation base.
+
 ## Band metadata
 By default rio cogeo DO NOT forward **band** metadata (e.g statistics) to the output dataset.
 

rio_cogeo/cogeo.py

@@ -98,6 +98,7 @@ def cog_translate(  # noqa: C901
     colormap: Optional[Dict] = None,
     additional_cog_metadata: Optional[Dict] = None,
     use_cog_driver: bool = False,
+    decimation_base: int = 2,
 ):
     """
     Create Cloud Optimized Geotiff.
@@ -171,6 +172,10 @@ def cog_translate(  # noqa: C901
         Additional dataset metadata to add to the COG.
     use_cog_driver: bool, optional (default: False)
         Use GDAL COG driver if set to True. COG driver is available starting with GDAL 3.1.
+    decimation_base: int, default: 2
+        How overviews are divided at each zoom level (default is 2). Must be greater than 1.
+        Also requires that `overview_level` is provided for `decimation_base` values greater
+        than 2.
 
     .. deprecated:: 5.1.2
         `web_optimized` is deprecated in favor of `tms`.
@@ -187,6 +192,15 @@ def cog_translate(  # noqa: C901
         )
         tms = tms or morecantile.tms.get("WebMercatorQuad")
 
+    if decimation_base <= 1:
+        raise ValueError(
+            "Decimation base must be greater than 1 for building overviews."
+        )
+    elif decimation_base > 2 and overview_level is None:
+        raise ValueError(
+            "Decimation base values greater than 2 require that overview_level is defined."
+        )
+
     dst_kwargs = dst_kwargs.copy()
 
     if isinstance(indexes, int):
@@ -343,7 +357,7 @@ def cog_translate(  # noqa: C901
                 if not quiet and overview_level:
                     click.echo("Adding overviews...", err=True)
 
-                overviews = [2**j for j in range(1, overview_level + 1)]
+                overviews = [decimation_base**j for j in range(1, overview_level + 1)]
                 tmp_dst.build_overviews(overviews, ResamplingEnums[overview_resampling])
 
                 if not quiet:

rio_cogeo/scripts/cli.py

@@ -122,6 +122,14 @@ def cogeo():
     "selected until the smallest overview is smaller than the value of the "
     "internal blocksize)",
 )
+@click.option(
+    "--decimation-base",
+    type=int,
+    help="Decimation base, how to sub-divide a raster for each overview level. "
+    "Also requires that --overview-level is provided.",
+    default=2,
+    show_default=True,
+)
 @click.option(
     "--overview-resampling",
     help="Overview creation resampling algorithm.",
@@ -230,6 +238,7 @@ def create(
     add_mask,
     blocksize,
     overview_level,
+    decimation_base,
     overview_resampling,
     overview_blocksize,
     web_optimized,
@@ -318,6 +327,7 @@ def create(
         tms=tilematrixset,
         use_cog_driver=use_cog_driver,
         quiet=quiet,
+        decimation_base=decimation_base,
     )
 
 

tests/test_cli.py

@@ -210,6 +210,28 @@ def test_cogeo_validOvrOption(runner):
             assert src.overviews(1) == [2, 4]
 
 
+def test_cogeo_decimation_base_option(runner):
+    """Should work as expected."""
+    with runner.isolated_filesystem():
+        result = runner.invoke(
+            cogeo,
+            [
+                "create",
+                raster_path_rgb,
+                "output.tif",
+                "--decimation-base",
+                3,
+                "--overview-level",
+                6,
+            ],
+        )
+        assert not result.exception
+        assert result.exit_code == 0
+        with rasterio.open("output.tif") as src:
+            assert len(src.overviews(1)) == 6
+            assert src.overviews(1)[0] == 3
+
+
 def test_cogeo_overviewTilesize(monkeypatch, runner):
     """Should work as expected."""
     with runner.isolated_filesystem():

tests/test_cogeo.py

@@ -724,3 +724,23 @@ def test_cog_translate_forward_ns_metadata(runner):
             with rasterio.open("cogeo.tif") as src:
                 assert src.tags(ns="IMD")
                 assert src.tags(ns="RPC")
+
+
+def test_cog_translate_decimation_base(runner):
+    """Should create proper overviews when using custom decimation base."""
+    with runner.isolated_filesystem():
+
+        base_level_pairs = [(3, 6), (4, 5), (5, 4)]
+
+        for decimation_base, overview_level in base_level_pairs:
+            cog_translate(
+                raster_path_rgb,
+                "cogeo.tif",
+                cog_profiles.get("deflate"),
+                decimation_base=decimation_base,
+                overview_level=overview_level,
+                quiet=True,
+            )
+
+            with rasterio.open("cogeo.tif") as src:
+                assert src.overviews(1)[0] == decimation_base