How to Write Partitioned GeoParquet and Query It with Filters

Problem statement

GeoParquet is dramatically smaller and faster than the alternatives. Measured on 63,517 real OpenStreetMap building polygons in central Manchester (Example 4):

GeoParquet               5.69 MB     read 0.05 s
GeoPackage              17.33 MB     read 0.11 s
GeoJSON                 24.25 MB     read 0.44 s
Shapefile (all parts)   32.60 MB     read 0.12 s

Almost six times smaller than a shapefile and about nine times faster to read than GeoJSON. But two features people expect from it need explicit setup, and without that setup they silently do nothing:

• Spatial filtering on read needs a covering-bbox column, or geopandas refuses.
• Row-group skipping needs more than one row group, and the default is often one.

Quick answer

import geopandas as gpd

# spatial filtering needs a bbox column written in
gdf.to_parquet("buildings.parquet", write_covering_bbox=True)

subset = gpd.read_parquet("buildings.parquet",
                          bbox=(-2.25, 53.47, -2.23, 53.49))

Without it:

ValueError: Specifying 'bbox' not supported for this Parquet file
            (it should either have a bbox covering column or use
             'point' encoding).

The covering-bbox column cost 42% more file size here — 8.09 MB against 5.69 MB — and that is the price of the feature.

1 row groups, 63,517 rows

Row groups, a bbox column, spatial sorting, partitions

1. Check how many row groups you actually have

import pyarrow.parquet as pq

f = pq.ParquetFile("buildings.parquet")
print(f"{f.metadata.num_row_groups} row groups, {f.metadata.num_rows:,} rows")

1 row groups, 63,517 rows

One row group. Predicate pushdown works by reading each row group's column statistics and skipping groups that cannot match — with one group there is nothing to skip, and every filter reads the whole file.

The default row-group size is around a million rows, so most single-file GeoParquet has exactly one.

gdf.to_parquet(path, row_group_size=50_000)

2. Add the covering-bbox column if you want spatial filtering

gdf.to_parquet(path, write_covering_bbox=True)

This adds four float columns — xmin, ymin, xmax, ymax — per row. Their per-row-group statistics let a reader skip groups whose bounding boxes miss the query.

The cost is real: 8.09 MB against 5.69 MB, a factor of 1.42. Worth it for a file that will be queried spatially many times; not worth it for one read end to end.

3. Sort spatially before writing

This is the step that makes the statistics useful. Row-group statistics only help if the rows in a group are spatially close — otherwise every group's bounding box covers the whole dataset and none can be skipped.

gdf = gdf.sort_values("hilbert")     # or geohash, or Morton order
gdf.to_parquet(path, row_group_size=50_000, write_covering_bbox=True)

Without a spatial sort, adding the bbox column and shrinking the row groups buys almost nothing.

4. Partition by a directory key for very large datasets

buildings/
  tile=530_-22/part.parquet
  tile=530_-23/part.parquet

Hive-style partitioning puts the key in the path, so a reader can skip whole files without opening them. Measured on the 63,517 buildings in 24 tiles of 0.025° (Example 4):

read the largest tile ( 8,890 rows)   0.021 s
read every tile      (63,517 rows)   0.444 s

Roughly proportional to rows, which is what you want — the partitioning removed the rest entirely.

Note that geopandas.to_parquet does not accept partition_cols; write the partitions yourself with a groupby, or use dask_geopandas.

5. Choose the partition key from the query

Partition by whatever you filter on most: a tile, a region, a date, a category. Partitioning by something you never filter on adds directories and no benefit.

Filterable and partitioned Manchester buildings

The data used on this page

Every building polygon in OpenStreetMap inside −2.31 to −2.18° E, 53.44 to 53.52° N — central Manchester and Salford — fetched with OSMnx, keeping the building type, name, building:levels and addr:street tags, as gdf. It is the same extract the web-map guides use; OpenStreetMap is edited daily, so the setup reads the answer both pages share, saved on 9 October 2026.

import io
import os
import requests
import time
import zipfile

import geopandas as gpd
import osmnx as ox

# OpenStreetMap as Overpass returned it on 9 October 2026, saved as osmnx's own cache so the counts
# below do not drift as OSM is edited; delete the cache folder to query live instead
if not os.path.exists("cache"):
    reply = requests.get("https://data.spatialworkflow.io/osm-cache-manchester-buildings-2026-10-09.zip", timeout=300)
    reply.raise_for_status()
    with zipfile.ZipFile(io.BytesIO(reply.content)) as archive:
        archive.extractall("cache")

ox.settings.overpass_url = "https://maps.mail.ru/osm/tools/overpass/api"
ox.settings.requests_timeout = 600
if not os.path.exists("osm_buildings.parquet"):
    raw = ox.features_from_bbox((-2.31, 53.44, -2.18, 53.52), tags={"building": True}).reset_index()
    raw = raw[raw.geom_type.isin(["Polygon", "MultiPolygon"])]
    for tag in ("name", "building:levels", "addr:street"):
        if tag not in raw:
            raw[tag] = None
    raw[["id", "building", "name", "building:levels", "addr:street", "geometry"]].rename(
        columns={"building:levels": "levels", "addr:street": "street"}).to_parquet("osm_buildings.parquet")
gdf = gpd.read_parquet("osm_buildings.parquet")
print(f"{len(gdf):,} buildings")

63,517 buildings

Example 1 — writing a file that can actually be filtered

import numpy as np
import geopandas as gpd
import pyarrow.parquet as pq


def write_queryable(gdf, path, row_group_size=50_000, sort=True):
    """Spatially sorted, sensible row groups, covering bbox."""
    if sort:
        bounds = gdf.geometry.bounds
        cx = ((bounds["minx"] + bounds["maxx"]) / 2).to_numpy()
        cy = ((bounds["miny"] + bounds["maxy"]) / 2).to_numpy()

        # Morton (Z-order) interleave of scaled integer coordinates
        def scale(v):
            lo, hi = v.min(), v.max()
            return ((v - lo) / max(hi - lo, 1e-12) * (2 ** 20 - 1)).astype(np.uint64)

        def interleave(a, b):
            out = np.zeros_like(a)
            for bit in range(20):
                out |= ((a >> bit) & 1) << (2 * bit + 1)
                out |= ((b >> bit) & 1) << (2 * bit)
            return out

        gdf = gdf.iloc[np.argsort(interleave(scale(cx), scale(cy)))]

    gdf.to_parquet(path, row_group_size=row_group_size,
                   write_covering_bbox=True)

    meta = pq.ParquetFile(path).metadata
    print(f"  {meta.num_rows:,} rows in {meta.num_row_groups} row groups")
    print(f"  {meta.serialized_size / 1e3:.0f} kB of metadata")
    if meta.num_row_groups == 1:
        print("  ! one row group — no statistics-based skipping is possible")
    return path


write_queryable(gdf, "queryable.parquet", row_group_size=5_000)

  63,517 rows in 13 row groups
  37 kB of metadata

The Morton sort is the part that makes everything else work. Twenty lines, once, and every spatial filter afterwards can skip most of the file.

Example 2 — writing Hive partitions

import os
import shutil
import numpy as np
import geopandas as gpd


def write_partitioned(gdf, out_dir, cell_deg=0.025, row_group_size=50_000):
    """One file per tile, in a directory layout a reader can prune."""
    shutil.rmtree(out_dir, ignore_errors=True)
    os.makedirs(out_dir)

    centroids = gdf.geometry.representative_point()
    tile = (np.floor(centroids.x / cell_deg).astype(int).astype(str) + "_" +
            np.floor(centroids.y / cell_deg).astype(int).astype(str))

    sizes = []
    for key, part in gdf.groupby(tile.values):
        directory = os.path.join(out_dir, f"tile={key}")
        os.makedirs(directory, exist_ok=True)
        path = os.path.join(directory, "part.parquet")
        part.to_parquet(path, row_group_size=row_group_size,
                        write_covering_bbox=True)
        sizes.append((key, len(part), os.path.getsize(path)))

    total = sum(s for _, _, s in sizes)
    print(f"  {len(sizes)} partitions, {total / 1e6:.1f} MB total")
    print(f"  rows per partition: min {min(n for _, n, _ in sizes):,}, "
          f"median {int(np.median([n for _, n, _ in sizes])):,}, "
          f"max {max(n for _, n, _ in sizes):,}")
    if max(n for _, n, _ in sizes) > 20 * max(1, min(n for _, n, _ in sizes)):
        print("  ! very uneven partitions — consider a quadtree rather than "
              "a fixed grid")
    return out_dir


write_partitioned(gdf, "buildings_parts", cell_deg=0.025, row_group_size=5_000)

  24 partitions, 8.7 MB total
  rows per partition: min 264, median 1,758, max 8,890
  ! very uneven partitions — consider a quadtree rather than a fixed grid

Even partition sizes matter more than the exact scheme. A fixed grid over clustered data produces one enormous partition and many tiny ones, which is the worst of both — no pruning benefit and per-file overhead.

Example 3 — reading with the filters that exist

import geopandas as gpd
import pyarrow.dataset as ds


def read_filtered(path, bbox=None, columns=None, where=None):
    """Spatial and attribute filtering, with a clear fallback."""
    extra = {"filters": where} if where else {}        # filters=None is rejected
    try:
        gdf = gpd.read_parquet(path, bbox=bbox, columns=columns, **extra)
        print(f"  {len(gdf):,} rows via pushdown")
        return gdf
    except ValueError as exc:
        if "bbox" not in str(exc):
            raise
        print(f"  ! {exc}")
        print("  falling back to reading everything and clipping")

    gdf = gpd.read_parquet(path, columns=columns, **extra)
    if bbox:
        gdf = gdf.cx[bbox[0]:bbox[2], bbox[1]:bbox[3]]
    print(f"  {len(gdf):,} rows after clipping")
    return gdf


box = (-2.25, 53.47, -2.23, 53.49)
gdf.to_parquet("plain.parquet")
print("queryable file:")
a = read_filtered("queryable.parquet", bbox=box)
print("plain file:")
b = read_filtered("plain.parquet", bbox=box)

queryable file:
  1,856 rows via pushdown
plain file:
  ! Specifying 'bbox' not supported for this Parquet file (it should either have a bbox covering column or use 'point' encoding).
  falling back to reading everything and clipping
  1,856 rows after clipping

Making the fallback loud is the point. A silent fallback turns "spatial pushdown" into "read the whole file and filter in memory", which works and does not scale, and nobody notices until the dataset grows.

Example 4 — the figures at the top of this page

import glob


def timed_read(fn, repeat=3):
    best = float("inf")
    for _ in range(repeat):
        started = time.perf_counter()
        fn()
        best = min(best, time.perf_counter() - started)
    return best


gdf.to_parquet("f.parquet")
gdf.to_file("f.gpkg")
gdf.to_file("f.geojson", driver="GeoJSON")
gdf.to_file("f.shp")
formats = {"GeoParquet": (["f.parquet"], lambda: gpd.read_parquet("f.parquet")),
           "GeoPackage": (["f.gpkg"], lambda: gpd.read_file("f.gpkg")),
           "GeoJSON": (["f.geojson"], lambda: gpd.read_file("f.geojson")),
           "Shapefile (all parts)": (glob.glob("f.shp") + glob.glob("f.[sdcp][hbhrp][xfpjg]") + glob.glob("f.cpg"),
                                     lambda: gpd.read_file("f.shp"))}
print(f"{len(gdf):,} buildings")
for name, (files, reader) in formats.items():
    size = sum(os.path.getsize(f) for f in set(files))
    print(f"{name:22} {size / 1e6:6.2f} MB     read {timed_read(reader):.2f} s")

gdf.to_parquet("covered.parquet", write_covering_bbox=True)
print(f"\ncovering bbox: {os.path.getsize('covered.parquet') / 1e6:.2f} MB against "
      f"{os.path.getsize('f.parquet') / 1e6:.2f} MB (x{os.path.getsize('covered.parquet') / os.path.getsize('f.parquet'):.2f})")
import shapely

print(f"local read: whole plain file {timed_read(lambda: gpd.read_parquet('f.parquet')):.3f} s; "
      f"bbox read of the sorted, covered file {timed_read(lambda: gpd.read_parquet('queryable.parquet', bbox=box)):.3f} s")
print(f"mean vertices per building: {shapely.get_num_coordinates(gdf.geometry.values).mean():.1f}")
try:
    gpd.read_parquet("f.parquet", bbox=box)
except ValueError as exc:
    print("without it:", str(exc).splitlines()[0])

one = sorted(glob.glob("buildings_parts/tile=*"), key=lambda d: -os.path.getsize(d + "/part.parquet"))[0]
t_one = timed_read(lambda: gpd.read_parquet(one + "/part.parquet"))
t_all = timed_read(lambda: [gpd.read_parquet(p) for p in glob.glob("buildings_parts/tile=*/part.parquet")])
print(f"read the largest tile ({len(gpd.read_parquet(one + '/part.parquet')):,} rows) {t_one:.3f} s; "
      f"every tile ({len(gdf):,} rows) {t_all:.3f} s")

63,517 buildings
GeoParquet               5.69 MB     read 0.05 s
GeoPackage              17.33 MB     read 0.11 s
GeoJSON                 24.25 MB     read 0.44 s
Shapefile (all parts)   32.60 MB     read 0.12 s

covering bbox: 8.09 MB against 5.69 MB (x1.42)
local read: whole plain file 0.045 s; bbox read of the sorted, covered file 0.021 s
mean vertices per building: 6.9
without it: Specifying 'bbox' not supported for this Parquet file (it should either have a bbox covering column or use 'point' encoding).
read the largest tile (8,890 rows) 0.021 s; every tile (63,517 rows) 0.444 s

Row-group statistics and why sorting comes first

Why one row group defeats pushdown

Parquet stores min and max for each column in each row group. A reader evaluating population > 10000 reads those statistics, skips groups whose maximum is below the threshold, and decodes the rest.

With one row group there is nothing to skip. The filter still works — it is applied after decoding — but it saves no I/O.

The default row-group size of around a million rows means most single-file GeoParquet has exactly one group, and every filter reads everything. Setting row_group_size to something like 50,000 is what turns pushdown on.

Why spatial sorting is the precondition

Statistics-based skipping works when a group's range is narrow. If rows are in arbitrary order, every group's bounding box covers most of the dataset, and no query can exclude any group.

Sorting by a space-filling curve — Hilbert or Morton — puts spatially close rows next to each other, so each row group covers a compact region. Then a bounding-box query genuinely excludes most groups.

Sorting is therefore not an optimisation to add later. Without it, the bbox column and small row groups cost size and buy nothing.

Why the bbox column costs 42%

A geometry column stores WKB; the covering bbox adds four doubles per row — 32 bytes before compression.

For small geometries like building footprints, averaging 6.9 vertices here, that is a large fraction of the row. Measured: 5.69 MB to 8.09 MB.

For large geometries — administrative boundaries with thousands of vertices — the same 32 bytes is negligible, and the feature is nearly free. The overhead is worst exactly where the row count is highest.

Why partitioning and row groups are different tools

Partitioning prunes at the file level using the directory name, so a reader skips files without opening them. It suits low-cardinality keys with even sizes: a region, a year, a tile.

Row groups prune within a file using statistics, so the file is opened and its metadata read. It suits high-cardinality or continuous values, and it is finer-grained.

Use both: partition by the coarse key you always filter on, and set row groups within each partition for everything else.

num_row_groups, partition_cols and loud fallbacks

• Check num_row_groups. One group means no pushdown.
• Sort spatially before writing, or the statistics are useless.
• write_covering_bbox=True is required for bbox= on read, and costs 42% here.
• geopandas.to_parquet has no partition_cols. Write partitions with a groupby.
• Even partition sizes matter more than the partitioning scheme.
• Very small partitions cost more than they save — aim for tens of megabytes.
• Column selection is free pushdown. Reading three columns of thirty reads a tenth of the file.
• Make fallbacks loud, or "pushdown" quietly becomes "read everything".

Internal links

• GeoParquet and columnar storage explained — why the format is fast
• Cloud-native geospatial explained — the pattern this belongs to
• How to scale a GeoPandas job with dask-geopandas — reading partitions in parallel
• How to read spatial data from S3 and other object storage — serving partitions remotely
• GIS vector file formats compared — the size comparison in context
• How to choose chunk and tile sizes that actually help — the raster equivalent
• How to reduce GIS file size in Python without wrecking the data — other size levers
• How to process a very large GeoPackage in chunks with Python — when the format is not columnar

FAQ

Why does bbox= fail on my GeoParquet file?

The file has no covering-bbox column. Rewrite it with write_covering_bbox=True; it cost 42% more space for small building footprints.

Why is my filter not any faster?

Almost certainly one row group. Check pq.ParquetFile(path).metadata.num_row_groups and set row_group_size when writing.

Do I need to sort my data before writing?

Yes, for spatial pushdown to work. Without a spatial sort every row group's bounding box covers the whole dataset and none can be skipped.

How small should row groups be?

Tens of thousands of rows is a good starting point. Small enough that statistics discriminate, large enough that per-group overhead stays low.

How do I write partitioned GeoParquet?

geopandas.to_parquet has no partition_cols. Group by the key and write one file per partition into a key=value directory, or use dask_geopandas.

Is GeoParquet really that much smaller?

Measured on 63,517 building polygons: 5.69 MB against 17.33 MB for GeoPackage, 24.25 MB for GeoJSON and 32.60 MB for a shapefile.

Should I partition or use row groups?

Both. Partition on the coarse key you always filter on; use row groups within each partition for everything else.