Add script to clip the global layer to an area of interest - #6
Open
Caleb Robinson (calebrob6) wants to merge 3 commits into
Open
Add script to clip the global layer to an area of interest#6Caleb Robinson (calebrob6) wants to merge 3 commits into
Caleb Robinson (calebrob6) wants to merge 3 commits into
Conversation
The global layer ships as thousands of Planet L15 quads, so using it for a country or region means locating the right quads, mosaicking them and clipping the result. This adds scripts/analysis/clip-to-boundary.py, which does that in one step and writes a single two-band COG. Only the tiles overlapping the area of interest are read, streamed over HTTP via /vsicurl, so nothing but the tile index is downloaded. The output keeps the source CRS, resolution and pixel grid and uses nearest-neighbour resampling on a grid-aligned extent, so its pixels are an exact copy of the source plus boundary masking. The area of interest can be any OGR-readable vector file or a country ISO3 code, in which case the ADM0 boundary is fetched from fieldmaps.io. Available quarters are discovered from the tile index schema rather than hard coded, so new releases work without changing the script. An optional --buffer-m keeps a margin of real data around the boundary so that focal and zonal statistics near the edge are not computed against NoData. The margin is a ground distance: EPSG:3857 is not conformal with respect to the WGS84 ellipsoid, so the distance is scaled by the local meridian scale factor and the result is verified geodesically. Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
Address review feedback: - Replace the ground-distance buffer with a plain --buffer in projected units. The previous version scaled the distance by the local Mercator meridian scale factor and verified it geodesically, which was a lot of machinery for an option that is off by default. Buffering is now a single call and no buffer is applied unless asked for. - Use rasterio to read the mosaic's grid, set band names and metadata, and write the COG, instead of the osgeo bindings. This drops the osgeo import entirely, and with it the gdal.UseExceptions() call that was only needed to make those bindings raise on failure. gdalbuildvrt and gdalwarp are still invoked as command line tools so the mosaic is never held in memory. - Trim explanatory comments down to what is needed to follow the code. Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
Findings from a review pass, in rough order of severity: - gdalbuildvrt skips tiles it cannot open and still exits 0, so a transient HTTP failure would have produced a COG with a NoData hole in it and no indication anything was wrong. Build the mosaic with -strict instead. - The tile index was downloaded by truncating the destination and then fetching, so an interrupted download left an empty file that later runs happily accepted. Stream it to a temporary file and move it into place once it is complete. This also avoids holding the whole index in memory. - Extent snapping used bare floor/ceil, so a coordinate landing on a grid line could come out as 2560.0000000000177 and gain a spurious extra row or column. Round to the nearest whole pixel when within a small tolerance. - Check the mosaic really is north-up, unrotated, square-pixelled and in the expected CRS before relying on the output being pixel-exact. Future quarters are discovered automatically, so these assumptions are worth asserting. - Warn when quads have no tile for the requested quarter, and when the area of interest reaches past the available tiles, rather than quietly returning a smaller raster than asked for. - Reject areas of interest spanning more than 180 degrees of longitude, which cannot be represented as one contiguous EPSG:3857 extent. - Reject non-finite --buffer values, which previously passed validation and then behaved like no buffer at all. - Write the output to a temporary path and move it into place, so a failure part way through does not leave a half written raster. - Add line continuations to the README examples so they can be pasted into a shell as-is. Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
What
Adds
scripts/analysis/clip-to-boundary.py, which turns the global tile index into a singletwo-band COG for an area of interest, plus a README section documenting it.
The global layer ships as thousands of Planet L15 quads, so using it for a country or region
currently means locating the right quads, mosaicking them and clipping the result. This does
that in one step:
Design notes
HTTP with
/vsicurl/. Just the tile index is fetched to disk (into the gitignoreddata/).extent is snapped outward onto the source grid, and resampling is nearest neighbour — so the
result is a copy of the source pixels plus boundary masking, with no half-pixel shift.
data_*columns) rather than hard coded,so when new quarters are published the script picks them up unchanged.
--list-quartersprints what's currently available.
--boundary, optional--layer) or acountry code (
--iso3), which pulls the geoBoundaries ADM0 from fieldmaps.io.--bufferoptionally keeps a margin of data around the AOI, in EPSG:3857 units, so thatfocal or zonal statistics near the edge aren't computed against NoData. Off by default.
COG.
gdalbuildvrtandgdalwarpare invoked as command line tools so the mosaic is neverheld in memory. Note that metadata has to be set on the clipped VRT rather than on the
finished file, because the COG driver is copy-only and reopening the output would break its
layout.
Follows the conventions in
scripts/: MIT header,argparse, module docstring with usageexamples,
loguru. No new dependencies beyond whatrequirements-scripts.txtalready lists,plus the GDAL command line tools.
Testing
a COG (
validate_cloud_optimized_geotiff.py), with two correctly named bands, NoData -1,overviews and metadata tags present.
same footprint in the output: identical values, max absolute difference 0.0 over 131,070
pixels inside the cutline.
--buffer 1000— tile count grows 139 -> 142 and the grid by ~26 px (1000 / 76.437), asexpected. Negative values are rejected.
--layerselection, mutual exclusion of--boundary/--iso3, unknownquarter rejection, a single-tile AOI, and an AOI with no coverage (clean error, no traceback).
Notes for reviewers
--iso3convenience is the only third-party network dependency (fieldmaps.io) and isisolated to one function. Happy to drop it and keep
--boundaryonly if you'd prefer therepo not reach out to a third-party host.
published tiles are 256x256 per quad at ~76.4 m/px. Left alone here — happy to fix it in a
separate PR if that wording should change.