diff --git a/.github/workflows/actions.yml b/.github/workflows/actions.yml index 478f353..5070e7a 100644 --- a/.github/workflows/actions.yml +++ b/.github/workflows/actions.yml @@ -49,7 +49,8 @@ jobs: - run: python .github/check_examples.py - run: cd js && npm ci && npm test - run: cd js && npm run biome-check - - run: cd python && python -m unittest test/test_* + - run: cd python/pmtiles && python -m unittest test/test_* + - run: cd python/rio-pmtiles && pip install -r requirements-dev.txt && pip install -e . && PYTHONPATH=. pytest tests - run: cd cpp && make - run: cd serverless/cloudflare && npm ci && npm test - run: cd serverless/vtfilter && npm ci && npm test diff --git a/REUSE.toml b/REUSE.toml index 692542d..2bff79c 100644 --- a/REUSE.toml +++ b/REUSE.toml @@ -1,10 +1,15 @@ version = 1 [[annotations]] -path = ["app/**","cpp/**","js/**","openlayers/**","python/**","serverless/**"] +path = ["app/**","cpp/**","js/**","openlayers/**","python/pmtiles/**","serverless/**"] SPDX-FileCopyrightText = "2021 and later, Protomaps LLC and contributors" SPDX-License-Identifier = "BSD-3-Clause" +[[annotations]] +path = ["python/rio-pmtiles/**"] +SPDX-FileCopyrightText = "2025 and later, Protomaps LLC and contributors; 2024-2021 Mapbox" +SPDX-License-Identifier = "MIT" + [[annotations]] path = "spec/**/*.md" SPDX-FileCopyrightText = "2021 and later, Protomaps LLC and contributors" diff --git a/python/README.md b/python/pmtiles/README.md similarity index 100% rename from python/README.md rename to python/pmtiles/README.md diff --git a/python/bin/pmtiles-convert b/python/pmtiles/bin/pmtiles-convert similarity index 100% rename from python/bin/pmtiles-convert rename to python/pmtiles/bin/pmtiles-convert diff --git a/python/bin/pmtiles-serve b/python/pmtiles/bin/pmtiles-serve similarity index 100% rename from python/bin/pmtiles-serve rename to python/pmtiles/bin/pmtiles-serve diff --git a/python/bin/pmtiles-show b/python/pmtiles/bin/pmtiles-show similarity index 100% rename from python/bin/pmtiles-show rename to python/pmtiles/bin/pmtiles-show diff --git a/python/examples/create_raster_example.py b/python/pmtiles/examples/create_raster_example.py similarity index 100% rename from python/examples/create_raster_example.py rename to python/pmtiles/examples/create_raster_example.py diff --git a/python/pmtiles/__init__.py b/python/pmtiles/pmtiles/__init__.py similarity index 100% rename from python/pmtiles/__init__.py rename to python/pmtiles/pmtiles/__init__.py diff --git a/python/pmtiles/convert.py b/python/pmtiles/pmtiles/convert.py similarity index 100% rename from python/pmtiles/convert.py rename to python/pmtiles/pmtiles/convert.py diff --git a/python/pmtiles/reader.py b/python/pmtiles/pmtiles/reader.py similarity index 100% rename from python/pmtiles/reader.py rename to python/pmtiles/pmtiles/reader.py diff --git a/python/pmtiles/tile.py b/python/pmtiles/pmtiles/tile.py similarity index 100% rename from python/pmtiles/tile.py rename to python/pmtiles/pmtiles/tile.py diff --git a/python/pmtiles/v2.py b/python/pmtiles/pmtiles/v2.py similarity index 100% rename from python/pmtiles/v2.py rename to python/pmtiles/pmtiles/v2.py diff --git a/python/pmtiles/writer.py b/python/pmtiles/pmtiles/writer.py similarity index 100% rename from python/pmtiles/writer.py rename to python/pmtiles/pmtiles/writer.py diff --git a/python/setup.py b/python/pmtiles/setup.py similarity index 100% rename from python/setup.py rename to python/pmtiles/setup.py diff --git a/python/test/__init__.py b/python/pmtiles/test/__init__.py similarity index 100% rename from python/test/__init__.py rename to python/pmtiles/test/__init__.py diff --git a/python/test/test_convert.py b/python/pmtiles/test/test_convert.py similarity index 100% rename from python/test/test_convert.py rename to python/pmtiles/test/test_convert.py diff --git a/python/test/test_reader_writer.py b/python/pmtiles/test/test_reader_writer.py similarity index 100% rename from python/test/test_reader_writer.py rename to python/pmtiles/test/test_reader_writer.py diff --git a/python/test/test_tile.py b/python/pmtiles/test/test_tile.py similarity index 100% rename from python/test/test_tile.py rename to python/pmtiles/test/test_tile.py diff --git a/python/rio-pmtiles/CHANGELOG.md b/python/rio-pmtiles/CHANGELOG.md new file mode 100644 index 0000000..07a21d6 --- /dev/null +++ b/python/rio-pmtiles/CHANGELOG.md @@ -0,0 +1,17 @@ +0.0.6 +------ + +This is a port of the original rio-mbtiles developed by @sgillies at https://github.com/mapbox/rio-mbtiles. + +Differences are: + +* output clustered, compressed PMTiles archives instead of MBTiles +* require Python 3.7 or above +* remove `multiprocessing` implementation in favor of `concurrent.futures`; remove `--implementation` option +* remove `--covers` quadkey feature +* replace progress estimation with exact counts; add `pyroaring` dependency +* remove `--image-dump` feature +* remove `--append/--overwrite` since PMTiles is not a database +* update dependencies + +Otherwise, the behavior should be the same as `rio-mbtiles` v1.6.0. \ No newline at end of file diff --git a/python/rio-pmtiles/LICENSE b/python/rio-pmtiles/LICENSE new file mode 100644 index 0000000..3113611 --- /dev/null +++ b/python/rio-pmtiles/LICENSE @@ -0,0 +1,22 @@ +The MIT License (MIT) + +Copyright © 2025 Protomaps LLC, 2014-2021 Mapbox + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. + diff --git a/python/rio-pmtiles/README.rst b/python/rio-pmtiles/README.rst new file mode 100644 index 0000000..bc2fd85 --- /dev/null +++ b/python/rio-pmtiles/README.rst @@ -0,0 +1,113 @@ +rio-pmtiles +=========== + +A plugin for the +`Rasterio CLI `__ +that exports a raster dataset to the PMTiles (version 3) format. Features +include automatic reprojection and concurrent tile generation. + +This is derived from the original rio-mbtiles CLI by Sean Gillies at Mapbox. See CHANGELOG.md for differences. + +Usage +----- + +.. code-block:: console + + Usage: rio pmtiles [OPTIONS] INPUT [OUTPUT] + + Export a dataset to PMTiles (version 3). + + The input dataset may have any coordinate reference system. It must have + at least three bands, which will be become the red, blue, and green bands + of the output image tiles. + + An optional fourth alpha band may be copied to the output tiles by using + the --rgba option in combination with the PNG or WEBP formats. This option + requires that the input dataset has at least 4 bands. + + The default quality for JPEG and WEBP output (possible range: 10-100) is + 75. This value can be changed with the use of the QUALITY creation option, + e.g. `--co QUALITY=90`. The default zlib compression level for PNG output + (possible range: 1-9) is 6. This value can be changed like `--co + ZLEVEL=8`. Lossless WEBP can be chosen with `--co LOSSLESS=TRUE`. + + If no zoom levels are specified, the defaults are the zoom levels nearest + to the one at which one tile may contain the entire source dataset. + + If a title or description for the output file are not provided, they will + be taken from the input dataset's filename. + + This command is suited for small to medium (~1 GB) sized sources. + + Python package: rio-pmtiles (https://github.com/protomaps/PMTiles). + + Options: + -o, --output PATH Path to output file (optional alternative to + a positional arg). + + --title TEXT PMTiles dataset title. + --description TEXT PMTiles dataset description. + --overlay Export as an overlay (the default). + --baselayer Export as a base layer. + -f, --format [JPEG|PNG|WEBP] Tile image format. + --tile-size INTEGER Width and height of individual square tiles + to create. [default: 256] + + --zoom-levels MIN..MAX A min...max range of export zoom levels. The + default zoom level is the one at which the + dataset is contained within a single tile. + + -j INTEGER Number of workers (default: number of + computer's processors). + + --src-nodata FLOAT Manually override source nodata + --dst-nodata FLOAT Manually override destination nodata + --resampling [nearest|bilinear|cubic|cubic_spline|lanczos|average|mode|gauss|max|min|med|q1|q3|rms] + Resampling method to use. [default: + nearest] + + --version Show the version and exit. + --rgba Select RGBA output. For PNG or WEBP only. + + -#, --progress-bar Display progress bar. + --cutline PATH Path to a GeoJSON FeatureCollection to be + used as a cutline. Only source pixels within + the cutline features will be exported. + + --oo NAME=VALUE Format driver-specific options to be used + when accessing the input dataset. See the + GDAL format driver documentation for more + information. + + --co, --profile NAME=VALUE Driver specific creation options. See the + documentation for the selected output driver + for more information. + + --wo NAME=VALUE See the GDAL warp options documentation for + more information. + + --exclude-empty-tiles / --include-empty-tiles + Whether to exclude or include empty tiles + from the output. + + --help Show this message and exit. + +Performance +----------- + +The rio-pmtiles command is suited for small to medium (~1 GB) raster sources. +On a MacBook Pro M1, the 1:10M scale Natural Earth raster +(a 21,600 x 10,800 pixel, 700 MB TIFF) exports to PMTiles (levels 1 through 5) +in 15 seconds. + +.. code-block:: console + + $ time GDAL_CACHEMAX=256 rio pmtiles NE1_HR_LC.tif \ + > -o ne.pmtiles --zoom-levels 1..5 -j 4 + + 14.87s user 10.40s system 258% cpu 9.787 total + +Installation +------------ + +``pip install rio-pmtiles`` diff --git a/python/rio-pmtiles/requirements-dev.txt b/python/rio-pmtiles/requirements-dev.txt new file mode 100644 index 0000000..00cbcf7 --- /dev/null +++ b/python/rio-pmtiles/requirements-dev.txt @@ -0,0 +1,14 @@ +affine==2.2.1 +cligj==0.5.0 +coveralls==1.5.1 +mercantile==1.1.6 +pmtiles~=3.0 +pyroaring~=1.0 +pytest~=8.3.5 +pytest-cov==2.5.1 +coverage==4.5.1 +rasterio~=1.0 +shapely~=2.0.0 +supermercado==0.2.0 +snuggs==1.4.2 +tqdm==4.48.2 diff --git a/python/rio-pmtiles/rio_pmtiles/__init__.py b/python/rio-pmtiles/rio_pmtiles/__init__.py new file mode 100644 index 0000000..7b74956 --- /dev/null +++ b/python/rio-pmtiles/rio_pmtiles/__init__.py @@ -0,0 +1,5 @@ +"""rio-pmtiles package""" + +import sys + +__version__ = "0.0.6" diff --git a/python/rio-pmtiles/rio_pmtiles/scripts/__init__.py b/python/rio-pmtiles/rio_pmtiles/scripts/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/python/rio-pmtiles/rio_pmtiles/scripts/cli.py b/python/rio-pmtiles/rio_pmtiles/scripts/cli.py new file mode 100644 index 0000000..56e1c38 --- /dev/null +++ b/python/rio-pmtiles/rio_pmtiles/scripts/cli.py @@ -0,0 +1,463 @@ +"""rio-pmtiles CLI""" + +import functools +import logging +import math +import os +import sqlite3 +import sys + +import click +from cligj.features import iter_features +import concurrent.futures +import gzip +from itertools import islice +import json +import mercantile +from pyroaring import BitMap64 +import rasterio +from rasterio.enums import Resampling +from rasterio.rio.options import creation_options, output_opt, _cb_key_val +from rasterio.warp import transform, transform_geom +import shapely.affinity +from shapely.geometry import mapping, shape +from shapely.ops import unary_union +import shapely.wkt +import supermercado.burntiles +from tqdm import tqdm +from pmtiles.tile import TileType, Entry, zxy_to_tileid, Compression, serialize_header, tileid_to_zxy +from pmtiles.writer import Writer, optimize_directories + +from rio_pmtiles import __version__ as rio_pmtiles_version +from rio_pmtiles.worker import init_worker, process_tile + + +DEFAULT_NUM_WORKERS = None +RESAMPLING_METHODS = [method.name for method in Resampling] +TILES_CRS = "EPSG:3857" + +log = logging.getLogger(__name__) + + +def resolve_inout( + input=None, output=None, files=None, num_inputs=None +): + """Resolves inputs and outputs from standard args and options. + + Parameters + ---------- + input : str + A single input filename, optional. + output : str + A single output filename, optional. + files : str + A sequence of filenames in which the last is the output filename. + num_inputs : int + Raise exceptions if the number of resolved input files is higher + or lower than this number. + + Returns + ------- + tuple (str, list of str) + The resolved output filename and input filenames as a tuple of + length 2. + + If provided, the output file may be overwritten. An output + file extracted from files will not be overwritten unless + overwrite is True. + + Raises + ------ + click.BadParameter + + """ + resolved_output = output or (files[-1] if files else None) + resolved_inputs = ( + [input] + if input + else [] + list(files[: -1 if not output else None]) + if files + else [] + ) + + if num_inputs is not None: + if len(resolved_inputs) < num_inputs: + raise click.BadParameter("Insufficient inputs") + elif len(resolved_inputs) > num_inputs: + raise click.BadParameter("Too many inputs") + + return resolved_output, resolved_inputs + + +def extract_features(ctx, param, value): + if value is not None: + with click.open_file(value, encoding="utf-8") as src: + return list(iter_features(iter(src))) + else: + return None + + +@click.command(short_help="Export a dataset to PMTiles.") +@click.argument( + "files", + nargs=-1, + type=click.Path(resolve_path=True), + required=True, + metavar="INPUT [OUTPUT]", +) +@output_opt +@click.option("--title", help="PMTiles dataset title.") +@click.option("--description", help="PMTiles dataset description.") +@click.option( + "--overlay", + "layer_type", + flag_value="overlay", + default=True, + help="Export as an overlay (the default).", +) +@click.option( + "--baselayer", "layer_type", flag_value="baselayer", help="Export as a base layer." +) +@click.option( + "-f", + "--format", + "img_format", + type=click.Choice(["JPEG", "PNG", "WEBP"]), + default="JPEG", + help="Tile image format.", +) +@click.option( + "--tile-size", + default=256, + show_default=True, + type=int, + help="Width and height of individual square tiles to create.", +) +@click.option( + "--zoom-levels", + default=None, + metavar="MIN..MAX", + help="A min...max range of export zoom levels. " + "The default zoom level " + "is the one at which the dataset is contained within " + "a single tile.", +) +@click.option( + "-j", + "num_workers", + type=int, + default=DEFAULT_NUM_WORKERS, + help="Number of workers (default: number of computer's processors).", +) +@click.option( + "--src-nodata", + default=None, + show_default=True, + type=float, + help="Manually override source nodata", +) +@click.option( + "--dst-nodata", + default=None, + show_default=True, + type=float, + help="Manually override destination nodata", +) +@click.option( + "--resampling", + type=click.Choice(RESAMPLING_METHODS), + default="nearest", + show_default=True, + help="Resampling method to use.", +) +@click.version_option(version=rio_pmtiles_version, message="%(version)s") +@click.option( + "--rgba", default=False, is_flag=True, help="Select RGBA output. For PNG or WEBP only." +) +@click.option( + "--progress-bar", "-#", default=False, is_flag=True, help="Don't display progress bar." +) +@click.option( + "--cutline", + type=click.Path(exists=True), + callback=extract_features, + default=None, + help="Path to a GeoJSON FeatureCollection to be used as a cutline. Only source pixels within the cutline features will be exported.", +) +@click.option( + "--oo", + "open_options", + metavar="NAME=VALUE", + multiple=True, + callback=_cb_key_val, + help="Format driver-specific options to be used when accessing the input dataset. See the GDAL format driver documentation for more information.", +) +@creation_options +@click.option( + "--wo", + "warp_options", + metavar="NAME=VALUE", + multiple=True, + callback=_cb_key_val, + help="See the GDAL warp options documentation for more information.", +) +@click.option( + "--exclude-empty-tiles/--include-empty-tiles", + default=True, + is_flag=True, + help="Whether to exclude or include empty tiles from the output.", +) +@click.pass_context +def pmtiles( + ctx, + files, + output, + title, + description, + layer_type, + img_format, + tile_size, + zoom_levels, + num_workers, + src_nodata, + dst_nodata, + resampling, + rgba, + progress_bar, + cutline, + open_options, + creation_options, + warp_options, + exclude_empty_tiles, +): + """Export a dataset to PMTiles. + + The input dataset may have any coordinate reference system. It must + have at least three bands, which will be become the red, blue, and + green bands of the output image tiles. + + An optional fourth alpha band may be copied to the output tiles by + using the --rgba option in combination with the PNG or WEBP formats. + This option requires that the input dataset has at least 4 bands. + + The default quality for JPEG and WEBP output (possible range: + 10-100) is 75. This value can be changed with the use of the QUALITY + creation option, e.g. `--co QUALITY=90`. The default zlib + compression level for PNG output (possible range: 1-9) is 6. This + value can be changed like `--co ZLEVEL=8`. Lossless WEBP can be + chosen with `--co LOSSLESS=TRUE`. + + If no zoom levels are specified, the defaults are the zoom levels + nearest to the one at which one tile may contain the entire source + dataset. + + If a title or description for the output file are not provided, + they will be taken from the input dataset's filename. + + This command is suited for small to medium (~1 GB) sized sources. + + Python package: rio-pmtiles (https://github.com/protomaps/PMTiles). + + """ + log = logging.getLogger(__name__) + + output, files = resolve_inout( + files=files, output=output, num_inputs=1, + ) + inputfile = files[0] + + with ctx.obj["env"]: + + # Read metadata from the source dataset. + with rasterio.open(inputfile, **open_options) as src: + + if dst_nodata is not None and ( + src_nodata is None and src.profile.get("nodata") is None + ): + raise click.BadParameter( + "--src-nodata must be provided because " "dst-nodata is not None." + ) + base_kwds = {"dst_nodata": dst_nodata, "src_nodata": src_nodata} + + if src_nodata is not None: + base_kwds.update(nodata=src_nodata) + + if dst_nodata is not None: + base_kwds.update(nodata=dst_nodata) + + # Name and description. + title = title or os.path.basename(src.name) + description = description or src.name + + # Compute the geographic bounding box of the dataset. + (west, east), (south, north) = transform( + src.crs, "EPSG:4326", src.bounds[::2], src.bounds[1::2] + ) + + # cutlines must be transformed from CRS84 to src pixel/line + # coordinates and then formatted as WKT. + if cutline is not None: + geoms = [shape(f["geometry"]) for f in cutline] + union = unary_union(geoms) + if union.geom_type not in ("MultiPolygon", "Polygon"): + raise click.ClickException("Unexpected cutline geometry type") + west, south, east, north = union.bounds + cutline_src = shape( + transform_geom("OGC:CRS84", src.crs, mapping(union)) + ) + invtransform = ~src.transform + shapely_matrix = ( + invtransform.a, + invtransform.b, + invtransform.d, + invtransform.e, + invtransform.xoff, + invtransform.yoff, + ) + cutline_rev = shapely.affinity.affine_transform( + cutline_src, shapely_matrix + ) + warp_options["cutline"] = shapely.wkt.dumps(cutline_rev) + + # Resolve the minimum and maximum zoom levels for export. + if zoom_levels: + minzoom, maxzoom = map(int, zoom_levels.split("..")) + else: + zw = int(round(math.log(360.0 / (east - west), 2.0))) + zh = int(round(math.log(170.1022 / (north - south), 2.0))) + minzoom = min(zw, zh) + maxzoom = max(zw, zh) + + log.debug("Zoom range: %d..%d", minzoom, maxzoom) + + if rgba: + if img_format == "JPEG": + raise click.BadParameter( + "RGBA output is not possible with JPEG format." + ) + else: + count = 4 + else: + count = src.count + + # Parameters for creation of tile images. + base_kwds.update( + { + "driver": img_format.upper(), + "dtype": "uint8", + "nodata": 0, + "height": tile_size, + "width": tile_size, + "count": count, + "crs": TILES_CRS, + } + ) + + # Constrain bounds. + EPS = 1.0e-10 + west = max(-180 + EPS, west) + south = max(-85.051129, south) + east = min(180 - EPS, east) + north = min(85.051129, north) + + outfile = open(output, "wb") + outfile.write(b"\x00" * 16384) + entries = [] + + metadata = gzip.compress(json.dumps({'name':title,'type':layer_type,'description':description,'writer':f'rio-pmtiles {rio_pmtiles_version}'}).encode()) + outfile.write(metadata) + + header = {} + header["version"] = 3 + header["root_offset"] = 127 + header["metadata_offset"] = 16384 + header["metadata_length"] = len(metadata) + header["tile_data_offset"] = 16384 + len(metadata) + + if img_format == "JPEG": + header["tile_type"] = TileType.JPEG + elif img_format == "WEBP": + header["tile_type"] = TileType.WEBP + elif img_format == "PNG": + header["tile_type"] = TileType.PNG + + header["tile_compression"] = Compression.NONE + header["internal_compression"] = Compression.GZIP + header["clustered"] = True + header["min_zoom"] = minzoom + header["max_zoom"] = maxzoom + header["min_lon_e7"] = int(west * 10000000) + header["min_lat_e7"] = int(south * 10000000) + header["max_lon_e7"] = int(east * 10000000) + header["max_lat_e7"] = int(north * 10000000) + header["center_zoom"] = minzoom + header["center_lon_e7"] = int((west + east) / 2 * 10000000) + header["center_lat_e7"] = int((south + north) / 2 * 10000000) + + tiles = BitMap64() + + if cutline: + for zk in range(minzoom, maxzoom + 1): + for arr in supermercado.burntiles.burn(cutline, zk): + tiles.add(zxy_to_tileid(arr[2],arr[0],arr[1])) + else: + for tile in mercantile.tiles(west, south, east, north, range(minzoom, maxzoom + 1)): + tiles.add(zxy_to_tileid(tile.z,tile.x,tile.y)) + + def unwrap_tiles(bmap): + for tile_id in bmap: + z, x, y = tileid_to_zxy(tile_id) + yield mercantile.Tile(x,y,z) + + if progress_bar: + pbar = tqdm(total=len(tiles)) + else: + pbar = None + + tile_data_offset = 0 + + """Warp imagery into tiles and write to pmtiles archive. + """ + with concurrent.futures.ProcessPoolExecutor( + max_workers=num_workers, + initializer=init_worker, + initargs=( + inputfile, + base_kwds, + resampling, + open_options, + warp_options, + creation_options, + exclude_empty_tiles, + ), + ) as executor: + for tile, contents in executor.map(process_tile, unwrap_tiles(tiles)): + if pbar is not None: + pbar.update(1) + if contents is None: + log.info("Tile %r is empty and will be skipped", tile) + continue + log.info("Inserting tile: tile=%r", tile) + + entries.append(Entry(zxy_to_tileid(tile.z,tile.x,tile.y), tile_data_offset, len(contents), 1)) + outfile.write(contents) + tile_data_offset += len(contents) + + header["addressed_tiles_count"] = len(entries) + header["tile_entries_count"] = len(entries) + header["tile_contents_count"] = len(entries) + header["tile_data_length"] = tile_data_offset + + root, leaves, num_leaves = optimize_directories(entries, 16384-127) + header["root_length"] = len(root) + if len(leaves) > 0: + outfile.write(leaves) + header["leaf_directory_offset"] = 16384 + len(metadata) + tile_data_offset + header["leaf_directory_length"] = len(leaves) + else: + header["leaf_directory_offset"] = header["tile_data_offset"] + header["leaf_directory_length"] = 0 + + outfile.seek(0) + outfile.write(serialize_header(header)) + outfile.write(root) diff --git a/python/rio-pmtiles/rio_pmtiles/worker.py b/python/rio-pmtiles/rio_pmtiles/worker.py new file mode 100644 index 0000000..bf9cf1b --- /dev/null +++ b/python/rio-pmtiles/rio_pmtiles/worker.py @@ -0,0 +1,138 @@ +"""rio-pmtiles processing worker""" + +import logging +import warnings + +from rasterio.enums import Resampling +from rasterio.io import MemoryFile +from rasterio.transform import from_bounds as transform_from_bounds +from rasterio.warp import reproject, transform_bounds +from rasterio.windows import Window +from rasterio.windows import from_bounds as window_from_bounds +import mercantile +import rasterio + +TILES_CRS = "EPSG:3857" + +log = logging.getLogger(__name__) + + +def init_worker( + path, + profile, + resampling_method, + open_opts=None, + warp_opts=None, + creation_opts=None, + exclude_empties=True, +): + global base_kwds, filename, resampling, open_options, warp_options, creation_options, exclude_empty_tiles + resampling = Resampling[resampling_method] + base_kwds = profile.copy() + filename = path + open_options = open_opts.copy() if open_opts is not None else {} + warp_options = warp_opts.copy() if warp_opts is not None else {} + creation_options = creation_opts.copy() if creation_opts is not None else {} + exclude_empty_tiles = exclude_empties + + +def process_tile(tile): + """Process a single PMTiles tile + + Parameters + ---------- + tile : mercantile.Tile + warp_options : Mapping + GDAL warp options as keyword arguments. + + Returns + ------- + + tile : mercantile.Tile + The input tile. + bytes : bytearray + Image bytes corresponding to the tile. + + """ + global base_kwds, resampling, filename, open_options, warp_options, creation_options, exclude_empty_tiles + + with rasterio.open(filename, **open_options) as src: + + bbox = mercantile.xy_bounds(tile) + + kwds = base_kwds.copy() + kwds.update(**creation_options) + kwds["transform"] = transform_from_bounds( + bbox.left, bbox.bottom, bbox.right, bbox.top, kwds["width"], kwds["height"] + ) + src_nodata = kwds.pop("src_nodata", None) + dst_nodata = kwds.pop("dst_nodata", None) + + src_alpha = None + dst_alpha = None + bindexes = None + + if kwds["count"] == 4: + bindexes = [1, 2, 3] + dst_alpha = 4 + + if src.count == 4: + src_alpha = 4 + else: + kwds["count"] = 4 + else: + bindexes = list(range(1, kwds["count"] + 1)) + + warnings.simplefilter("ignore") + + log.info("Reprojecting tile: tile=%r", tile) + + with MemoryFile() as memfile: + + with memfile.open(**kwds) as tmp: + + # determine window of source raster corresponding to the tile + # image, with small buffer at edges + try: + west, south, east, north = transform_bounds( + TILES_CRS, src.crs, bbox.left, bbox.bottom, bbox.right, bbox.top + ) + tile_window = window_from_bounds( + west, south, east, north, transform=src.transform + ) + adjusted_tile_window = Window( + tile_window.col_off - 1, + tile_window.row_off - 1, + tile_window.width + 2, + tile_window.height + 2, + ) + tile_window = adjusted_tile_window.round_offsets().round_shape() + + # if no data in window, skip processing the tile + if ( + exclude_empty_tiles + and not src.read_masks(1, window=tile_window).any() + ): + return tile, None + + except ValueError: + log.info( + "Tile %r will not be skipped, even if empty. This is harmless.", + tile, + ) + + num_threads = int(warp_options.pop("num_threads", 2)) + + reproject( + rasterio.band(src, bindexes), + rasterio.band(tmp, bindexes), + src_nodata=src_nodata, + dst_nodata=dst_nodata, + src_alpha=src_alpha, + dst_alpha=dst_alpha, + num_threads=num_threads, + resampling=resampling, + **warp_options + ) + + return tile, memfile.read() diff --git a/python/rio-pmtiles/setup.py b/python/rio-pmtiles/setup.py new file mode 100644 index 0000000..b5b6476 --- /dev/null +++ b/python/rio-pmtiles/setup.py @@ -0,0 +1,50 @@ +from codecs import open as codecs_open +from setuptools import setup, find_packages + + +# Parse the version from the rio_pmtiles module. +with open('rio_pmtiles/__init__.py') as f: + for line in f: + if line.find("__version__") >= 0: + version = line.split("=")[1].strip() + version = version.strip('"') + version = version.strip("'") + break + +# Get the long description from the relevant file +with codecs_open('README.rst', encoding='utf-8') as f: + long_description = f.read() + + +setup( + name="rio-pmtiles", + version=version, + description=u"A Rasterio plugin command that exports PMTiles", + long_description=long_description, + classifiers=[], + keywords="", + author=u"Brandon Liu", + author_email="brandon@protomaps.com", + url="https://github.com/protomaps/PMTiles", + license="MIT", + packages=find_packages(exclude=["ez_setup", "examples", "tests"]), + include_package_data=True, + zip_safe=False, + python_requires=">=3.7.0", + install_requires=[ + "click", + "cligj>=0.5", + "mercantile", + "pmtiles~=3.0", + "pyroaring~=1.0", + "rasterio~=1.0", + "shapely~=2.0.0", + "supermercado", + "tqdm~=4.0", + ], + extras_require={"test": ["coveralls", "pytest", "pytest-cov"]}, + entry_points=""" + [rasterio.rio_plugins] + pmtiles=rio_pmtiles.scripts.cli:pmtiles + """ + ) diff --git a/python/rio-pmtiles/tests/conftest.py b/python/rio-pmtiles/tests/conftest.py new file mode 100644 index 0000000..f6d6824 --- /dev/null +++ b/python/rio-pmtiles/tests/conftest.py @@ -0,0 +1,63 @@ +import functools +import operator +import os +import shutil +import sys + +import py +import pytest +import rasterio + +from unittest import mock + +test_files = [ + os.path.join(os.path.dirname(__file__), p) + for p in [ + "data/RGB.byte.tif", + "data/RGBA.byte.tif", + "data/rgb-193f513.vrt", + "data/rgb-fa48952.vrt", + ] +] + + +def pytest_cmdline_main(config): + # Bail if the test raster data is not present. Test data is not + # distributed with sdists since 0.12. + if functools.reduce(operator.and_, map(os.path.exists, test_files)): + print("Test data present.") + else: + print("Test data not present. See download directions in tests/README.txt") + sys.exit(1) + + +@pytest.fixture(scope="function") +def data(tmpdir): + """A temporary directory containing a copy of the files in data.""" + datadir = tmpdir.ensure("tests/data", dir=True) + for filename in test_files: + shutil.copy(filename, str(datadir)) + return datadir + + +@pytest.fixture(scope="function") +def empty_data(tmpdir): + """A temporary directory containing a folder with an empty data file.""" + filename = test_files[0] + out_filename = os.path.join(str(tmpdir), "empty.tif") + with rasterio.open(filename, "r") as src: + with rasterio.open(out_filename, "w", **src.meta) as dst: + pass + return out_filename + + +@pytest.fixture() +def rgba_cutline_path(): + """Path to a GeoJSON rhombus within the extents of RGBA.byte.tif""" + return os.path.join(os.path.dirname(__file__), "data/rgba_cutline.geojson") + + +@pytest.fixture() +def rgba_points_path(): + """Path to a pair of GeoJSON points. This is not a valid cutline.""" + return os.path.join(os.path.dirname(__file__), "data/rgba_points.geojson") diff --git a/python/rio-pmtiles/tests/data/RGB.byte.tif b/python/rio-pmtiles/tests/data/RGB.byte.tif new file mode 100644 index 0000000..1efaf4a Binary files /dev/null and b/python/rio-pmtiles/tests/data/RGB.byte.tif differ diff --git a/python/rio-pmtiles/tests/data/RGBA.byte.tif b/python/rio-pmtiles/tests/data/RGBA.byte.tif new file mode 100644 index 0000000..6aedc37 Binary files /dev/null and b/python/rio-pmtiles/tests/data/RGBA.byte.tif differ diff --git a/python/rio-pmtiles/tests/data/rgb-193f513.vrt b/python/rio-pmtiles/tests/data/rgb-193f513.vrt new file mode 100644 index 0000000..634027c --- /dev/null +++ b/python/rio-pmtiles/tests/data/rgb-193f513.vrt @@ -0,0 +1,40 @@ + + PROJCS["UTM Zone 18, Northern Hemisphere",GEOGCS["Unknown datum based upon the WGS 84 ellipsoid",DATUM["Not_specified_based_on_WGS_84_spheroid",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-75],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]]] + -1.0618653799999999e+06, 3.0003792667509481e+02, 0.0000000000000000e+00, 3.5412780899999999e+06, 0.0000000000000000e+00, -3.0004178272980499e+02 + + 0 + Red + + RGB.byte.tif + 1 + + + + 0 + + + + 0 + Green + + RGB.byte.tif + 2 + + + + 0 + + + + 0 + Blue + + RGB.byte.tif + 3 + + + + 0 + + + diff --git a/python/rio-pmtiles/tests/data/rgb-fa48952.vrt b/python/rio-pmtiles/tests/data/rgb-fa48952.vrt new file mode 100644 index 0000000..494604e --- /dev/null +++ b/python/rio-pmtiles/tests/data/rgb-fa48952.vrt @@ -0,0 +1,40 @@ + + PROJCS["UTM Zone 18, Northern Hemisphere",GEOGCS["Unknown datum based upon the WGS 84 ellipsoid",DATUM["Not_specified_based_on_WGS_84_spheroid",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-75],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]]] + 1.1256674000000001e+05, 3.0003792667509481e+02, 0.0000000000000000e+00, 3.5598010899999999e+06, 0.0000000000000000e+00, -3.0004178272980499e+02 + + 0 + Red + + RGB.byte.tif + 1 + + + + 0 + + + + 0 + Green + + RGB.byte.tif + 2 + + + + 0 + + + + 0 + Blue + + RGB.byte.tif + 3 + + + + 0 + + + diff --git a/python/rio-pmtiles/tests/data/rgba_cutline.geojson b/python/rio-pmtiles/tests/data/rgba_cutline.geojson new file mode 100644 index 0000000..3e8ed40 --- /dev/null +++ b/python/rio-pmtiles/tests/data/rgba_cutline.geojson @@ -0,0 +1 @@ +{"type":"FeatureCollection","features":[{"type":"Feature","properties":{},"geometry":{"type":"Polygon","coordinates":[[[-76.47171020507812,37.2198034112712],[-76.47891998291016,37.196834728499866],[-76.46072387695312,37.16770367048253],[-76.40974044799805,37.161000570006095],[-76.39875411987305,37.15224460472995],[-76.34931564331053,37.1284341983056],[-76.31103515625,37.139382442337094],[-76.30794525146484,37.16291580223116],[-76.32253646850586,37.197108206316365],[-76.37094497680664,37.23798199321937],[-76.44441604614258,37.22827818273987],[-76.47171020507812,37.2198034112712]]]}},{"type":"Feature","properties":{},"geometry":{"type":"Polygon","coordinates":[[[-78.585205078125,24.67946552658519],[-78.22265625,24.410889551000935],[-77.69119262695312,24.69942955501979],[-77.9754638671875,24.992281691278635],[-78.585205078125,24.67946552658519]]]}}]} \ No newline at end of file diff --git a/python/rio-pmtiles/tests/data/rgba_points.geojson b/python/rio-pmtiles/tests/data/rgba_points.geojson new file mode 100644 index 0000000..df9da62 --- /dev/null +++ b/python/rio-pmtiles/tests/data/rgba_points.geojson @@ -0,0 +1 @@ +{"type":"FeatureCollection","features":[{"type":"Feature","properties":{},"geometry":{"type":"Point","coordinates":[-77.82028198242188,24.720637830132038]}},{"type":"Feature","properties":{},"geometry":{"type":"Point","coordinates":[-77.607421875,24.165549146828848]}}]} \ No newline at end of file diff --git a/python/rio-pmtiles/tests/test_cli.py b/python/rio-pmtiles/tests/test_cli.py new file mode 100644 index 0000000..adac969 --- /dev/null +++ b/python/rio-pmtiles/tests/test_cli.py @@ -0,0 +1,290 @@ +import os +import sqlite3 +import sys +import warnings + +import click +from click.testing import CliRunner +import pytest +import rasterio +from rasterio.rio.main import main_group + +from pmtiles.reader import Reader, MmapSource, all_tiles +import rio_pmtiles.scripts.cli + +from conftest import mock + +class Output: + def __init__(self, fname): + self.file = open(fname, 'rb') + + def __enter__(self): + return Reader(MmapSource(self.file)) + + def __exit__(self, exc_type, exc_val, exc_tb): + self.file.close() + +def test_cli_help(): + runner = CliRunner() + result = runner.invoke(main_group, ["pmtiles", "--help"]) + assert result.exit_code == 0 + assert "Export a dataset to PMTiles." in result.output + + +@mock.patch("rio_pmtiles.scripts.cli.rasterio") +def test_dst_nodata_validation(rio): + """--dst-nodata requires source nodata in some form""" + rio.open.return_value.__enter__.return_value.profile.get.return_value = None + runner = CliRunner() + result = runner.invoke( + main_group, ["pmtiles", "--dst-nodata", "0", "in.tif", "out.pmtiles"] + ) + assert result.exit_code == 2 + + +@pytest.mark.parametrize("filename", ["RGB.byte.tif", "RGBA.byte.tif"]) +def test_export_metadata(tmpdir, data, filename): + inputfile = str(data.join(filename)) + outputfile = str(tmpdir.join("export.pmtiles")) + runner = CliRunner() + result = runner.invoke(main_group, ["pmtiles", inputfile, outputfile]) + assert result.exit_code == 0 + + with Output(outputfile) as p: + assert p.metadata()['name'] == filename + +def test_export_metadata_output_opt(tmpdir, data): + inputfile = str(data.join("RGB.byte.tif")) + outputfile = str(tmpdir.join("export.pmtiles")) + runner = CliRunner() + result = runner.invoke(main_group, ["pmtiles", inputfile, "-o", outputfile]) + assert result.exit_code == 0 + + with Output(outputfile) as p: + assert p.metadata()['name'] == "RGB.byte.tif" + +def test_export_tiles(tmpdir, data): + inputfile = str(data.join("RGB.byte.tif")) + outputfile = str(tmpdir.join("export.pmtiles")) + runner = CliRunner() + result = runner.invoke(main_group, ["pmtiles", inputfile, outputfile]) + assert result.exit_code == 0 + + with open(outputfile, 'rb') as f: + src = MmapSource(f) + assert len(list(all_tiles(src))) == 6 + +def test_export_zoom(tmpdir, data): + inputfile = str(data.join("RGB.byte.tif")) + outputfile = str(tmpdir.join("export.pmtiles")) + runner = CliRunner() + result = runner.invoke( + main_group, ["pmtiles", inputfile, outputfile, "--zoom-levels", "6..7"] + ) + assert result.exit_code == 0 + + with open(outputfile, 'rb') as f: + src = MmapSource(f) + assert len(list(all_tiles(src))) == 6 + +def test_export_jobs(tmpdir, data): + inputfile = str(data.join("RGB.byte.tif")) + outputfile = str(tmpdir.join("export.pmtiles")) + runner = CliRunner() + result = runner.invoke(main_group, ["pmtiles", inputfile, outputfile, "-j", "4"]) + assert result.exit_code == 0 + + with open(outputfile, 'rb') as f: + src = MmapSource(f) + assert len(list(all_tiles(src))) == 6 + + +def test_export_src_nodata(tmpdir, data): + inputfile = str(data.join("RGB.byte.tif")) + outputfile = str(tmpdir.join("export.pmtiles")) + runner = CliRunner() + result = runner.invoke( + main_group, + ["pmtiles", inputfile, outputfile, "--src-nodata", "0", "--dst-nodata", "0"], + ) + assert result.exit_code == 0 + + with open(outputfile, 'rb') as f: + src = MmapSource(f) + assert len(list(all_tiles(src))) == 6 + + +def test_export_bilinear(tmpdir, data): + inputfile = str(data.join("RGB.byte.tif")) + outputfile = str(tmpdir.join("export.pmtiles")) + runner = CliRunner() + result = runner.invoke( + main_group, ["pmtiles", inputfile, outputfile, "--resampling", "bilinear"] + ) + assert result.exit_code == 0 + + with open(outputfile, 'rb') as f: + src = MmapSource(f) + assert len(list(all_tiles(src))) == 6 + + +def test_skip_empty(tmpdir, empty_data): + """This file has the same shape as RGB.byte.tif, but no data.""" + inputfile = empty_data + outputfile = str(tmpdir.join("export.pmtiles")) + runner = CliRunner() + result = runner.invoke(main_group, ["pmtiles", inputfile, outputfile]) + assert result.exit_code == 0 + + with open(outputfile, 'rb') as f: + src = MmapSource(f) + assert len(list(all_tiles(src))) == 0 + + +def test_include_empty(tmpdir, empty_data): + """This file has the same shape as RGB.byte.tif, but no data.""" + inputfile = empty_data + outputfile = str(tmpdir.join("export.pmtiles")) + runner = CliRunner() + result = runner.invoke( + main_group, ["pmtiles", "--include-empty-tiles", inputfile, outputfile] + ) + assert result.exit_code == 0 + + with open(outputfile, 'rb') as f: + src = MmapSource(f) + assert len(list(all_tiles(src))) == 6 + + +def test_invalid_format_rgba(tmpdir, empty_data): + """--format JPEG --rgba is not allowed""" + inputfile = empty_data + outputfile = str(tmpdir.join("export.pmtiles")) + runner = CliRunner() + result = runner.invoke( + main_group, ["pmtiles", "--format", "JPEG", "--rgba", inputfile, outputfile] + ) + assert result.exit_code == 2 + + +@pytest.mark.parametrize("filename", ["RGBA.byte.tif", "RGB.byte.tif"]) +def test_rgba_png(tmpdir, data, filename): + inputfile = str(data.join(filename)) + outputfile = str(tmpdir.join("export.pmtiles")) + runner = CliRunner() + result = runner.invoke( + main_group, ["pmtiles", "--rgba", "--format", "PNG", inputfile, outputfile] + ) + with Output(outputfile) as p: + assert p.metadata()['name'] == filename + + +@pytest.mark.parametrize( + "minzoom,maxzoom,exp_num_tiles,source", + [ + (4, 10, 70, "RGB.byte.tif"), + (6, 7, 6, "RGB.byte.tif"), + (4, 10, 12, "rgb-193f513.vrt"), + (4, 10, 70, "rgb-fa48952.vrt"), + ], +) +def test_export_count(tmpdir, data, minzoom, maxzoom, exp_num_tiles, source): + inputfile = str(data.join(source)) + outputfile = str(tmpdir.join("export.pmtiles")) + runner = CliRunner() + result = runner.invoke( + main_group, + [ + "pmtiles", + "--zoom-levels", + "{}..{}".format(minzoom, maxzoom), + inputfile, + outputfile, + ], + ) + assert result.exit_code == 0 + with open(outputfile, 'rb') as f: + src = MmapSource(f) + assert len(list(all_tiles(src))) == exp_num_tiles + + +@pytest.mark.parametrize("filename", ["RGBA.byte.tif"]) +def test_progress_bar(tmpdir, data, filename): + inputfile = str(data.join(filename)) + outputfile = str(tmpdir.join("export.pmtiles")) + runner = CliRunner() + result = runner.invoke( + main_group, + [ + "pmtiles", + "-#", + "--zoom-levels", + "4..11", + "--rgba", + "--format", + "PNG", + inputfile, + outputfile, + ], + ) + assert result.exit_code == 0 + assert "100%" in result.output + + +@pytest.mark.parametrize("inputfiles", [[], ["a.tif", "b.tif"]]) +def test_input_required(inputfiles): + """We require exactly one input file""" + runner = CliRunner() + result = runner.invoke(main_group, ["pmtiles"] + inputfiles + ["foo.pmtiles"]) + assert result.exit_code == 2 + + +@pytest.mark.parametrize("filename", ["RGBA.byte.tif"]) +def test_cutline_progress_bar(tmpdir, data, rgba_cutline_path, filename): + """rio-pmtiles accepts and uses a cutline""" + inputfile = str(data.join(filename)) + outputfile = str(tmpdir.join("export.pmtiles")) + runner = CliRunner() + result = runner.invoke( + main_group, + [ + "pmtiles", + "-#", + "--zoom-levels", + "4..11", + "--rgba", + "--format", + "PNG", + "--cutline", + rgba_cutline_path, + inputfile, + outputfile, + ], + ) + assert result.exit_code == 0 + assert "100%" in result.output + + +@pytest.mark.parametrize("filename", ["RGBA.byte.tif"]) +def test_invalid_cutline(tmpdir, data, rgba_points_path, filename): + """Points cannot serve as a cutline""" + inputfile = str(data.join(filename)) + outputfile = str(tmpdir.join("export.pmtiles")) + runner = CliRunner() + result = runner.invoke( + main_group, + [ + "pmtiles", + "-#", + "--zoom-levels", + "4..11", + "--rgba", + "--format", + "PNG", + "--cutline", + rgba_points_path, + inputfile, + outputfile, + ], + ) + assert result.exit_code == 1 diff --git a/python/rio-pmtiles/tests/test_mod.py b/python/rio-pmtiles/tests/test_mod.py new file mode 100644 index 0000000..c33ba55 --- /dev/null +++ b/python/rio-pmtiles/tests/test_mod.py @@ -0,0 +1,31 @@ +"""Module tests""" + +from mercantile import Tile +import pytest + +import rio_pmtiles.worker + + +@pytest.mark.parametrize("tile", [Tile(36, 73, 7), Tile(0, 0, 0), Tile(1, 1, 1)]) +@pytest.mark.parametrize("filename", ["RGB.byte.tif", "RGBA.byte.tif"]) +def test_process_tile(data, filename, tile): + sourcepath = str(data.join(filename)) + rio_pmtiles.worker.init_worker( + sourcepath, + { + "driver": "PNG", + "dtype": "uint8", + "nodata": 0, + "height": 256, + "width": 256, + "count": 3, + "crs": "EPSG:3857", + }, + "nearest", + {}, + {}, + ) + t, contents = rio_pmtiles.worker.process_tile(tile) + assert t.x == tile.x + assert t.y == tile.y + assert t.z == tile.z