Cloud-Native Geospatial Formats

Cloud-native formats are designed to be read in place over HTTP, a few bytes at a time, instead of downloaded whole. Cloud Optimized GeoTIFF (COG) does this for rasters; GeoParquet and FlatGeobuf do it for vectors; PMTiles does it for tiles. Together they let a pipeline query terabytes on object storage while pulling only the windows it needs. This is the storage-and-delivery stage of Geospatial Data Ingestion & Processing Workflows — the format you commit to after files land through Shapefile & GeoJSON Parsing and are normalized by Coordinate Reference System Transformations. It also feeds the engines in DuckDB Spatial Analytics and the delivery layer in Web Mapping & Interactive Visualization.

Cloud-native read model An object store holds COG, GeoParquet, and PMTiles files with internal indexes; clients issue HTTP range requests to read only the windows, row groups, or tiles they need. Read the bytes you need, not the whole file Object store (S3/R2) COG — raster windows GeoParquet — row groups PMTiles — map tiles Python client rasterio / duckdb reads a window HTTP range request internal index → byte offset
Each cloud-native format carries an internal index so a client can range-request exactly the window, row group, or tile it needs.

Architecture & Data Structures

Every cloud-native format solves the same problem in the same way: an internal index describes where each chunk of data physically lives, and clients fetch chunks with HTTP range requests (Range: bytes=8192-16383) rather than a full GET. The file lives on plain object storage — S3, Cloudflare R2, Google Cloud Storage, Azure Blob — with no database process and no tile server in front of it. All the intelligence about which bytes to ask for lives in the client library.

The four formats differ only in how they lay out data and index it:

Two more layouts round out the family. Zarr (and its geospatial profile, GeoZarr) breaks an n-dimensional array into a directory of independently-addressable chunks with a small JSON metadata document at the root — the natural home for time-series and multi-variable data cubes where COG's two-dimensional tiling runs out of dimensions, and the format behind most of the workflows in Xarray & rioxarray Raster Cubes. And STAC is not a data format at all but the catalogue layer above them: a JSON description of which assets exist, where, and over what footprint and time range, so a client can decide which COG to open before it opens anything. In a mature cloud-native stack, STAC answers "which files", the internal index answers "which bytes".

The invariant behind all of them is worth stating plainly, because it is what breaks when a deployment misbehaves: the format guarantees that a small, bounded read from a known location yields enough information to compute the byte range of everything else. That contract needs exactly two things from the storage layer — support for HTTP Range requests, and a 206 Partial Content response rather than a helpful 200 OK with the whole object. Nothing else about the server matters. A CDN that buffers and re-serves complete objects turns every cloud-native format back into a download, and no amount of client tuning recovers it.

A minimal COG open demonstrates the model — nothing downloads until you ask for pixels, and even then only the overlapping tiles transfer:

import rasterio
from rasterio.windows import Window

# Read a single 512x512 window from a COG on S3 — only that window transfers
cog_url = "https://example-bucket.s3.amazonaws.com/sentinel_ortho_cog.tif"
with rasterio.open(cog_url) as src:
    print(src.profile["driver"], src.block_shapes[0])  # GTiff (512, 512) — internally tiled
    window = Window(col_off=4096, row_off=4096, width=512, height=512)
    patch = src.read(1, window=window)
    print(patch.shape)   # (512, 512) — the full image was never localized

The reciprocal for vectors is a row-group-aware read: the reader inspects Parquet footer statistics, discards groups whose bounding box misses the query, and decodes only the survivors.

Cloud-native format comparison matrix A four-by-four matrix comparing COG, GeoParquet, FlatGeobuf, and PMTiles across data unit, index structure, primary use, and typical reader. COG stores raster tiles indexed by a tile offset table, read with rasterio. GeoParquet stores row groups indexed by row-group bounding-box statistics, read with duckdb and pyarrow. FlatGeobuf stores features indexed by a packed R-tree, read with fiona. PMTiles stores map tiles indexed by a z/x/y directory, read with maplibre. Same model, four layouts: what each format indexes Format → COG GeoParquet FlatGeobuf PMTiles Data unit Index structure Primary use Typical reader tile row group feature map tile tile offset table row-group bbox stats packed R-tree z/x/y directory raster analytical vector streaming vector rendered tiles rasterio duckdb + pyarrow fiona maplibre
All four formats share the range-request model; they differ only in the unit they store, the index that maps it to byte offsets, and the reader that speaks their layout.

Environment Configuration & Dependency Resolution

Cloud-native reads lean on GDAL's virtual filesystem layer (/vsis3/, /vsicurl/, /vsigs/, /vsiaz/), so the binding versions and their bundled GDAL matter more than usual. Install from conda-forge to keep GDAL, PROJ, and the Python bindings ABI-compatible:

conda install -c conda-forge \
    "rasterio=1.3.*" "geopandas=0.14.*" "pyogrio=0.7.*" \
    "pyarrow=15.*" "gdal=3.8.*" "duckdb=0.10.*"
pip install "pmtiles>=3.2"          # PMTiles reader/writer, no GDAL dependency

Three things commonly bite here. GDAL below 3.1 cannot write COGs with the dedicated COG driver; GeoParquet round trips need pyarrow>=7 on both the writer and reader, and version skew between them is the single most common read failure; and pyogrio (not the older fiona path) is what pushes bounding-box filters down into GeoParquet and FlatGeobuf reads. The same rasterio install covered in Raster Data Handling with Rasterio provides COG read and write support out of the box.

For authenticated cloud reads, set the GDAL environment knobs so it issues efficient range requests instead of listing buckets or re-reading directories:

import rasterio

gdal_env = {
    "AWS_S3_ENDPOINT": "s3.amazonaws.com",
    "GDAL_DISABLE_READDIR_ON_OPEN": "EMPTY_DIR",     # never LIST the bucket on open
    "CPL_VSIL_CURL_ALLOWED_EXTENSIONS": ".tif",      # skip sidecar probes (.ovr, .aux.xml)
    "GDAL_HTTP_MULTIPLEX": "YES",                     # HTTP/2 multiplexed range requests
    "VSI_CACHE": "TRUE",                             # cache fetched blocks in memory
}
with rasterio.Env(**gdal_env):
    with rasterio.open("/vsis3/example-bucket/sentinel_ortho_cog.tif") as src:
        overview_factors = src.overviews(1)          # e.g. [2, 4, 8, 16] — decimation levels

GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR is the highest-impact setting: without it, GDAL LISTs the whole prefix on every open, which on a bucket with millions of objects turns a sub-second read into minutes. Credentials come from the standard AWS_ACCESS_KEY_ID/AWS_SECRET_ACCESS_KEY (or an attached IAM role); for public buckets set AWS_NO_SIGN_REQUEST=YES, and for buckets that bill the reader set AWS_REQUEST_PAYER=requester — omitting it on a requester-pays bucket produces a bare 403 that looks identical to a credentials problem.

Which capabilities you actually have is a property of the GDAL build, not just its version number, and this is where environments diverge most:

Capability Needs How to check
Write a real COG in one pass GDAL 3.1+ (COG driver) gdalinfo --formats | grep COG
Read/write Parquet and Arrow GDAL 3.5+ built with Arrow ogrinfo --formats | grep -i parquet
LERC raster compression GDAL built with liblerc gdalinfo --format GTiff | grep LERC
/vsiaz/ Azure blob access GDAL 3.0+ gdalinfo --formats and a test open
Spatial pruning in DuckDB spatial + httpfs extensions SELECT * FROM duckdb_extensions()

The Parquet row is the one that surprises people: a Linux distribution package of GDAL 3.8 may have no Parquet driver at all, while a conda-forge build of 3.6 does. The version number tells you nothing on its own, which is why the check belongs in your environment smoke test rather than in a README.

Two more knobs are worth knowing before you need them. VSI_CACHE_SIZE sets the per-file-handle byte cache (default 25 MB) and is what makes repeated reads of the same header free. CPL_VSIL_CURL_USE_HEAD=NO suppresses the HEAD request GDAL issues to learn an object's size, which is worth setting when your storage layer answers HEAD slowly or not at all — some signed-URL schemes sign the GET only.

HTTP calls GDAL makes when opening a remote COG, before and after tuning Two side-by-side call sequences for a single open of a Cloud Optimized GeoTIFF on S3. With default settings GDAL makes five round trips: a LIST of the whole bucket prefix, two HEAD probes for sidecar files that return 404, then a GET of the header and a GET of one tile. With directory listing disabled and a curl extension allowlist set, the first three calls are skipped and only the two GETs that carry data remain. One open() call: what GDAL fetches before the first pixel Default settings 5 round trips, 3 of them wasted Listing disabled + extension allowlist 2 round trips, both carry data LIST prefix — every key in the bucket HEAD ortho_cog.tif.ovr → 404 HEAD ortho_cog.tif.aux.xml → 404 GET bytes 0-16383 — header + IFDs GET one 512-pixel tile byte range skipped — READDIR_ON_OPEN is EMPTY_DIR skipped — extension allowlist is .tif skipped — extension allowlist is .tif GET bytes 0-16383 — header + IFDs GET one 512-pixel tile byte range On a bucket holding millions of objects the LIST alone turns a sub-second window read into minutes.
The two environment knobs remove three latency-only round trips from every remote open, leaving just the header read and the tile read.

Vectorized Operations & Core Workflow

The everyday shape of a cloud-native pipeline: keep the canonical analytical dataset as GeoParquet, query bounding-box windows with GeoPandas or DuckDB in place, and derive COGs and PMTiles as rendering artifacts when you need them. The windowed-raster recipe is expanded in Windowed Reads from Cloud Optimized GeoTIFF.

Writing GeoParquet is a one-liner that preserves the CRS and produces row-group statistics automatically:

import geopandas as gpd

# Normalize to a documented CRS, then write GeoParquet
parcels = gpd.read_file("parcels.gpkg").to_crs(epsg=4326)
parcels.to_parquet(
    "parcels.parquet",
    compression="zstd",       # smaller than snappy at similar decode speed
    geometry_encoding="WKB",  # default; GeoArrow available in geopandas>=1.0
    row_group_size=50_000,    # tune so a group is ~64-128 MB for skip efficiency
)

# Read back only a bounding-box window — pyogrio/pyarrow prunes row groups by bbox
aoi = (7.60, 45.00, 7.80, 45.10)   # xmin, ymin, xmax, ymax in EPSG:4326
window = gpd.read_parquet("parcels.parquet", bbox=aoi)
print(len(window), "features touched the window")

The same window read works against a remote URL — gpd.read_parquet("s3://bucket/parcels.parquet", bbox=aoi) — with s3fs installed, and only the intersecting row groups plus the footer transfer. For SQL-style analytics over the same file, DuckDB reads Parquet natively and pushes spatial predicates down; the join and aggregation patterns live in DuckDB Spatial Analytics:

import duckdb

con = duckdb.connect()
con.execute("INSTALL spatial; LOAD spatial; INSTALL httpfs; LOAD httpfs;")

# Query a remote GeoParquet without ever downloading the whole file
rows = con.execute("""
    SELECT parcel_id, ST_Area(ST_Transform(geometry, 'EPSG:4326', 'EPSG:32632')) AS area_m2
    FROM 's3://example-bucket/parcels.parquet'
    WHERE ST_Intersects(geometry, ST_MakeEnvelope(7.60, 45.00, 7.80, 45.10))
    ORDER BY area_m2 DESC
    LIMIT 10
""").fetchall()

Once a store grows past a few tens of gigabytes, stop writing one file and write a partitioned dataset instead: a directory tree of Parquet parts under Hive-style key=value subdirectories. Partition pruning happens from the directory names alone, before any footer is parsed, so it is strictly cheaper than row-group skipping and composes with it. Partition on whatever your queries filter by most — region, acquisition date, administrative level — and keep each part in the hundreds of megabytes.

import geopandas as gpd

# Read one partition out of s3://example-bucket/parcels/region=*/  without listing the rest
piemonte = gpd.read_parquet(
    "s3://example-bucket/parcels/",
    filters=[("region", "==", "piemonte")],   # partition pruning, evaluated on paths
    columns=["parcel_id", "land_use", "geometry"],   # column pruning, evaluated on the footer
)
print(piemonte.crs, len(piemonte))

The vector counterpart for feature-at-a-time consumers is FlatGeobuf, whose packed R-tree answers a bounding-box query without a columnar engine at all. pyogrio pushes the box down into GDAL, so only the intersecting features are decoded:

import pyogrio

# The R-tree is read first; only features whose entry overlaps the box are fetched
sensors = pyogrio.read_dataframe(
    "https://example-bucket.s3.amazonaws.com/sensors.fgb",
    bbox=(7.60, 45.00, 7.80, 45.10),   # in the file's own CRS
)

Note the comment: FlatGeobuf's bbox argument, like a raster Window, is interpreted in the file's CRS, not in whatever CRS you happen to be thinking in. That asymmetry with GeoParquet — where the covering columns are stored in the geometry column's declared CRS — is a recurring source of empty result sets.

Geometry / Data Processing Details

GeoParquet keeps geometry as WKB with CRS metadata (a PROJJSON block) in the file's schema, so a round trip is lossless — unlike Shapefile, which truncates field names to 10 characters, splits large layers at 2 GB, and demotes the CRS to a sidecar .prj that is easily lost. That fidelity is why GeoParquet is the right interchange format between processing stages rather than at the edges only; the full trade-off against legacy formats is worked through in GeoParquet vs Shapefile for Storage.

Two processing details determine whether cloud-native storage actually performs:

Spatial ordering drives skip efficiency. Row-group and R-tree pruning only help when nearby features sit near each other in the file. If features are written in random or insertion order, every row group's bounding box spans the whole dataset and no group can be skipped. Sort by a space-filling curve before writing:

import geopandas as gpd
from shapely import STRtree

buildings = gpd.read_parquet("buildings.parquet")
# Hilbert-curve ordering groups spatially-adjacent rows into the same row group
buildings = buildings.sort_values(
    by="geometry",
    key=lambda geom: geom.hilbert_distance(total_bounds=buildings.total_bounds),
)
buildings.to_parquet("buildings_sorted.parquet", row_group_size=50_000)

Overviews decide COG read cost at low zoom. A COG without overviews forces a full-resolution read even when the client only needs a thumbnail, so build a decimated pyramid at encode time. The rio cogeo tool does both the tiling and the overviews in one pass:

# Re-encode a plain GeoTIFF into a genuine COG: internal tiles + overviews
rio cogeo create dem_raw.tif dem_cog.tif \
    --cog-profile deflate --overview-resampling average --blocksize 512
rio cogeo validate dem_cog.tif   # confirms tiling + overviews, not just a renamed .tif

Compression is a per-format decision, not a global one. For GeoParquet, ZSTD at its default level is close to strictly better than Snappy: noticeably smaller files at comparable decode speed, and the decode is rarely the bottleneck when bytes travel over a network. For rasters, DEFLATE is the safe interoperable choice, LERC (with a controlled error bound) is dramatically smaller for continuous data like elevation, and JPEG is appropriate only for visual RGB imagery where lossy artefacts do not propagate into analysis. The failure mode is picking a raster codec the consumer's GDAL was not built with — a file that opens fine on your machine and raises a cryptic band-read error on theirs.

PMTiles is readable from Python, not only from a browser. That matters for validating a published archive without loading a web map:

from pmtiles.reader import Reader, MmapSource

with open("buildings.pmtiles", "rb") as f:
    reader = Reader(MmapSource(f))
    header = reader.header()
    print(header["min_zoom"], header["max_zoom"])   # 0 14
    tile = reader.get(12, 2200, 1345)               # one z/x/y tile as raw bytes
    print(len(tile) if tile else "no tile at that address")

A None here is diagnostic rather than an error: it means the archive genuinely has no tile at that address, usually because the source data did not cover it or the zoom range is narrower than the map requests.

For very large vector sources — global building footprints, road networks, administrative boundaries — never materialize the whole thing. Stream it feature-by-feature or partition-by-partition; the pattern for pulling continental extracts without exhausting memory is in Streaming Overture Maps Data with DuckDB.

CRS Alignment & Projection Pipeline

Cloud-native vector formats embed CRS metadata, so the discipline is to tag data correctly at write time and reproject deliberately at read time — never to leave a file's CRS ambiguous. GeoParquet stores a PROJJSON CRS in the geo metadata; COGs store theirs in GeoTIFF geokeys; both survive a round trip. Align projections with Coordinate Systems with PyProj, and keep one rule front of mind: reproject the windowed result, not the whole dataset, so you never pay to transform data you did not read.

import geopandas as gpd

fields = gpd.read_parquet("fields.parquet")
print(fields.crs)                       # EPSG:4326 — read straight from file metadata

# Read a small window in the file's native CRS, then reproject ONLY that subset
subset = gpd.read_parquet("fields.parquet", bbox=(10.0, 45.0, 10.2, 45.2))
subset_utm = subset.to_crs(epsg=32632)  # UTM 32N — a metric CRS, not Web Mercator
subset_utm["area_ha"] = subset_utm.geometry.area / 1e4   # areas are meaningful now

Two CRS traps specific to cloud reads. First, when you convert a geographic bounding box into a raster Window, the bounds must be in the raster's CRS — reproject them first with a pyproj.Transformer using always_xy=True, or the window lands in the wrong place. Second, do not use Web Mercator (EPSG:3857) for any metric result: its scale distortion grows with latitude, so an area or distance computed there is wrong by tens of percent away from the equator. Reproject windowed results into a local projected CRS (a UTM zone or a national grid) before measuring. PMTiles is the one place Web Mercator is correct — it is a rendering format whose tile grid is defined in EPSG:3857, and it is a display artifact, not an analytical one.

Production Export & Integration

Pick the output format by role, not by habit:

A short pre-flight checklist before you publish:

The cost model is requests plus egress, not storage

Object storage is cheap; the two lines that grow are per-request charges and per-gigabyte egress, and cloud-native formats move you decisively from the second to the first. A windowed COG read is a handful of GETs carrying megabytes where a download was one GET carrying gigabytes — an enormous egress win and a modest request-count cost. The pathological shape is the opposite: thousands of tiny scattered reads, each fetching far less than the minimum useful chunk, where you pay per request for bytes you barely use. That is exactly what an unaligned window loop or a 16 KB curl chunk size produces.

Three decisions follow from that. Prefer fewer, larger, block-aligned reads over many small ones. Put a CDN in front of anything read repeatedly, with long Cache-Control: max-age and a version in the filename so a rebuild never mixes old and new bytes in a client's cache — the deployment pattern described in deploying a Python map app behind a CDN. And check whether your storage provider charges egress at all: zero-egress object stores change the arithmetic enough that a format decision made on an egress-billed provider may not be the right one elsewhere.

Windows / Platform Edge Cases & Debugging

Frequently Asked Questions

Do I need all four formats, or can I standardise on one? Standardise on two: GeoParquet for analytical vector data and COG for rasters. Those cover the overwhelming majority of pipelines and both are readable by everything in the Python stack. Add PMTiles only when you are actually serving a web map, and FlatGeobuf only when a consumer needs streaming feature-at-a-time access without a columnar engine. Adding a format costs you a conversion step, a validation step, and a new failure mode, so each one should earn its place.

Is cloud-native worth it if my data lives on a local disk? Partly. The internal indexing still pays off — column pruning, row-group skipping, and overview reads are just as valuable against a local SSD, and GeoParquet is dramatically faster to read than a Shapefile of the same content regardless of where it sits. What you do not get is the range-request story, and on a fast local disk the difference between a windowed read and a full read narrows considerably. The strongest local argument is fidelity and size rather than speed.

How do I know whether a file someone sent me is genuinely cloud-optimized? Validate it, do not trust the extension. rio cogeo validate reports whether a GeoTIFF is internally tiled with overviews or merely renamed. For GeoParquet, read the footer and check both the row-group count and whether the geo metadata declares a covering column — a single row group means no skipping is possible, and a missing covering declaration means a reader has nothing to prune with. Both checks belong in the ingestion step, not in a code review.

Should the storage CRS be Web Mercator so web maps are fast? No, except for tiles. Store analytical data in a CRS suited to measurement and let the rendering layer reproject; tile pyramids are the one artefact whose grid is defined in EPSG:3857, and they are derived outputs. Storing your canonical data in Web Mercator to save a reprojection means every area, length, and buffer computed from it is distorted, and the distortion grows with latitude. Reproject at the edge, not at the source — the projection discipline covered in Coordinate Systems with PyProj.

When does the whole cloud-native model stop being the right answer? When the workload is transactional rather than analytical. These formats are immutable: appending or editing means rewriting a file or a partition, which is fine for datasets refreshed on a schedule and hopeless for data edited row-by-row by concurrent users. At that point you want a database with real transactions and indexes — the trade-off examined in PostGIS Integration with Python. The two coexist happily: transact in the database, publish snapshots as cloud-native files.

Can I convert between these formats without losing anything? Between GeoParquet and FlatGeobuf, essentially yes — both keep full field names, real types, and an embedded CRS. Converting to tiles or to a raster is lossy by definition: PMTiles generation simplifies geometry per zoom level and drops attributes you did not ask it to keep, and rasterising vectors discards topology entirely. Treat those as one-way derivations from a canonical GeoParquet or COG source, regenerated rather than round-tripped, which is how the pipeline in Generating PMTiles from GeoParquet is structured.

Related