⚙️ How-To DuckDB & analytic SQL Files & formats DuckDB for Spatial Data · step 4 of 20 Updated

How to Read Shapefiles, GeoJSON and GeoParquet in DuckDB

Problem statement

The data is already on disk: a folder of shapefiles, a GeoPackage from a colleague, a GeoParquet export, a CSV of coordinates. The usual next step is to load all of it into GeoPandas, which means holding it all in memory before you know which parts you need.

DuckDB reads these formats in place and lets the query decide what to load. That changes the shape of the work: instead of read_file() then filter, you filter in the query and only the result crosses into Python.

It also has two readers with different characteristics, and choosing the wrong one is the usual first mistake. ST_Read goes through GDAL and handles everything GDAL handles; read_parquet is DuckDB's native columnar reader and is far faster, but only for Parquet.

Quick answer 

import duckdb

con = duckdb.connect()
con.execute("install spatial; load spatial;")

# GDAL-backed: shapefile, GeoPackage, GeoJSON, KML, GML, and everything else
con.execute("select count(*) from st_read('provinces.shp')").fetchone()
con.execute("select count(*) from st_read('data.gpkg', layer='roads')").fetchone()

# native columnar: Parquet and GeoParquet
con.execute("select count(*) from read_parquet('places.parquet')").fetchone()

# plain tabular with coordinate columns
con.execute("""
    select st_point(lon, lat) as geom, name
    from read_csv('points.csv')
""").df()

Measured on Natural Earth: ST_Read on 7,342 populated places took 0.132 s against GeoPandas' 0.276 s; on 4,596 provinces, 0.162 s against 0.210 s.

Decision diagram routing three file kinds to read_parquet, st_read or read_csv.
ST_Read gives universal format support; read_parquet gives pruning.

Step-by-step solution

1. Use ST_Read for anything GDAL supports

ST_Read is a table function: it appears in the FROM clause and yields rows with a geom column plus the source's attributes.

select * from st_read('boundaries.geojson') limit 5;
select * from st_read('data.gpkg', layer = 'roads');
select * from st_read('archive.zip/inner/file.shp');       -- GDAL vsizip paths work

To see what a file contains before reading it:

select * from st_read_meta('data.gpkg');

which lists the layers, their geometry types, feature counts and CRS.

2. Use read_parquet for Parquet, always

Parquet through DuckDB's native reader is a different class of fast, because it can skip columns and row groups. Measured on 13,464,017 rows:

                                     TSV       Parquet
count(distinct country)             1.30 s      0.03 s
group by floor(lon), floor(lat)     1.01 s      0.06 s

GeoParquet is Parquet with a WKB geometry column and a geo metadata key, so read_parquet reads it. Convert the WKB to a geometry when you need spatial functions:

select st_geomfromwkb(geometry) as geom, name
from read_parquet('places.parquet');

3. Read many files at once with a glob

Both readers accept wildcards, which replaces the usual "loop over the folder and concatenate" pattern:

select count(*) from read_parquet('data/*.parquet');
select count(*) from read_parquet('data/**/*.parquet');
select county, count(*) from st_read('counties/*.shp') group by 1;

read_parquet also understands Hive-style directories, turning path components into columns:

select country, count(*) from read_parquet('by_country/**/*.parquet',
                                           hive_partitioning = 1)
group by 1;

4. Filter in the query, not afterwards

The point of reading in place is that the query decides what is materialised:

select name, geom
from st_read('provinces.shp')
where admin = 'United States of America';

For Parquet this pushes down to column and row-group pruning; for ST_Read the filter still happens after GDAL has produced the rows, but nothing is copied into Python that the query did not ask for.

5. Watch what happens to the CRS

ST_Read attaches the source CRS to the column type:

select typeof(geom) from st_read('provinces.shp') limit 1;
-- GEOMETRY('EPSG:4326')

Two things to know. A shapefile with no .prj yields a plain GEOMETRY with no CRS and no warning. And ST_Transform returns a plain GEOMETRY too, so the CRS is dropped by the transform that changed it.

6. Write back with the right driver

COPY writes Parquet natively; the GDAL driver writes everything else:

copy places to 'places.parquet' (format parquet, compression zstd);
copy places to 'places.gpkg' with (format gdal, driver 'GPKG');
copy places to 'places.geojson' with (format gdal, driver 'GeoJSON', srs 'EPSG:4326');

A plain COPY ... TO 'x.csv' writes the geometry as WKT text, which round-trips as text and loses the CRS. That is fine as an interchange format and wrong as an archive.

Bar chart of read times for two shapefiles in DuckDB ST_Read and GeoPandas.
Both go through GDAL; the difference is what happens to the rows afterwards.

Code examples

Example 1 — inspecting a file before reading it 

import duckdb


def describe_spatial_file(path, con=None):
    """Layers, geometry types, counts and CRS — without loading the features."""
    con = con or duckdb.connect()
    con.execute("load spatial")

    meta = con.execute(f"select * from st_read_meta('{path}')").df()
    for layer in meta["layers"][0]:
        fields = layer.get("geometry_fields", [{}])[0]
        crs = fields.get("crs") or {}
        print(f"layer {layer['name']!r}")
        print(f"  features   {layer.get('feature_count', '?')}")
        print(f"  geometry   {fields.get('type', '?')}")
        print(f"  crs        {crs.get('name', 'NONE')} "
              f"({crs.get('auth_name', '')}:{crs.get('auth_code', '')})")
        print(f"  attributes {len(layer.get('fields', []))}")
    return meta

Example 2 — a reader that picks the right function 

import os
import duckdb

GDAL_EXTENSIONS = {".shp", ".gpkg", ".geojson", ".json", ".kml", ".gml",
                   ".gdb", ".fgb", ".tab", ".mif"}


def read_spatial(con, path, layer=None, columns="*", where=None):
    """One entry point; the right reader chosen from the extension."""
    ext = os.path.splitext(path)[1].lower()

    if ext == ".parquet":
        source = f"read_parquet('{path}')"
    elif ext in (".csv", ".tsv", ".txt"):
        source = f"read_csv('{path}')"
    elif ext in GDAL_EXTENSIONS or "*" in path:
        args = f"'{path}'" + (f", layer='{layer}'" if layer else "")
        source = f"st_read({args})"
    else:
        raise ValueError(f"no reader for {ext!r}; try st_read and see what GDAL says")

    sql = f"select {columns} from {source}"
    if where:
        sql += f" where {where}"
    return con.execute(sql)

Example 3 — converting a folder of shapefiles into one GeoParquet file 

def shapefiles_to_geoparquet(con, pattern, target, sort_by=None,
                             compression="zstd"):
    """The conversion that makes every later query faster."""
    import os
    import time

    started = time.perf_counter()
    order = f"order by {sort_by}" if sort_by else ""
    con.execute(f"""
        copy (select * exclude geom, st_aswkb(geom) as geometry
              from st_read('{pattern}') {order})
        to '{target}' (format parquet, compression {compression})
    """)
    rows = con.execute(f"select count(*) from read_parquet('{target}')").fetchone()[0]
    print(f"{rows:,} features → {target} "
          f"({os.path.getsize(target) / 1e6:,.1f} MB) "
          f"in {time.perf_counter() - started:.1f}s")
    if sort_by:
        print(f"  sorted by {sort_by} — filters on it can skip row groups")

The sort_by argument is worth using. Measured on the same 13.5-million-row dataset, a filtered query read 0.26 MB from a sorted file and 10.25 MB from a shuffled one, and the sorted file was 40% smaller.

Explanation

Why ST_Read and read_parquet are different tools

ST_Read calls GDAL, which opens the file, decodes features one at a time and hands them to DuckDB. That gives universal format support and a per-feature cost: GDAL is doing the parsing, and DuckDB cannot push predicates into it beyond what the driver supports.

read_parquet is DuckDB's own reader operating on a format designed for exactly this. It knows the column layout and the row-group statistics, so it can read three columns of nineteen and skip whole groups whose statistics exclude the filter.

The rule of thumb: use ST_Read to ingest, read_parquet to work.

Why converting to GeoParquet is usually worth it once

The conversion is a single COPY, and it changes every subsequent query's economics. Measured on the same data: 1.79 GB of TSV became 446 MB of Parquet, and a grid aggregation went from 1.01 s to 0.06 s.

For a dataset that is queried more than about twice, the conversion pays for itself. For a one-off read of a shapefile somebody just sent, it does not.

Why the CRS behaviour needs watching

ST_Read reads the source's CRS and attaches it to the column type — GEOMETRY('EPSG:4326') — which is a genuine improvement over the older behaviour of a bare geometry.

But the type is not enforced across operations. Joining a 4326 column with a 3857 column returns zero rows and no error, and ST_Transform returns a geometry with no CRS at all. So the type is a helpful label rather than a guarantee, and the discipline of tracking CRS yourself still applies.

Why globs replace the folder loop

A folder of a hundred shapefiles is normally handled by a Python loop that reads, appends and concatenates — a hundred GeoDataFrames materialised so that one aggregate can be computed.

st_read('counties/*.shp') treats the whole folder as one table. The engine streams it, the aggregation happens as it goes, and only the result exists in memory. That is the same shift the rest of this article is about, applied to files instead of rows.

Grid of four DuckDB export targets and what happens to geometry and CRS in each.
A plain COPY to CSV succeeds and produces a file no GIS tool recognises.

Edge cases or notes

  • st_read_meta before st_read for an unfamiliar file — it lists layers, counts and CRS without reading features.
  • A shapefile without a .prj reads fine and yields a geometry with no CRS. No warning.
  • A shapefile without a .dbf still reads — measured, 7,342 rows came back with the geometry and no attributes.
  • COPY TO 'x.csv' writes geometry as WKT, losing the CRS. Use the GDAL driver for real formats.
  • COPY ... (FORMAT GDAL, DRIVER 'Parquet') fails — GDAL's Parquet driver is not in the bundled build; use the native FORMAT PARQUET.
  • The native Parquet writer preserves the CRS in the geo metadata key, and GeoPandas reads it back.
  • Glob patterns with ** need read_parquet's recursive form or a shell that expands them.
  • Sort before writing if you know the filter column; it is measurably cheaper afterwards.

FAQ

How do I read a shapefile in DuckDB?

select * from st_read('file.shp') after loading the spatial extension. It goes through GDAL, so GeoPackages, GeoJSON, KML and the rest work the same way.

Should I use st_read or read_parquet for GeoParquet?

read_parquet. It is DuckDB's native columnar reader and can skip columns and row groups; measured, that made a one-column aggregate 43 times faster than reading a row-oriented equivalent.

How do I read a whole folder of files?

Use a glob: st_read('counties/*.shp') or read_parquet('data/**/*.parquet'). The folder becomes one table, and only the aggregate crosses into Python.

Does DuckDB keep the CRS?

ST_Read attaches it to the column type — GEOMETRY('EPSG:4326') — but it is a label, not a guarantee. Joining mismatched CRSs returns zero rows silently, and ST_Transform drops it.

How do I write results back to a GeoPackage?

copy tbl to 'out.gpkg' with (format gdal, driver 'GPKG'). A plain COPY to CSV writes the geometry as WKT text and loses the CRS.

Is converting my shapefiles to GeoParquet worth it?

If the data is queried more than a couple of times, yes. In measurement, 1.79 GB of row-oriented text became 446 MB of Parquet and a repeated aggregation went from 1.01 s to 0.06 s.