diff --git a/CHANGES.txt b/CHANGES.txt
index 22b6258..540010b 100644
--- 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
index d523525..c6be445 100644
--- 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,6 +218,22 @@ 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
============
@@ -215,6 +241,9 @@ 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** __` rasterio plugin to explore COG locally in your web browser.
+Checkout `rio-glui `__ rasterio plugin to explore COG locally in your web browser.
diff --git a/rio_cogeo/cogeo.py b/rio_cogeo/cogeo.py
index 9da270d..d9e73d2 100644
--- 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()
)
diff --git a/rio_cogeo/scripts/cli.py b/rio_cogeo/scripts/cli.py
index 361fb84..1889570 100644
--- a/rio_cogeo/scripts/cli.py
+++ b/rio_cogeo/scripts/cli.py
@@ -1,12 +1,13 @@
"""Rio_cogeo.scripts.cli."""
import os
+import warnings
import click
import numpy
from rasterio.rio import options
-from rasterio.enums import Resampling
+from rasterio.enums import Resampling as ResamplingEnums
from rio_cogeo.cogeo import cog_translate, cog_validate
from rio_cogeo.profiles import cog_profiles
@@ -76,33 +77,53 @@ 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",
- help="Resampling algorithm.",
+ help="Overview creation resampling algorithm.",
type=click.Choice(
- [it.name for it in Resampling if it.value in [0, 1, 2, 3, 4, 5, 6, 7]]
+ [it.name for it in ResamplingEnums if it.value in [0, 1, 2, 3, 4, 5, 6, 7]]
),
default="nearest",
)
@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(
+ "--resampling",
+ "-r",
+ help="Resampling algorithm.",
+ type=click.Choice(
+ [it.name for it in ResamplingEnums if it.value in [0, 1, 2, 3, 4, 5, 6, 7]]
+ ),
+ default="nearest",
)
@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 +135,20 @@ def create(
overview_level,
overview_resampling,
overview_blocksize,
+ web_optimized,
+ latitude_adjustment,
+ resampling,
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 +169,9 @@ def create(
add_mask,
overview_level,
overview_resampling,
+ web_optimized,
+ latitude_adjustment,
+ resampling,
config,
quiet,
)
diff --git a/rio_cogeo/utils.py b/rio_cogeo/utils.py
index 19eed74..8c3aeb1 100644
--- a/rio_cogeo/utils.py
+++ b/rio_cogeo/utils.py
@@ -1,8 +1,97 @@
"""rio_cogeo.utils: Utility functions."""
+import math
+
+from rasterio.warp import calculate_default_transform
from rasterio.enums import MaskFlags, ColorInterp
+def _meters_per_pixel(zoom, lat=0.0, tilesize=256):
+ """
+ Return the pixel resolution for a given mercator tile zoom and lattitude.
+
+ Parameters
+ ----------
+ zoom: int
+ Mercator zoom level
+ lat: float, optional
+ Latitude in decimal degree (default: 0)
+ tilesize: int, optional
+ Mercator tile size (default: 256).
+
+ Returns
+ -------
+ Pixel resolution in meters
+
+ """
+ return (math.cos(lat * math.pi / 180.0) * 2 * math.pi * 6378137) / (
+ tilesize * 2 ** zoom
+ )
+
+
+def zoom_for_pixelsize(pixel_size, max_z=24, tilesize=256):
+ """
+ Get mercator zoom level corresponding to a pixel resolution.
+
+ Freely adapted from
+ https://github.com/OSGeo/gdal/blob/b0dfc591929ebdbccd8a0557510c5efdb893b852/gdal/swig/python/scripts/gdal2tiles.py#L294
+
+ Parameters
+ ----------
+ pixel_size: float
+ Pixel size
+ max_z: int, optional (default: 24)
+ Max mercator zoom level allowed
+ tilesize: int, optional
+ Mercator tile size (default: 256).
+
+ Returns
+ -------
+ Mercator zoom level corresponding to the pixel resolution
+
+ """
+ for z in range(max_z):
+ if pixel_size > _meters_per_pixel(z, 0, tilesize=tilesize):
+ return max(0, z - 1) # We don't want to scale up
+
+ return max_z - 1
+
+
+def get_max_zoom(src_dst, lat=0.0, tilesize=256):
+ """
+ Calculate raster max zoom level.
+
+ Parameters
+ ----------
+ src: rasterio.io.DatasetReader
+ Rasterio io.DatasetReader object
+ lat: float, optional
+ Center latitude of the dataset. This is only needed in case you want to
+ apply latitude correction factor to ensure consitent maximum zoom level
+ (default: 0.0).
+ tilesize: int, optional
+ Mercator tile size (default: 256).
+
+ Returns
+ -------
+ max_zoom: int
+ Max zoom level.
+
+ """
+ dst_affine, w, h = calculate_default_transform(
+ src_dst.crs, "epsg:3857", src_dst.width, src_dst.height, *src_dst.bounds
+ )
+
+ native_resolution = max(abs(dst_affine[0]), abs(dst_affine[4]))
+
+ # Correction factor for web-mercator projection latitude distortion
+ latitude_correction_factor = math.cos(math.radians(lat))
+ corrected_resolution = native_resolution * latitude_correction_factor
+
+ max_zoom = zoom_for_pixelsize(corrected_resolution, tilesize=tilesize)
+ return max_zoom
+
+
def get_maximum_overview_level(src_dst, minsize=512):
"""
Calculate the maximum overview level.
@@ -20,13 +109,9 @@ def get_maximum_overview_level(src_dst, minsize=512):
overview level.
"""
- width = src_dst.width
- height = src_dst.height
-
nlevel = 0
overview = 1
-
- while min(width // overview, height // overview) > minsize:
+ while min(src_dst.width // overview, src_dst.height // overview) > minsize:
overview *= 2
nlevel += 1
diff --git a/setup.py b/setup.py
index e2c5cdd..08bb985 100644
--- a/setup.py
+++ b/setup.py
@@ -1,4 +1,6 @@
"""Setup."""
+
+import sys
from setuptools import setup, find_packages
with open("rio_cogeo/__init__.py") as f:
@@ -14,13 +16,23 @@
readme = f.read()
# Runtime requirements.
-inst_reqs = ["click", "rasterio[s3]>=1.0.9", "numpy~=1.15"]
+inst_reqs = [
+ "click",
+ "rasterio[s3]>=1.0.9",
+ "numpy~=1.15",
+ "supermercado",
+ "mercantile",
+]
extra_reqs = {
- "test": ["pytest", "pytest-cov"],
- "dev": ["pytest", "pytest-cov", "pre-commit"],
+ "test": ["pytest", "pytest-cov", "rio-tiler"],
+ "dev": ["pytest", "pytest-cov", "rio-tiler", "pre-commit"],
}
+if sys.version_info >= (3, 6):
+ extra_reqs["test"].append("cogdumper")
+ extra_reqs["dev"].append("cogdumper")
+
setup(
name="rio-cogeo",
version=version,
@@ -30,14 +42,15 @@
"Intended Audience :: Information Technology",
"Intended Audience :: Science/Research",
"License :: OSI Approved :: BSD License",
- "Programming Language :: Python :: 3.6",
"Programming Language :: Python :: 2.7",
+ "Programming Language :: Python :: 3.6",
+ "Programming Language :: Python :: 3.7",
"Topic :: Scientific/Engineering :: GIS",
],
keywords="COGEO CloudOptimized Geotiff rasterio",
author=u"Vincent Sarago",
- author_email="vincent@mapbox.com",
- url="https://github.com/mapbox/rio-cogeo",
+ author_email="vincent@developmentseed.com",
+ url="https://github.com/cogeotiff/rio-cogeo",
license="BSD-3",
packages=find_packages(exclude=["ez_setup", "examples", "tests"]),
install_requires=inst_reqs,
diff --git a/tests/fixtures/image_north.tif b/tests/fixtures/image_north.tif
new file mode 100644
index 0000000..ae28013
Binary files /dev/null and b/tests/fixtures/image_north.tif differ
diff --git a/tests/fixtures/image_web.tif b/tests/fixtures/image_web.tif
new file mode 100644
index 0000000..f76597a
Binary files /dev/null and b/tests/fixtures/image_web.tif differ
diff --git a/tests/test_cli.py b/tests/test_cli.py
index 0c00104..1394157 100644
--- a/tests/test_cli.py
+++ b/tests/test_cli.py
@@ -279,6 +279,32 @@ def test_cogeo_overviewTilesize(monkeypatch):
assert src.block_shapes[0] == (64, 64)
+def test_cogeo_web():
+ """Should work as expected."""
+ runner = CliRunner()
+ with runner.isolated_filesystem():
+ with pytest.warns(UserWarning):
+ result = runner.invoke(
+ cogeo,
+ ["create", raster_path_rgb, "output.tif", "--latitude-adjustment"],
+ )
+ assert not result.exception
+ assert result.exit_code == 0
+
+ with pytest.warns(UserWarning):
+ result = runner.invoke(
+ cogeo, ["create", raster_path_rgb, "output.tif", "--global-maxzoom"]
+ )
+ assert not result.exception
+ assert result.exit_code == 0
+
+ result = runner.invoke(
+ cogeo, ["create", raster_path_rgb, "output.tif", "--web-optimized"]
+ )
+ assert not result.exception
+ assert result.exit_code == 0
+
+
def test_cogeo_validgdalBlockOption():
"""Should work as expected."""
runner = CliRunner()
@@ -401,9 +427,6 @@ def test_cogeo_validate():
result = runner.invoke(
cogeo, ["create", raster_path_rgb, "output.tif", "--quiet"]
)
- with rasterio.open("output.tif") as src_dst:
- print(src_dst.meta)
-
result = runner.invoke(cogeo, ["validate", "output.tif"])
assert "is a valid cloud optimized GeoTIFF" in result.output
assert not result.exception
diff --git a/tests/test_cogeo.py b/tests/test_cogeo.py
index a79af55..eeae9aa 100644
--- a/tests/test_cogeo.py
+++ b/tests/test_cogeo.py
@@ -33,7 +33,6 @@
os.path.dirname(__file__), "fixtures", "image_rgb_mask.tif"
)
-
jpeg_profile = cog_profiles.get("jpeg")
jpeg_profile.update({"blockxsize": 64, "blockysize": 64})
webp_profile = cog_profiles.get("webp")
@@ -76,6 +75,8 @@ def test_cog_translate_valid():
)
with rasterio.open("cogeo.tif") as src:
assert has_mask_band(src)
+ with rasterio.open("cogeo.tif", OVERVIEW_LEVEL=1) as src:
+ assert src.block_shapes[0] == (64, 64)
def test_cog_translate_NodataLossyWarning():
@@ -92,6 +93,20 @@ def test_cog_translate_NodataLossyWarning():
assert not has_mask_band(src)
+def test_cog_translate_optionWarnings():
+ """Should work as expected but warns about invalid options."""
+ runner = CliRunner()
+ with runner.isolated_filesystem():
+ with pytest.warns(UserWarning):
+ cog_translate(
+ raster_path_rgb, "cogeo.tif", jpeg_profile, nodata=0, quiet=True
+ )
+ with rasterio.open("cogeo.tif") as src:
+ assert src.nodata == 0
+ assert src.compression.value == "JPEG"
+ assert not has_mask_band(src)
+
+
def test_cog_translate_NodataMask():
"""Should work as expected (create cogeo and translate nodata to mask)."""
runner = CliRunner()
@@ -109,14 +124,6 @@ def test_cog_translate_NodataMask():
assert has_mask_band(src)
assert not src.dataset_mask().all()
- cog_translate(
- raster_path_nodata, "cogeo.tif", deflate_profile, add_mask=True, quiet=True
- )
- with rasterio.open("cogeo.tif") as src:
- assert src.nodata is None
- assert has_mask_band(src)
- assert not src.dataset_mask().all()
-
def test_cog_translate_validRaw():
"""Should work as expected (create cogeo file)."""
diff --git a/tests/test_web.py b/tests/test_web.py
new file mode 100644
index 0000000..1a1f727
--- /dev/null
+++ b/tests/test_web.py
@@ -0,0 +1,250 @@
+"""tests rio_cogeo web."""
+
+import os
+import sys
+import struct
+
+import pytest
+from click.testing import CliRunner
+
+import numpy
+
+import mercantile
+import rasterio
+from rasterio.warp import transform_bounds
+
+from rio_cogeo.utils import get_max_zoom
+from rio_cogeo.cogeo import cog_translate
+from rio_cogeo.profiles import cog_profiles
+
+from rio_tiler.utils import tile_read
+
+
+raster_path_web = os.path.join(os.path.dirname(__file__), "fixtures", "image_web.tif")
+raster_path_north = os.path.join(
+ os.path.dirname(__file__), "fixtures", "image_north.tif"
+)
+
+
+def test_cog_translate_webZooms():
+ """
+ Test Web-Optimized COG.
+
+ - Test COG size is a multiple of 256 (mercator tile size)
+ - Test COG bounds are aligned with mercator grid at max zoom
+ - Test high resolution internal tiles are equal to mercator tile using
+ cogdumper and rio-tiler
+ - Test overview internal tiles are equal to mercator tile using
+ cogdumper and rio-tiler
+ """
+
+ runner = CliRunner()
+ with runner.isolated_filesystem():
+ web_profile = cog_profiles.get("raw")
+ web_profile.update({"blockxsize": 256, "blockysize": 256})
+ config = dict(GDAL_TIFF_OVR_BLOCKSIZE="128")
+
+ cog_translate(
+ raster_path_north,
+ "cogeo.tif",
+ web_profile,
+ quiet=True,
+ web_optimized=True,
+ config=config,
+ )
+ with rasterio.open("cogeo.tif") as out_dst:
+ assert get_max_zoom(out_dst) == 8
+
+ cog_translate(
+ raster_path_north,
+ "cogeo.tif",
+ web_profile,
+ quiet=True,
+ web_optimized=True,
+ latitude_adjustment=False,
+ config=config,
+ )
+ with rasterio.open("cogeo.tif") as out_dst:
+ assert get_max_zoom(out_dst) == 10
+
+
+def test_cog_translate_web():
+ """
+ Test Web-Optimized COG.
+
+ - Test COG size is a multiple of 256 (mercator tile size)
+ - Test COG bounds are aligned with mercator grid at max zoom
+ """
+ runner = CliRunner()
+ with runner.isolated_filesystem():
+
+ web_profile = cog_profiles.get("raw")
+ web_profile.update({"blockxsize": 256, "blockysize": 256})
+ config = dict(GDAL_TIFF_OVR_BLOCKSIZE="128")
+
+ cog_translate(
+ raster_path_web,
+ "cogeo.tif",
+ web_profile,
+ quiet=True,
+ web_optimized=True,
+ config=config,
+ )
+ with rasterio.open(raster_path_web) as src_dst:
+ with rasterio.open("cogeo.tif") as out_dst:
+ blocks = list(set(out_dst.block_shapes))
+ assert len(blocks) == 1
+ ts = blocks[0][0]
+ assert not out_dst.width % ts
+ assert not out_dst.height % ts
+ max_zoom = get_max_zoom(out_dst)
+
+ bounds = list(
+ transform_bounds(
+ *[src_dst.crs, "epsg:4326"] + list(src_dst.bounds),
+ densify_pts=21
+ )
+ )
+
+ leftTile = mercantile.tile(bounds[0], bounds[3], max_zoom)
+ tbounds = mercantile.xy_bounds(leftTile)
+ west, north = tbounds.left, tbounds.top
+ assert out_dst.transform.xoff == west
+ assert out_dst.transform.yoff == north
+
+ rightTile = mercantile.tile(bounds[2], bounds[1], max_zoom)
+ tbounds = mercantile.xy_bounds(rightTile)
+ east, south = tbounds.right, tbounds.bottom
+
+ lrx = round(
+ out_dst.transform.xoff + out_dst.transform.a * out_dst.width, 6
+ )
+ lry = round(
+ out_dst.transform.yoff + out_dst.transform.e * out_dst.height, 6
+ )
+ assert lrx == round(east, 6)
+ assert lry == round(south, 6)
+
+
+@pytest.mark.skipif(sys.version_info < (3, 6), reason="requires python3.6 or higher")
+def test_cog_translate_Internal():
+ """
+ Test Web-Optimized COG.
+
+ - Test COG size is a multiple of 256 (mercator tile size)
+ - Test COG bounds are aligned with mercator grid at max zoom
+ - Test high resolution internal tiles are equal to mercator tile using
+ cogdumper and rio-tiler
+ - Test overview internal tiles are equal to mercator tile using
+ cogdumper and rio-tiler
+ """
+ from cogdumper.cog_tiles import COGTiff
+ from cogdumper.filedumper import Reader as FileReader
+
+ runner = CliRunner()
+ with runner.isolated_filesystem():
+
+ web_profile = cog_profiles.get("raw")
+ web_profile.update({"blockxsize": 256, "blockysize": 256})
+ config = dict(GDAL_TIFF_OVR_BLOCKSIZE="128")
+
+ cog_translate(
+ raster_path_web,
+ "cogeo.tif",
+ web_profile,
+ quiet=True,
+ web_optimized=True,
+ config=config,
+ )
+ with rasterio.open(raster_path_web) as src_dst:
+ with rasterio.open("cogeo.tif") as out_dst:
+ blocks = list(set(out_dst.block_shapes))
+ assert len(blocks) == 1
+ ts = blocks[0][0]
+ assert not out_dst.width % ts
+ assert not out_dst.height % ts
+ max_zoom = get_max_zoom(out_dst)
+
+ bounds = list(
+ transform_bounds(
+ *[src_dst.crs, "epsg:4326"] + list(src_dst.bounds),
+ densify_pts=21
+ )
+ )
+
+ leftTile = mercantile.tile(bounds[0], bounds[3], max_zoom)
+ tbounds = mercantile.xy_bounds(leftTile)
+ west, north = tbounds.left, tbounds.top
+ assert out_dst.transform.xoff == west
+ assert out_dst.transform.yoff == north
+
+ rightTile = mercantile.tile(bounds[2], bounds[1], max_zoom)
+ tbounds = mercantile.xy_bounds(rightTile)
+ east, south = tbounds.right, tbounds.bottom
+
+ lrx = round(
+ out_dst.transform.xoff + out_dst.transform.a * out_dst.width, 6
+ )
+ lry = round(
+ out_dst.transform.yoff + out_dst.transform.e * out_dst.height, 6
+ )
+ assert lrx == round(east, 6)
+ assert lry == round(south, 6)
+
+ with open("cogeo.tif", "rb") as out_body:
+ reader = FileReader(out_body)
+ cog = COGTiff(reader.read)
+
+ # High resolution
+ # Top Left tile
+ mime_type, tile = cog.get_tile(0, 0, 0)
+ tile_length = 256 * 256 * 3
+ t = struct.unpack_from("{}b".format(tile_length), tile)
+ arr = numpy.array(t).reshape(256, 256, 3).astype(numpy.uint8)
+ arr = numpy.transpose(arr, [2, 0, 1])
+
+ tbounds = mercantile.xy_bounds(leftTile)
+ data, mask = tile_read(
+ "cogeo.tif", tbounds, 256, resampling_method="nearest"
+ )
+ numpy.testing.assert_array_equal(data, arr)
+
+ # Bottom right tile
+ mime_type, tile = cog.get_tile(4, 3, 0)
+ tile_length = 256 * 256 * 3
+ t = struct.unpack_from("{}b".format(tile_length), tile)
+ arr = numpy.array(t).reshape(256, 256, 3).astype(numpy.uint8)
+ arr = numpy.transpose(arr, [2, 0, 1])
+
+ tbounds = mercantile.xy_bounds(rightTile)
+ data, mask = tile_read(
+ "cogeo.tif", tbounds, 256, resampling_method="nearest"
+ )
+ numpy.testing.assert_array_equal(data, arr)
+
+ # Low resolution (overview 1)
+ # Top Left tile
+ # NOTE: overview internal tile size is 128px
+ # We need to stack two internal tiles to compare with
+ # the 256px mercator tile fetched by rio-tiler
+ # ref: https://github.com/cogeotiff/rio-cogeo/issues/60
+ mime_type, tile = cog.get_tile(1, 0, 1)
+ tile_length = 128 * 128 * 3
+ t = struct.unpack_from("{}b".format(tile_length), tile)
+ arr1 = numpy.array(t).reshape(128, 128, 3).astype(numpy.uint8)
+ arr1 = numpy.transpose(arr1, [2, 0, 1])
+
+ mime_type, tile = cog.get_tile(2, 0, 1)
+ tile_length = 128 * 128 * 3
+ t = struct.unpack_from("{}b".format(tile_length), tile)
+ arr2 = numpy.array(t).reshape(128, 128, 3).astype(numpy.uint8)
+ arr2 = numpy.transpose(arr2, [2, 0, 1])
+ arr = numpy.dstack((arr1, arr2))
+
+ lowTile = mercantile.Tile(118594, 60034, 17)
+ tbounds = mercantile.xy_bounds(lowTile)
+ data, mask = tile_read(
+ "cogeo.tif", tbounds, 256, resampling_method="nearest"
+ )
+ data = data[:, 128:, :]
+ numpy.testing.assert_array_equal(data, arr)
diff --git a/tox.ini b/tox.ini
index c315ff8..987a910 100644
--- a/tox.ini
+++ b/tox.ini
@@ -1,5 +1,5 @@
[tox]
-envlist = py37,py36,py27
+envlist = py27,py37,py36
[testenv]
extras = test