Web-Optimized COG: COGEO optimized for dynamic tiling

CHANGES.txt

@@ -1,3 +1,10 @@
+Next (TBD)
+----------
+
+- add `--web-optimized` option to create a web optimized COG (#10)
+- add `--latitude-adjustment/--global-maxzoom` option to adjust MAX_ZOOM for global datasets
+- Web-optimized tests needs python3.6 (cogdumper)
+
 
 1.0b3 (2019-03-30)
 ------------------
@@ -41,6 +48,7 @@ Breacking Changes:
 - internal nodata or alpha channel can be forwarded to the output dataset.
 - removed default overview blocksize to be equal to the raw data blocksize (#60)
 
+
 1.0dev10 (2019-02-12)
 ---------------------
 - allow non integer nodata value (#51)

README.rst

@@ -50,51 +50,55 @@ Usage
 .. code-block:: console
 
   $ rio cogeo --help
-  Usage: rio cogeo [OPTIONS] COMMAND [ARGS]...
+    Usage: rio cogeo [OPTIONS] COMMAND [ARGS]...
 
     Rasterio cogeo subcommands.
 
-  Options:
-    --help  Show this message and exit.
+    Options:
+      --help  Show this message and exit.
 
-  Commands:
-    create    Create COGEO
-    validate  Validate COGEO
+    Commands:
+      create    Create COGEO
+      validate  Validate COGEO
 
 - Create a Cloud Optimized Geotiff.
 
 .. code-block:: console
 
-  $ rio cogeo --help
-  Usage: rio cogeo [OPTIONS] INPUT OUTPUT
-
-    Create Cloud Optimized Geotiff.
-
-  Options:
-    -b, --bidx BIDX                 Band indexes to copy.
-    -p, --cog-profile [jpeg|webp|zstd|lzw|deflate|packbits|raw] CloudOptimized GeoTIFF profile (default: deflate).
-    --nodata NUMBER|nan             Set nodata masking values for input dataset.
-    --add-mask                      Force output dataset creation with an internal mask (convert alpha band or nodata to mask).
-    --overview-level INTEGER        Overview level (if not provided, appropriate overview level will be selected until the
-                                    smallest overview is smaller than the internal block size).
-    --overview-resampling [nearest|bilinear|cubic|cubic_spline|lanczos|average|mode|gauss] Resampling algorithm.
-    --overview-blocksize TEXT       Overview's internal tile size (default defined by GDAL_TIFF_OVR_BLOCKSIZE env or 128)
-    --threads INTEGER
-    --co, --profile NAME=VALUE      Driver specific creation options.See the documentation for the selected output driver for more information.
-    -q, --quiet                     Suppress progress bar and other non-error output.
-    --help                          Show this message and exit.
+    $ rio cogeo create --help
+      Usage: rio cogeo create [OPTIONS] INPUT OUTPUT
+
+      Create Cloud Optimized Geotiff.
+
+      Options:
+        -b, --bidx BIDX                 Band indexes to copy.
+        -p, --cog-profile [jpeg|webp|zstd|lzw|deflate|packbits|raw] CloudOptimized GeoTIFF profile (default: deflate).
+        --nodata NUMBER|nan             Set nodata masking values for input dataset.
+        --add-mask                      Force output dataset creation with an internal mask (convert alpha band or nodata to mask).
+        --overview-level INTEGER        Overview level (if not provided, appropriate overview level will be selected
+                                        until the smallest overview is smaller than the value of the internal blocksize)
+        --overview-resampling [nearest|bilinear|cubic|cubic_spline|lanczos|average|mode|gauss] Resampling algorithm.
+        --overview-blocksize TEXT       Overview's internal tile size (default defined by GDAL_TIFF_OVR_BLOCKSIZE env or 128)
+        -w, --web-optimized             Create COGEO optimized for Web.
+        --latitude-adjustment / --global-maxzoom
+                                        Use dataset native mercator resolution for MAX_ZOOM calculation (linked to dataset center latitude, default)
+                                        or ensure MAX_ZOOM equality for multiple dataset accross latitudes.
+        --threads INTEGER
+        --co, --profile NAME=VALUE      Driver specific creation options.See the documentation for the selected output driver for more information.
+        -q, --quiet                     Remove progressbar and other non-error output.
+        --help                          Show this message and exit.
 
 - Check if a Cloud Optimized Geotiff is valid.
 
 .. code-block:: console
 
   $ rio cogeo validate --help
-  Usage: rio cogeo validate [OPTIONS] INPUT
+    Usage: rio cogeo validate [OPTIONS] INPUT
 
     Validate Cloud Optimized Geotiff.
 
-  Options:
-    --help  Show this message and exit.
+    Options:
+      --help  Show this message and exit.
 
 
 Examples
@@ -167,20 +171,25 @@ Profiles can be extended by providing '--co' option in command line
     $ rio cogeo create mydataset.tif mydataset_raw.tif --co BLOCKXSIZE=1024 --co BLOCKYSIZE=1024 --cog-profile raw --overview-blocksize 256
 
 
-Overview levels
-===============
+Web-Optimized COG
+=================
 
-By default rio cogeo will calculate the optimal overview level based on dataset
-size and internal tile size (overview should not be smaller than internal tile
-size (e.g 512px). Overview level will be translated to decimation level of
-power of two:
+rio-cogeo provide a *--web-optimized* option which aims to create a web-tiling friendly COG.
 
-.. code-block:: python
+Output dataset features:
+
+- bounds and internal tiles aligned with web-mercator grid.
+- raw data and overviews resolution match mercator zoom level resolution.
+
+**Important**
+
+Because the mercator project does not respect the distance, when working with
+multiple images covering different latitudes, you may want to use the *--global-maxzoom* option
+to create output dataset having the same MAX_ZOOM (raw data resolution).
+
+Because it will certainly create a larger file, a nodata value or alpha band should
+be present in the input dataset. If not the original data will be surrounded by black (0) data.
 
-  overview_level = 3
-  overviews = [2 ** j for j in range(1, overview_level + 1)]
-  print(overviews)
-  [2, 4, 8]
 
 Internal tile size
 ==================
@@ -208,13 +217,32 @@ GDAL configuration to merge consecutive range requests
     GDAL_HTTP_VERSION=2
 
 
+Overview levels
+===============
+
+By default rio cogeo will calculate the optimal overview level based on dataset
+size and internal tile size (overview should not be smaller than internal tile
+size (e.g 512px). Overview level will be translated to decimation level of
+power of two:
+
+.. code-block:: python
+
+  overview_level = 3
+  overviews = [2 ** j for j in range(1, overview_level + 1)]
+  print(overviews)
+  [2, 4, 8]
+
+
 GDAL Version
 ============
 
 It is recommanded to use GDAL > 2.3.2. Previous version might not be able to
 create proper COGs (ref: https://github.com/OSGeo/gdal/issues/754).
 
 
+More info in https://github.com/cogeotiff/rio-cogeo/issues/60
+
+
 Nodata, Alpha and Mask
 ======================
 
@@ -288,9 +316,10 @@ This repo is set to use `pre-commit` to run *flake8*, *pydocstring* and *black*
 
   $ pre-commit install
 
+
 Extras
 ======
 
 Blog post on good and bad COG formats: https://medium.com/@_VincentS_/do-you-really-want-people-using-your-data-ec94cd94dc3f
 
-Checkout `**rio-glui** <https://github.com/mapbox/rio-glui/>__` rasterio plugin to explore COG locally in your web browser.
+Checkout `rio-glui <https://github.com/mapbox/rio-glui/>`__ rasterio plugin to explore COG locally in your web browser.

rio_cogeo/cogeo.py

@@ -10,11 +10,22 @@
 from rasterio.io import MemoryFile
 from rasterio.env import GDALVersion
 from rasterio.vrt import WarpedVRT
+from rasterio.warp import transform_bounds
 from rasterio.enums import Resampling
 from rasterio.shutil import copy
+from rasterio.transform import Affine
+
+import mercantile
+from supermercado.burntiles import tile_extrema
 
 from rio_cogeo.errors import LossyCompression
-from rio_cogeo.utils import get_maximum_overview_level, has_alpha_band, has_mask_band
+from rio_cogeo.utils import (
+    get_maximum_overview_level,
+    has_alpha_band,
+    has_mask_band,
+    get_max_zoom,
+    _meters_per_pixel,
+)
 
 
 def cog_translate(
@@ -26,6 +37,8 @@ def cog_translate(
     add_mask=None,
     overview_level=None,
     overview_resampling="nearest",
+    web_optimized=False,
+    latitude_adjustment=True,
     config=None,
     quiet=False,
 ):
@@ -40,7 +53,7 @@ def cog_translate(
         An output dataset path or or PathLike object.
         Will be opened in "w" mode.
     dst_kwargs: dict
-        output dataset creation options.
+        Output dataset creation options.
     indexes : tuple, int, optional
         Raster band indexes to copy.
     nodata, int, optional
@@ -51,6 +64,10 @@ def cog_translate(
         COGEO overview (decimation) level
     overview_resampling : str, optional (default: "nearest")
         Resampling algorithm for overviews
+    web_optimized: bool, option (default: False)
+        Create web-optimized cogeo.
+    latitude_adjustment: bool, option (default: True)
+        Use mercator meters for zoom calculation or ensure max zoom equality.
     config : dict
         Rasterio Env options.
     quiet: bool, optional (default: False)
@@ -77,12 +94,6 @@ def cog_translate(
                     LossyCompression,
                 )
 
-            if overview_level is None:
-                overview_level = get_maximum_overview_level(
-                    src_dst,
-                    min(int(dst_kwargs["blockxsize"]), int(dst_kwargs["blockysize"])),
-                )
-
             vrt_params = dict(add_alpha=True)
 
             if nodata is not None:
@@ -93,6 +104,41 @@ def cog_translate(
             if alpha:
                 vrt_params.update(dict(add_alpha=False))
 
+            if web_optimized:
+                bounds = list(
+                    transform_bounds(
+                        *[src_dst.crs, "epsg:4326"] + list(src_dst.bounds),
+                        densify_pts=21
+                    )
+                )
+                center = [(bounds[0] + bounds[2]) / 2, (bounds[1] + bounds[3]) / 2]
+
+                tilesize = min(
+                    int(dst_kwargs["blockxsize"]), int(dst_kwargs["blockysize"])
+                )
+                lat = 0 if latitude_adjustment else center[1]
+                max_zoom = get_max_zoom(src_dst, lat=lat, tilesize=tilesize)
+
+                extrema = tile_extrema(bounds, max_zoom)
+                w, n = mercantile.xy(
+                    *mercantile.ul(extrema["x"]["min"], extrema["y"]["min"], max_zoom)
+                )
+                vrt_res = _meters_per_pixel(max_zoom, 0, tilesize=tilesize)
+                vrt_transform = Affine(vrt_res, 0, w, 0, -vrt_res, n)
+
+                vrt_width = (extrema["x"]["max"] - extrema["x"]["min"]) * tilesize
+                vrt_height = (extrema["y"]["max"] - extrema["y"]["min"]) * tilesize
+
+                vrt_params.update(
+                    dict(
+                        crs="epsg:3857",
+                        transform=vrt_transform,
+                        width=vrt_width,
+                        height=vrt_height,
+                        resampling=Resampling[overview_resampling],
+                    )
+                )
+
             with WarpedVRT(src_dst, **vrt_params) as vrt_dst:
                 meta = vrt_dst.meta
                 meta["count"] = len(indexes)
@@ -125,6 +171,16 @@ def cog_translate(
 
                         if not quiet:
                             click.echo("Adding overviews...", err=True)
+
+                        if overview_level is None:
+                            overview_level = get_maximum_overview_level(
+                                vrt_dst,
+                                min(
+                                    int(dst_kwargs["blockxsize"]),
+                                    int(dst_kwargs["blockysize"]),
+                                ),
+                            )
+
                         overviews = [2 ** j for j in range(1, overview_level + 1)]
                         mem.build_overviews(overviews, Resampling[overview_resampling])
 

rio_cogeo/scripts/cli.py

@@ -1,6 +1,7 @@
 """Rio_cogeo.scripts.cli."""
 
 import os
+import warnings
 
 import click
 import numpy
@@ -76,12 +77,15 @@ def cogeo():
 @click.option(
     "--add-mask",
     is_flag=True,
-    help="Force output dataset creation with an internal mask (convert alpha band or nodata to mask).",
+    help="Force output dataset creation with an internal mask (convert alpha "
+    "band or nodata to mask).",
 )
 @click.option(
     "--overview-level",
     type=int,
-    help="Overview level (if not provided, appropriate overview level will be selected until the smallest overview is smaller than the internal block size).",
+    help="Overview level (if not provided, appropriate overview level will be "
+    "selected until the smallest overview is smaller than the value of the "
+    "internal blocksize)",
 )
 @click.option(
     "--overview-resampling",
@@ -94,15 +98,23 @@ def cogeo():
 @click.option(
     "--overview-blocksize",
     default=lambda: os.environ.get("GDAL_TIFF_OVR_BLOCKSIZE", 128),
-    help="Overview's internal tile size (default defined by GDAL_TIFF_OVR_BLOCKSIZE env or 128)",
+    help="Overview's internal tile size (default defined by "
+    "GDAL_TIFF_OVR_BLOCKSIZE env or 128)",
+)
+@click.option(
+    "--web-optimized", "-w", is_flag=True, help="Create COGEO optimized for Web."
+)
+@click.option(
+    "--latitude-adjustment/--global-maxzoom",
+    default=None,
+    help="Use dataset native mercator resolution for MAX_ZOOM calculation "
+    "(linked to dataset center latitude, default) or ensure MAX_ZOOM equality for multiple "
+    "dataset accross latitudes.",
 )
 @click.option("--threads", type=int, default=8)
 @options.creation_options
 @click.option(
-    "--quiet",
-    "-q",
-    help="Suppress progress bar and other non-error output.",
-    is_flag=True,
+    "--quiet", "-q", help="Remove progressbar and other non-error output.", is_flag=True
 )
 def create(
     input,
@@ -114,11 +126,19 @@ def create(
     overview_level,
     overview_resampling,
     overview_blocksize,
+    web_optimized,
+    latitude_adjustment,
     threads,
     creation_options,
     quiet,
 ):
     """Create Cloud Optimized Geotiff."""
+    if latitude_adjustment is not None and not web_optimized:
+        warnings.warn(
+            "'latitude_adjustment' option has to be used with --web-optimized options. "
+            "Will be ignored."
+        )
+
     output_profile = cog_profiles.get(cogeo_profile)
     output_profile.update(dict(BIGTIFF=os.environ.get("BIGTIFF", "IF_SAFER")))
     if creation_options:
@@ -139,6 +159,8 @@ def create(
         add_mask,
         overview_level,
         overview_resampling,
+        web_optimized,
+        latitude_adjustment,
         config,
         quiet,
     )