GeoParquet and Columnar Storage Explained: Why It Reads So Much Faster

Problem statement

The same 232,717 census blocks — every 2020 census block in Georgia, 17 attribute columns and a polygon each — saved in four formats (the setup is under The data used on this page):

import time
import geopandas as gpd

for path in ["blocks.shp", "blocks.gpkg", "blocks.geojson", "blocks.parquet"]:
    times = []
    for _ in range(2):                                   # the faster of two reads
        t0 = time.perf_counter()
        gdf = gpd.read_parquet(path) if path.endswith("parquet") else gpd.read_file(path)
        times.append(time.perf_counter() - t0)
    print(f"{path:<16} {min(times):>6.2f} s")

blocks.shp         1.73 s
blocks.gpkg        1.77 s
blocks.geojson    15.52 s
blocks.parquet     0.62 s

And then this:

import pandas as pd

t0 = time.perf_counter()
gdf = gpd.read_parquet("blocks.parquet", columns=["GEOID20", "POP20", "geometry"])
print(f"three columns with geometry:    {time.perf_counter() - t0:.2f} s")

t0 = time.perf_counter()
df = pd.read_parquet("blocks.parquet", columns=["GEOID20", "POP20", "HOUSING20"])
print(f"three columns without geometry: {time.perf_counter() - t0:.3f} s")

three columns with geometry:    0.59 s
three columns without geometry: 0.021 s

The second read never touches the geometry column, which is most of the file, so it skips most of the bytes.

That last result is not an optimisation applied to a row-based format — it is impossible in one. The difference is structural, and understanding it tells you when GeoParquet is the right choice and when it is not.

Quick answer

import geopandas as gpd

gdf = gpd.read_file("blocks.gpkg")
gdf.to_parquet("blocks.parquet", compression="zstd")      # write once

gdf = gpd.read_parquet("blocks.parquet")                  # read many times
gdf = gpd.read_parquet("blocks.parquet", columns=["GEOID20", "POP20", "geometry"])

All from the runs on this page. GeoParquet's advantage on a full read is real but modest for polygon-heavy data, because geometry is nearly all the bytes and every format has to turn it into Shapely objects. It is overwhelming when you can leave the geometry out.

Use GeoParquet for analytical reads and GeoPackage for editing and interchange. They are not competitors; they answer different questions.

What columnar means on disk, and writing it properly

1. Understand what "columnar" means physically

A table has to be laid out linearly on disk. There are two ways to do it.

Row-oriented — everything about record 1, then everything about record 2:

[GEOID20=130010001001000|UR20=R|POP20=23|geom=...] [GEOID20=130010001001001|UR20=R|POP20=0|geom=...] ...

Reading POP20 alone means seeking past every other field of every row, or more realistically reading whole pages and discarding most of them. Shapefile, GeoPackage and most databases work this way, and it is the right choice when you routinely want whole records.

Column-oriented — every GEOID20, then every UR20, then every POP20:

[GEOID20: 130010001001000, 130010001001001, ...] [UR20: R, R, U, ...] [POP20: 23, 0, 41, ...] [geom: WKB, ...]

Reading POP20 alone means reading one contiguous run of bytes and nothing else. That is the whole trick, and everything else follows from it.

2. Follow the consequences

Column pruning. Read only the columns you name — as long as geometry is not one of them, which is the catch measured above.

Better compression. A column holds one type with repetitive values. A column of county codes drawn from Georgia's 159 counties compresses to almost nothing by dictionary encoding; in a row layout those values are separated by unrelated bytes, so the compressor cannot see the pattern. The file's own metadata shows it:

import pyarrow.parquet as pq

meta = pq.read_metadata("blocks.parquet")
rg = meta.row_group(0)
print(f"{meta.num_row_groups} row group(s), {rg.num_rows:,} rows")
raw = {"COUNTYFP20": 3, "UR20": 1, "GEOID20": 15, "INTPTLAT20": 11}   # characters per value as text
for j in range(rg.num_columns):
    col = rg.column(j)
    name = col.path_in_schema
    if name in ("COUNTYFP20", "UR20", "GEOID20", "INTPTLAT20", "POP20", "geometry"):
        per_row = col.total_compressed_size / rg.num_rows
        text = f"   (as text: {raw[name]} bytes)" if name in raw else ""
        print(f"  {name:<11} {col.total_compressed_size / 1e6:>6.1f} MB  {per_row:>8.2f} bytes a row{text}")

1 row group(s), 232,717 rows
  COUNTYFP20     0.2 MB      0.88 bytes a row   (as text: 3 bytes)
  GEOID20        1.9 MB      8.11 bytes a row   (as text: 15 bytes)
  UR20           0.0 MB      0.15 bytes a row   (as text: 1 bytes)
  INTPTLAT20     2.0 MB      8.73 bytes a row   (as text: 11 bytes)
  POP20          0.3 MB      1.35 bytes a row
  geometry     251.8 MB   1082.07 bytes a row

The county code is three characters as text and under a byte a row here, because dictionary encoding stores each of the 159 codes once and then a small index per row, and run-length encoding collapses runs of the same index. Pyarrow applies that encoding before the compressor, so the metadata's own "compression ratio" for these columns reads close to 1× — the saving has already happened. geometry is the opposite case: WKB is already dense binary, it barely compresses, and it is nearly all of the file.

Predicate pushdown. Each row group stores per-column min and max, so a reader can skip an entire group without decompressing it:

rural = gpd.read_parquet("blocks.parquet", filters=[("UR20", "==", "R")])
print(f"{len(rural):,} rural blocks")

122,668 rural blocks

Typed, vectorised reads. A column of int64 is already a contiguous typed array, so it becomes a NumPy array with a memory copy rather than a per-value parse.

3. Write GeoParquet properly

blocks.to_parquet(
    "blocks_zstd.parquet",
    compression="zstd",            # smaller than the snappy default
    compression_level=3,
    row_group_size=10_000,         # rows per group — the pushdown granularity
    write_covering_bbox=True,      # per-row bbox columns for spatial filtering
    geometry_encoding="WKB",       # the interoperable default
    index=False,
)
import os
for path in ["blocks.parquet", "blocks_zstd.parquet"]:
    print(f"{path:<20} {os.path.getsize(path) / 1e6:>6.0f} MB, "
          f"{pq.read_metadata(path).num_row_groups} row groups")

blocks.parquet          263 MB, 1 row groups
blocks_zstd.parquet     241 MB, 24 row groups

Row-group size is the one real tuning decision, and the default is usually wrong for spatial filtering: pyarrow puts up to about a million rows in a group, so this whole file was written as a single group, and no statistic can skip part of one group. Pick a size that gives tens of groups or more — 10,000 rows here, 24 groups.

4. Use spatial filtering

Parquet has no R-tree, but bounding-box statistics per row group achieve much of the same. What decides whether they help is the order the rows are in — so write three copies: shuffled, as the census file came (sorted by block code), and sorted along a Hilbert curve:

import numpy as np


def groups_touched(path, bbox):
    """How many row groups a bbox read cannot rule out, from the bbox column statistics."""
    meta, (xmin, ymin, xmax, ymax), n = pq.read_metadata(path), bbox, 0
    names = [meta.schema.column(j).path for j in range(meta.num_columns)]
    for i in range(meta.num_row_groups):
        st = {names[j]: meta.row_group(i).column(j).statistics for j in range(meta.num_columns)}
        n += not (st["bbox.xmin"].min > xmax or st["bbox.xmax"].max < xmin or
                  st["bbox.ymin"].min > ymax or st["bbox.ymax"].max < ymin)
    return n


orders = {
    "shuffled": blocks.sample(frac=1, random_state=0),
    "census order": blocks,
    "hilbert": blocks.iloc[blocks.hilbert_distance().argsort()],
}
atlanta = (1_060_000, 1_247_000, 1_070_000, 1_257_000)      # 10 km square, central Atlanta
for name, gdf in orders.items():
    path = f"blocks_{name.replace(' ', '_')}.parquet"
    gdf.to_parquet(path, row_group_size=10_000, write_covering_bbox=True, index=False)
    t0 = time.perf_counter()
    subset = gpd.read_parquet(path, bbox=atlanta)
    print(f"{name:<13} {time.perf_counter() - t0:.3f} s  {len(subset):,} blocks  "
          f"{groups_touched(path, atlanta):>2} of {pq.read_metadata(path).num_row_groups} row groups read")

shuffled      0.246 s  3,179 blocks  24 of 24 row groups read
census order  0.266 s  3,179 blocks  24 of 24 row groups read
hilbert       0.062 s  3,179 blocks   5 of 24 row groups read

The saving depends entirely on spatial locality within row groups. If the data is written in a spatially coherent order, a bounding box intersects few groups and most of the file is skipped. If the rows are in random order, every group's bbox covers the whole state and nothing can be skipped. The order the census file arrives in turned out to be no better than shuffled for this query: each 10,000-row group still spans enough of the state that none could be ruled out. Only the Hilbert sort let the reader skip most of the file.

The sort pays a second time, in file size. Neighbouring census blocks share their boundaries, so in Hilbert order the same coordinates recur a few rows apart — close enough for zstd, which looks back over a long window, to store them once:

for name, gdf in orders.items():
    for codec in ["snappy", "zstd"]:
        gdf.to_parquet("size_test.parquet", compression=codec, row_group_size=10_000, index=False)
        print(f"{name:<13} {codec:<7} {os.path.getsize('size_test.parquet') / 1e6:>4.0f} MB")

shuffled      snappy   266 MB
shuffled      zstd     238 MB
census order  snappy   265 MB
census order  zstd     237 MB
hilbert       snappy   241 MB
hilbert       zstd     157 MB

hilbert_distance() maps 2-D position to a 1-D ordering that preserves proximity, so spatially close features land in the same row group. It is a one-off cost at write time and it is what makes bbox= reads worthwhile. Note it returns a Series of distances to sort by — gdf.sort_values(by=gdf.hilbert_distance()) raises, because by= takes column names; iloc[... .argsort()] is the working form.

5. Know when GeoParquet is the wrong choice

Parquet is immutable and analytical. It has no support for:

• Updating a row. You rewrite the file, or write a new one.
• Appending efficiently. Appending means a new file; readers then treat a directory of files as one dataset.
• Random access by key. No index — finding one id means scanning, or relying on statistics.
• Editing in a desktop GIS. QGIS reads GeoParquet through GDAL's Parquet driver but is not an editing environment for it.
• Multiple layers in one file. One Parquet file is one table.

# ❌ this rewrites the whole file to change one value
gdf.loc[gdf["GEOID20"] == "130010001001000", "POP20"] = 24
gdf.to_parquet("blocks.parquet")

# ✅ mutable data belongs in GeoPackage or PostGIS

The practical pattern is GeoPackage or PostGIS as the source of record, GeoParquet as the analytical copy — written once when the source changes, read many times.

Georgia's census blocks as GeoParquet against row formats

The data used on this page

Every 2020 census block in Georgia, from the US Census Bureau's TIGER/Line files: 232,717 polygons with 17 attribute columns, reprojected to the Albers equal-area projection for the contiguous US (EPSG:5070) so coordinates and areas are in metres. The setup writes the four copies the problem statement reads.

import os
import time

import geopandas as gpd
import requests

URL = "https://www2.census.gov/geo/tiger/TIGER2020/TABBLOCK20/tl_2020_13_tabblock20.zip"
with open("tl_2020_13_tabblock20.zip", "wb") as f:
    f.write(requests.get(URL, timeout=600).content)
blocks = gpd.read_file("tl_2020_13_tabblock20.zip").to_crs(5070)

blocks.to_file("blocks.shp")
blocks.to_file("blocks.gpkg")
blocks.to_file("blocks.geojson")
blocks.to_parquet("blocks.parquet")
print(f"{len(blocks):,} blocks, {len(blocks.columns)} columns, EPSG:{blocks.crs.to_epsg()}")

232,717 blocks, 18 columns, EPSG:5070

Example 1: measuring the difference on your own data

import os
import time
from pathlib import Path

import geopandas as gpd


def compare_formats(gdf, out_dir="fmt_test", columns=None):
    out = Path(out_dir)
    out.mkdir(exist_ok=True)
    columns = columns or [c for c in gdf.columns if c != gdf.geometry.name][:3]

    targets = [
        ("GeoPackage", out / "t.gpkg", lambda p: gdf.to_file(p, driver="GPKG")),
        ("GeoJSON", out / "t.geojson", lambda p: gdf.to_file(p, driver="GeoJSON")),
        ("Parquet snappy", out / "t_snappy.parquet",
         lambda p: gdf.to_parquet(p, compression="snappy", write_covering_bbox=True, row_group_size=10_000)),
        ("Parquet zstd", out / "t_zstd.parquet",
         lambda p: gdf.to_parquet(p, compression="zstd", compression_level=3,
                                  write_covering_bbox=True, row_group_size=10_000)),
        ("FlatGeobuf", out / "t.fgb", lambda p: gdf.to_file(p, driver="FlatGeobuf")),
    ]

    print(f"{len(gdf):,} rows, {len(gdf.columns)} columns\n")
    print(f"{'format':<16}{'write':>9}{'size':>10}{'read all':>10}"
          f"{'read 3 cols':>13}{'bbox read':>11}")

    minx, miny, maxx, maxy = gdf.total_bounds
    dx, dy = (maxx - minx) * 0.1, (maxy - miny) * 0.1
    bbox = (minx + dx, miny + dy, minx + 3 * dx, miny + 3 * dy)

    for name, path, write in targets:
        t0 = time.perf_counter(); write(path); w = time.perf_counter() - t0
        size = os.path.getsize(path) / 1e6

        is_pq = path.suffix == ".parquet"
        reader = gpd.read_parquet if is_pq else gpd.read_file

        t0 = time.perf_counter(); reader(path); r_all = time.perf_counter() - t0
        t0 = time.perf_counter(); reader(path, columns=[*columns, "geometry"] if is_pq else columns)
        r_cols = time.perf_counter() - t0
        t0 = time.perf_counter(); reader(path, bbox=bbox)
        r_bbox = time.perf_counter() - t0

        print(f"{name:<16}{w:>8.2f}s{size:>9.0f}M{r_all:>9.2f}s{r_cols:>12.2f}s{r_bbox:>10.2f}s")


compare_formats(blocks)

232,717 rows, 18 columns

format              write      size  read all  read 3 cols  bbox read
GeoPackage          2.19s      360M     1.67s        0.95s      0.04s
GeoJSON            23.80s      940M    14.77s       14.50s     13.36s
Parquet snappy      1.22s      274M     0.73s        0.67s      0.28s
Parquet zstd        1.71s      241M     0.72s        0.64s      0.27s
FlatGeobuf          2.62s      333M     1.65s        0.96s      0.04s

Four readings worth pulling out. GeoJSON is an order of magnitude slower at everything, and naming columns or a bbox does not help it, because the whole file has to be parsed either way. Parquet reads everything fastest, but only by two to three times here, because geometry dominates. GeoPackage and FlatGeobuf have real spatial indexes and win the bbox read outright — hundredths of a second, against a quarter of a second for Parquet written in the census file's order — so for a workload that is only ever spatial filtering, they beat an unsorted Parquet file. And Parquet writes fastest, so the convert-once cost is small.

Much of GeoPackage's gap is GDAL's default read path rather than the format. Reading through Arrow, and leaving the geometry out, change the picture:

for use_arrow in (False, True):
    t0 = time.perf_counter()
    gpd.read_file("blocks.gpkg", use_arrow=use_arrow)
    print(f"GeoPackage, use_arrow={use_arrow!s:<5}         {time.perf_counter() - t0:.2f} s")
t0 = time.perf_counter()
gpd.read_file("blocks.gpkg", columns=["GEOID20", "POP20", "HOUSING20"], ignore_geometry=True)
print(f"GeoPackage, 3 columns, no geometry   {time.perf_counter() - t0:.2f} s")
t0 = time.perf_counter()
pd.read_parquet("blocks.parquet", columns=["GEOID20", "POP20", "HOUSING20"])
print(f"GeoParquet, 3 columns, no geometry   {time.perf_counter() - t0:.3f} s")

GeoPackage, use_arrow=False         2.24 s
GeoPackage, use_arrow=True          1.14 s
GeoPackage, 3 columns, no geometry   0.36 s
GeoParquet, 3 columns, no geometry   0.014 s

With Arrow, GeoPackage's full read roughly halves, closing most of the gap to GeoParquet. The attribute-only read is where the layouts part company: GeoPackage still walks every row to pick the three fields out, and Parquet reads three contiguous columns.

Example 2: a conversion utility with sensible defaults

from pathlib import Path
import time
import geopandas as gpd


def to_geoparquet(src, dst=None, *, sort_spatially=True, row_group_size=10_000,
                  compression="zstd", drop_columns=None, target_crs=None):
    """Convert any GeoPandas-readable file to a well-formed GeoParquet."""
    src = Path(src)
    dst = Path(dst) if dst else src.with_suffix(".parquet")

    t0 = time.perf_counter()
    gdf = gpd.read_file(src)
    read_s = time.perf_counter() - t0

    if drop_columns:
        gdf = gdf.drop(columns=[c for c in drop_columns if c in gdf.columns])
    if target_crs and gdf.crs != target_crs:
        gdf = gdf.to_crs(target_crs)
    if gdf.crs is None:
        raise ValueError(f"{src} has no CRS — GeoParquet requires one in its metadata")

    if sort_spatially and len(gdf) > row_group_size:
        t0 = time.perf_counter()
        gdf = gdf.iloc[gdf.hilbert_distance().argsort()].reset_index(drop=True)
        sort_s = time.perf_counter() - t0
    else:
        sort_s = 0.0

    t0 = time.perf_counter()
    tmp = dst.with_suffix(".tmp.parquet")
    gdf.to_parquet(tmp, compression=compression, compression_level=3,
                   row_group_size=row_group_size, write_covering_bbox=True,
                   index=False)
    tmp.replace(dst)
    write_s = time.perf_counter() - t0

    src_mb = src.stat().st_size / 1e6
    dst_mb = dst.stat().st_size / 1e6
    print(f"{src.name} → {dst.name}")
    print(f"  {len(gdf):,} rows, {len(gdf.columns)} columns, {gdf.crs}")
    print(f"  read {read_s:.1f}s  sort {sort_s:.1f}s  write {write_s:.1f}s")
    print(f"  {src_mb:,.0f} MB → {dst_mb:,.0f} MB ({100 * (1 - dst_mb/src_mb):.0f}% smaller)")
    return dst


to_geoparquet("blocks.gpkg", "blocks_clean.parquet", target_crs=5070,
              drop_columns=["GEOIDFQ20", "INTPTLAT20", "INTPTLON20"])

blocks.gpkg → blocks_clean.parquet
  232,717 rows, 15 columns, EPSG:5070
  read 1.8s  sort 0.1s  write 1.8s
  360 MB → 145 MB (60% smaller)

The Hilbert sort is a one-off cost and is what makes every later bbox= read fast. Without it, spatially adjacent features can be scattered across row groups and the bbox statistics cannot exclude anything.

Requiring a CRS is deliberate: GeoParquet stores it in the file's metadata, and a file written without one is valid Parquet but not valid GeoParquet — readers then have no way to interpret the coordinates.

The temp-and-rename means an interrupted conversion leaves no half-written file to be mistaken for a complete one.

Example 3: reading larger-than-memory data

Parquet's structure allows things a single-file read cannot:

from pathlib import Path

import pyarrow as pa
import pyarrow.parquet as pq
import pyarrow.dataset as ds
import geopandas as gpd

# ── (a) one row group at a time ─────────────────────────────────────────────
def iter_row_groups(path, columns=None):
    pf = pq.ParquetFile(path)
    for i in range(pf.num_row_groups):
        table = pf.read_row_group(i, columns=columns)
        yield gpd.GeoDataFrame.from_arrow(table)

total_area, rows = 0.0, 0
for chunk in iter_row_groups("blocks_clean.parquet", columns=["GEOID20", "geometry"]):
    total_area += chunk.area.sum()
    rows += len(chunk)
print(f"{rows:,} rows, {total_area / 1e6:,.0f} km², peak memory bounded by one group")

# ── (b) partitioned writing, so future reads skip whole directories ─────────
# GeoDataFrame.to_parquet writes one file, so write one per county yourself
for county, part in blocks.groupby("COUNTYFP20"):
    folder = Path(f"blocks_by_county/COUNTYFP20={county}")
    folder.mkdir(parents=True, exist_ok=True)
    part.drop(columns="COUNTYFP20").to_parquet(folder / "part-0.parquet", compression="zstd",
                                               write_covering_bbox=True)

# ── (c) a directory of files as one dataset ─────────────────────────────────
# declare the partition column a string, or "121" is read back as the integer 121
by_county = ds.partitioning(pa.schema([("COUNTYFP20", pa.string())]), flavor="hive")
fulton = gpd.read_parquet("blocks_by_county/", partitioning=by_county,
                          filters=[("COUNTYFP20", "==", "121"), ("UR20", "==", "U")])   # Fulton County, urban
files = ds.dataset("blocks_by_county/", format="parquet", partitioning=by_county).files
print(f"{len(fulton):,} rows read from a {len(files)}-file dataset")

232,717 rows, 153,905 km², peak memory bounded by one group
11,384 rows read from a 159-file dataset

Two things in (b) and (c) are not what the pyarrow documentation leads you to expect. GeoDataFrame.to_parquet has no partition_cols — passing one fails with IsADirectoryError — so the loop writes the Hive layout itself. And Hive partitioning infers the type of the directory values, so COUNTYFP20=121 comes back as int32, the string filter fails with ArrowNotImplementedError: Function 'equal' has no kernel matching input types (int32, string), and county 001 would silently become 1. Declaring the partition schema fixes both.

Partitioning by a column writes blocks_by_county/COUNTYFP20=121/part-0.parquet, and a filter on COUNTYFP20 then skips entire directories without opening a file. That is a different mechanism from row-group statistics and much more decisive, but it only helps for the column you partitioned on — so partition by whatever you filter on most.

Per-row-group iteration is the memory-bounded read: peak usage is one group, so a dataset far larger than memory is no harder than a small one. The same reasoning as processing a large GeoPackage in chunks, with the advantage that the chunk boundaries are already in the file.

Physical layout explains every difference

The performance difference between GeoParquet and a row-based format is not an implementation detail — it follows from a decision about physical layout, and every advantage and limitation traces back to it.

Storage is linear; a table is two-dimensional. Something has to be chosen, and the choice determines what is cheap. Row-oriented layout keeps a record's fields adjacent, so reading or writing a whole record touches one contiguous region — excellent for transactional work, where you fetch a record by key, edit it and write it back. Column-oriented layout keeps a column's values adjacent, so reading one column touches one contiguous region — excellent for analytical work, where you compute over a few columns across many rows.

Compression is where the difference compounds, and it is not obvious in advance. Compressors work by finding repetition. Georgia's 232,717 county codes are drawn from 159 values, and dictionary encoding reduces them to a lookup table plus a small index per row — under a byte a row, measured earlier, against three as text. The same values in a row layout are separated by unrelated bytes of other columns, so the pattern is invisible within any compression window. Columnar storage does not merely compress better; it makes a different class of encoding possible.

Parsing cost disappears for typed columns. A float64 column in Parquet is already a contiguous run of IEEE doubles, so becoming a NumPy array is a memory copy. In GeoJSON the same value is the text "-2.24261084" and must be converted digit by digit — which is why GeoJSON read about 24 times slower than Parquet on the blocks while being under four times the size: the gap is mostly parsing, not bytes.

Geometry is the exception, and it is worth knowing. WKB is already dense binary, so it barely compresses and it still has to be parsed into Shapely objects on read. In the measurements above, geometry is 252 MB of the 263 MB file — 1,082 bytes a row, against under a byte for the county code. So GeoParquet's advantage is largest on attribute-heavy data and smallest on geometry-heavy data with few attributes — and the columns= win is largest precisely when you can leave the geometry column out entirely.

The absence of a spatial index is the honest trade-off. GeoPackage carries an R-tree, so a bbox query goes straight to the matching rows regardless of storage order. Parquet has only per-row-group min/max statistics, which can exclude a group but cannot locate a row. This works well when the data is spatially sorted — a Hilbert ordering makes each group cover a compact region, so a small bbox intersects few groups — and not at all when rows are in arbitrary order, since then every group's bbox spans the whole extent. That is why the sort is not an optional refinement but the thing that makes spatial filtering work.

Finally, immutability is a design decision rather than a missing feature. Parquet's compression and encoding are computed per row group over the whole group, so changing one value requires re-encoding the group and rewriting the file. That is what buys the compression ratio and the read speed. It also means Parquet is the wrong home for data people edit — the natural pattern is a mutable source of record in GeoPackage or PostGIS, with GeoParquet as a derived, regenerated analytical copy.

CRS metadata, column skipping and spatial sorting

• GeoParquet requires a CRS in its metadata. Writing a frame with crs=None produces valid Parquet that is not valid GeoParquet.
• gpd.read_parquet accepts columns=; gpd.read_file does too but cannot skip them physically.
• bbox= only helps if the data is spatially sorted. Use hilbert_distance() before writing.
• write_covering_bbox=True adds bbox columns and is what makes bbox= reads possible at all.
• Row-group size is the main tuning knob. Pyarrow's default puts up to about a million rows in one group, which leaves nothing for statistics to skip in a file this size. Aim for tens of groups or more.
• zstd was about 10% smaller than snappy on the blocks in their original order and 35% smaller once Hilbert-sorted; snappy has the widest reader support.
• geometry_encoding="geoarrow" is faster than WKB and less widely readable. Prefer WKB for interchange.
• Partitioning by a column lets readers skip whole directories, which beats row-group statistics — but only for that column.
• QGIS reads GeoParquet when its GDAL has the Parquet driver (GDAL 3.5 or later), and can export to it; it cannot edit a Parquet file in place.
• A directory of Parquet files is one dataset to pyarrow.dataset, which is how appends are normally handled.

Internal links

• GIS vector file formats compared — where GeoParquet sits among the alternatives
• Why GeoPandas is slow: the four real bottlenecks — I/O as one of the four
• How to reduce GIS file size in Python — compression and simplification together
• How to process a very large GeoPackage in chunks — the row-group iteration idea, on GeoPackage
• Fixing memory errors in GeoPandas when working with large files — what bounded reads solve
• PostGIS explained: when a spatial database beats a folder of files — the mutable alternative
• WKT, WKB and GeoJSON: how geometry is actually stored — what is inside the geometry column
• How to read a large PostGIS table into Python — streaming from the other direction
• M values explained: geometry that carries a measure — the ordinate only Parquet keeps here

FAQ

Why is GeoParquet so much faster to read?

It only reads the columns you ask for, its columnar layout compresses attribute columns far better, and typed columns become NumPy arrays by memory copy rather than by parsing text. For polygon-heavy data the full-read gain is modest — geometry is most of the bytes in any format — and the large gains come from reading attributes without geometry.

Should I replace GeoPackage with GeoParquet?

Not replace — complement. GeoPackage for editing, interchange and desktop GIS; GeoParquet as a derived analytical copy that is regenerated when the source changes.

Does GeoParquet support spatial indexing?

Not an R-tree. It stores per-row-group bounding boxes, which can skip groups. That only works if the data is spatially sorted, so sort by hilbert_distance() before writing.

What row-group size should I use?

Enough rows that per-group metadata is negligible, few enough that the file has tens of groups or more — 10,000 gave 24 groups for Georgia's 232,717 blocks. Smaller gives finer filtering and more metadata overhead; larger gives better compression and coarser skipping. The pyarrow default of about a million rows is too coarse for spatial filtering on most single files.

Can I append to a GeoParquet file?

Not to one file. Write a new file into the same directory and read the directory as a dataset with pyarrow.dataset, which is the standard pattern.

Which compression should I choose?

zstd for the best size, snappy for the widest reader compatibility. Neither changes the data; both are lossless. zstd gains most on spatially sorted data, where shared boundaries repeat nearby.

Can QGIS open GeoParquet?

Yes, when its GDAL (3.5 or later) is built with the Arrow/Parquet driver, and it can export to GeoParquet through the same driver. A Parquet file cannot be edited in place, though, which is another reason to keep an editable GeoPackage alongside.