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

How to Read Shapefiles, GeoJSON and GeoParquet in DuckDB

On this page

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
print(con.execute("select count(*) from st_read('provinces.shp')").fetchone())
print(con.execute("select count(*) from st_read('data.gpkg', layer='roads')").fetchone())

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

# plain tabular with coordinate columns
iceland = con.execute("""
    select st_point(lon, lat) as geom, name
    from read_csv('points.csv')
""").df()
print(len(iceland), "Icelandic points")
(4596,)
(56600,)
(13472195,)
16296 Icelandic points

Measured on Natural Earth: ST_Read on 7,342 populated places took 0.073 s against GeoPandas' 0.178 s; on 4,596 provinces, 0.052 s against 0.143 s (measured here).

Choosing between ST_Read and read_parquet for a spatial file. two readers, different jobswhat kind of file is it?Parquetread_parquetcolumn + row-group pruninganything GDAL opensst_readuniversal, per-feature costCSV with lon/latread_csv + st_pointno geometry neededUse ST_Read to ingest, read_parquet to work spatialworkflow.io
ST_Read gives universal format support; read_parquet gives pruning.

ST_Read for GDAL formats, read_parquet for Parquet

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('/vsizip/archive.zip/inner/file.shp');  -- GDAL's zip path: the prefix is required

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,472,195 rows (here):

                                     TSV       Parquet
count(distinct country)             0.56 s      0.01 s
group by floor(lon), floor(lat)     0.42 s      0.07 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

read_parquet and read_csv accept wildcards, which replaces the usual "loop over the folder and concatenate" pattern. ST_Read does not: it hands the path to GDAL, and st_read('counties/*.shp') fails with IO Error: Could not open GDAL dataset at: counties/*.shp. For a folder of shapefiles, expand the pattern in Python and union the reads (Example 3):

select count(*) from read_parquet('data/*.parquet');
select count(*) from read_parquet('data/**/*.parquet');

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.

Reading Natural Earth files with ST_Read and with GeoPandas. reading the same shapefilesplaces · ST_Read0.073 splaces · GeoPandas0.178 sprovinces · ST_Read0.052 sprovinces · GeoPandas0.143 s00.3 s7,342 populated places and 4,596 provinces, from the same files on the same machine spatialworkflow.io
Both go through GDAL; the difference is what happens to the rows afterwards.

Inspecting, reading and converting GeoNames files

The data used on this page

The points are GeoNames' full allCountries dump — every feature, about 13.5 million, rebuilt daily at download.geonames.org, so counts drift a little between runs — kept both as the 1.79 GB text file and as places.parquet (also linked as geonames.parquet) with name, latitude, longitude, feature class, country and population. The polygons are Natural Earth's 1:10m admin-1 states and provinces, 4,596 of them, as provinces.shp.

import os
import zipfile

import duckdb
import requests


def fetch(url, path):
    if not os.path.exists(path):
        with requests.get(url, stream=True, timeout=600) as reply:
            reply.raise_for_status()
            with open(path, "wb") as out:
                for chunk in reply.iter_content(1 << 20):
                    out.write(chunk)


fetch("https://download.geonames.org/export/dump/allCountries.zip", "allCountries.zip")   # about 420 MB
if not os.path.exists("allCountries.txt"):
    zipfile.ZipFile("allCountries.zip").extract("allCountries.txt")
if not os.path.exists("places.parquet"):
    duckdb.sql("""
        copy (select column01 as name, column04 as lat, column05 as lon, column06 as feature_class,
                     column08 as country, column14 as population
              from read_csv('allCountries.txt', delim = '\\t', header = false, quote = '',
                            all_varchar = false, sample_size = -1))
        to 'places.parquet' (format parquet)""")
if not os.path.exists("geonames.parquet"):
    os.link("places.parquet", "geonames.parquet")
fetch("https://naciscdn.org/naturalearth/10m/cultural/ne_10m_admin_1_states_provinces.zip", "provinces.zip")
if not os.path.exists("provinces.shp"):
    with zipfile.ZipFile("provinces.zip") as z:
        for member in z.namelist():
            if member.startswith("ne_10m_admin_1_states_provinces."):
                open("provinces" + os.path.splitext(member)[1], "wb").write(z.read(member))
n = duckdb.sql("select count(*) from read_parquet('places.parquet')").fetchone()[0]
print(f"{n:,} GeoNames points; places.parquet {os.path.getsize('places.parquet') / 1e6:.0f} MB, "
      f"allCountries.txt {os.path.getsize('allCountries.txt') / 1e9:.2f} GB")
13,472,195 GeoNames points; places.parquet 306 MB, allCountries.txt 1.79 GB

The examples also need the kinds of file the page talks about, so the setup makes them from real data: Natural Earth's populated places as a shapefile, its 1:10m roads as the roads layer of data.gpkg, Iceland's GeoNames points as points.csv, a counties/ folder with one shapefile of provinces per Nordic country, and a Hive-partitioned by_country/ copy of two countries' points.

import geopandas as gpd

fetch("https://naciscdn.org/naturalearth/10m/cultural/ne_10m_populated_places.zip", "populated_places.zip")
fetch("https://naciscdn.org/naturalearth/10m/cultural/ne_10m_roads.zip", "roads.zip")
if not os.path.exists("populated_places.shp"):
    gpd.read_file("populated_places.zip")[["NAME", "ADM0_A3", "POP_MAX", "geometry"]].to_file("populated_places.shp")
if not os.path.exists("data.gpkg"):
    gpd.read_file("roads.zip")[["type", "name", "continent", "geometry"]].to_file("data.gpkg", layer="roads")
if not os.path.exists("points.csv"):
    duckdb.sql("copy (select name, lat, lon from read_parquet('places.parquet') where country = 'IS') to 'points.csv'")
if not os.path.exists("counties"):
    os.makedirs("counties")
    provinces = gpd.read_file("provinces.shp")
    for code in ("ISL", "NOR", "SWE", "FIN", "DNK"):
        provinces[provinces.adm0_a3 == code][["name", "adm0_a3", "geometry"]].to_file(f"counties/{code}.shp")
if not os.path.exists("by_country"):
    duckdb.sql("copy (select country, name, lat, lon from read_parquet('places.parquet') where country in ('IS', 'MT')) "
               "to 'by_country' (format parquet, partition_by (country))")
print(sorted(os.listdir("counties"))[:6], len(os.listdir("by_country")), "partitions")
['DNK.cpg', 'DNK.dbf', 'DNK.prj', 'DNK.shp', 'DNK.shx', 'FIN.cpg'] 2 partitions

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


meta = describe_spatial_file("data.gpkg")
layer 'roads'
  features   56600
  geometry   Multi Line String
  crs        WGS 84 (EPSG:4326)
  attributes 3

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 and "*" not in path:              # GDAL does not expand globs
        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)


for path, extra in (("counties/NOR.shp", {"columns": "name, adm0_a3"}), ("points.csv", {"where": "name like 'Reykjav%'"}),
                    ("places.parquet", {"columns": "count(*)", "where": "country = 'IS'"}),
                    ("data.gpkg", {"layer": "roads", "columns": "count(*)"})):
    print(path, read_spatial(con, path, **extra).fetchone())
counties/NOR.shp ('Troms', 'NOR')
points.csv ('Reykjavellir', 65.48818, -19.38765)
places.parquet (16296,)
data.gpkg (56600,)

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 glob
    import os
    import time

    started = time.perf_counter()
    files = sorted(glob.glob(pattern))                 # ST_Read cannot glob; union the files
    union = " union all by name ".join(f"select * from st_read('{f}')" for f in files)
    order = f"order by {sort_by}" if sort_by else ""
    con.execute(f"""
        copy (select * exclude geom, st_aswkb(geom) as geometry
              from ({union}) {order})
        to '{target}' (format parquet, compression {compression})
    """)
    rows = con.execute(f"select count(*) from read_parquet('{target}')").fetchone()[0]
    print(f"{len(files)} files, {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")


shapefiles_to_geoparquet(con, "counties/*.shp", "nordic.parquet", sort_by="adm0_a3")
5 files, 74 features → nordic.parquet (0.3 MB) in 0.0s
  sorted by adm0_a3 — 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.26 MB from a shuffled one, and the sorted file was 40% smaller (here).

Example 4 — the CRS and write-back claims, checked

import shutil

for ext in ("shp", "shx", "dbf"):
    shutil.copy(f"counties/ISL.{ext}", f"no_prj.{ext}")
for label, path in (("with .prj", "counties/ISL.shp"), ("without .prj", "no_prj.shp")):
    print(f"{label:13}", con.execute(f"select typeof(geom) from st_read('{path}') limit 1").fetchone()[0])
print("after ST_Transform:", con.execute(
    "select typeof(st_transform(geom, 'EPSG:4326', 'EPSG:3857', always_xy := true)) "
    "from st_read('counties/ISL.shp') limit 1").fetchone()[0])

con.execute("create or replace table isl as select name, geom from st_read('counties/ISL.shp')")
con.execute("copy isl to 'isl.csv'")
print("CSV geometry column starts:", open("isl.csv").read().splitlines()[1].split(",", 1)[1][:40])
con.execute("copy isl to 'isl.gpkg' with (format gdal, driver 'GPKG')")
print("GeoPackage read back:", con.execute("select count(*), typeof(any_value(geom)) from st_read('isl.gpkg')").fetchone())

for path in ("provinces.zip/ne_10m_admin_1_states_provinces.shp", "/vsizip/provinces.zip/ne_10m_admin_1_states_provinces.shp"):
    try:
        print(f"{path[:34]:34}", con.execute(f"select count(*) from st_read('{path}')").fetchone()[0], "features")
    except duckdb.Error as error:
        print(f"{path[:34]:34}", str(error).splitlines()[0][:60])
with .prj     GEOMETRY('EPSG:4326')
without .prj  GEOMETRY
after ST_Transform: GEOMETRY
CSV geometry column starts: "POLYGON ((-14.614613410999937 65.986273
GeoPackage read back: (9, "GEOMETRY('EPSG:4326')")
provinces.zip/ne_10m_admin_1_state IO Error: Could not open GDAL dataset at: provinces.zip/ne_1
/vsizip/provinces.zip/ne_10m_admin 4596 features

Two readers, one through GDAL

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 0.42 s to 0.07 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.

What each DuckDB export writes for a geometry column. four exports, four outcomesgeometry becomesCRS keptformat csvWKT text✕noformat jsonWKT text✕noformat parquetWKB + geo metadata✓yesgdal + GPKGnative geometry✓yesformat gdal, driver ‘Parquet’ does not exist in the bundled build — use the native writer spatialworkflow.io
A plain COPY to CSV succeeds and produces a file no GIS tool recognises.

Missing sidecars and exports that lose geometry

  • 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.

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 0.42 s to 0.07 s.