Merge pull request #62 from cogeotiff/webOpt

Web-Optimized COG: COGEO optimized for dynamic tiling
cogeotiff · Apr 16, 2019 · c5d7344 · c5d7344
2 parents 48b015b + cf8bdcb
commit c5d7344
Show file tree

Hide file tree

Showing 12 changed files with 595 additions and 84 deletions.
diff --git a/CHANGES.txt b/CHANGES.txt
@@ -1,3 +1,11 @@
+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)
+- add `--resampling` option to select the resampling algorithm when using `--web-optimized`
+
 
 1.0b3 (2019-03-30)
 ------------------
@@ -41,6 +49,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)

diff --git a/README.rst b/README.rst
@@ -50,51 +50,56 @@ 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] Overview creation 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.
+        -r, --resampling [nearest|bilinear|cubic|cubic_spline|lanczos|average|mode|gauss] Resampling algorithm.
+        --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 +172,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 +218,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 +317,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.
diff --git a/rio_cogeo/cogeo.py b/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.enums import Resampling
+from rasterio.warp import transform_bounds
+from rasterio.enums import Resampling as ResamplingEnums
 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,9 @@ def cog_translate(
     add_mask=None,
     overview_level=None,
     overview_resampling="nearest",
+    web_optimized=False,
+    latitude_adjustment=True,
+    resampling="nearest",
     config=None,
     quiet=False,
 ):
@@ -40,7 +54,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 +65,12 @@ 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.
+    resampling : str, optional (default: "nearest")
+        Resampling algorithm.
     config : dict
         Rasterio Env options.
     quiet: bool, optional (default: False)
@@ -77,12 +97,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 +107,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=ResamplingEnums[resampling],
+                    )
+                )
+
             with WarpedVRT(src_dst, **vrt_params) as vrt_dst:
                 meta = vrt_dst.meta
                 meta["count"] = len(indexes)
@@ -125,8 +174,20 @@ 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])
+                        mem.build_overviews(
+                            overviews, ResamplingEnums[overview_resampling]
+                        )
 
                         if not quiet:
                             click.echo("Updating dataset tags...", err=True)
@@ -137,7 +198,7 @@ def cog_translate(
                         tags = src_dst.tags()
                         tags.update(
                             dict(
-                                OVR_RESAMPLING_ALG=Resampling[
+                                OVR_RESAMPLING_ALG=ResamplingEnums[
                                     overview_resampling
                                 ].name.upper()
                             )