Fixing a Geometry Column Lost on Export from DuckDB

Problem statement

The query is correct, the geometry is right in DuckDB, and the exported file is not spatial:

copy places to 'places.csv' (format csv, header true);

NAME,geom
Vatican City,POINT (12.4533865 41.9032822)
San Marino,POINT (12.4417702 43.9360958)

The geometry survived — as text. It has no CRS, no binary encoding, and nothing downstream will recognise the file as spatial. The same happens with JSON:

{"NAME":"Vatican City","geom":"POINT (12.4533865 41.9032822)"}

And there is a second, sharper failure: the obvious way to write a spatial format does not work at all.

copy places to 'places.parquet' with (format gdal, driver 'Parquet');
-- Binder Error: Could not find GDAL driver: Parquet

Three exports, three different outcomes, and only one of them produces a file another GIS tool will open correctly. A fourth failure is quieter than all three: export a geometry you have transformed, and the file comes back labelled with the wrong coordinate system.

Quick answer

Use the native Parquet writer for Parquet, and the GDAL driver for everything else:

-- GeoParquet: native writer, keeps the CRS in the file's `geo` metadata
copy places to 'places.parquet' (format parquet, compression zstd);

-- GeoPackage, GeoJSON, shapefile: the GDAL driver, with the SRS stated
copy places to 'places.gpkg'    with (format gdal, driver 'GPKG');
copy places to 'places.geojson' with (format gdal, driver 'GeoJSON',
                                      srs 'EPSG:4326');

Verified on Natural Earth's 243 populated places: the native Parquet export round-tripped through GeoPandas with its CRS intact, and the GPKG, GeoJSON and FlatGeobuf exports came back as 243 features in EPSG:4326 (Example 2).

That holds only while the column's type still carries the CRS. After ST_Transform it does not, and the writers fill the gap with a wrong answer rather than an error — so restore it before writing:

copy (select name, st_setcrs(st_transform(geom, 'EPSG:4326', 'EPSG:3857', always_xy := true),
                             'EPSG:3857') as geom
      from places)
to 'places_3857.parquet' (format parquet);

For CSV, convert the geometry deliberately so the loss is a decision rather than an accident:

copy (select name, st_x(geom) as lon, st_y(geom) as lat from places)
to 'places.csv' (format csv, header true);

What each export writes, and how to keep the geometry

1. Understand what each export actually writes

The last row is the trap: the GDAL-driver syntax is the general answer for spatial formats, and Parquet is the one format where it does not apply because the bundled GDAL build has no Parquet driver.

2. Use the native Parquet writer, and check the metadata

copy places to 'places.parquet' (format parquet, compression zstd);

DuckDB writes the geometry as WKB and records the CRS in the file's geo key-value metadata, which is the GeoParquet convention:

import pyarrow.parquet as pq
import json

con.execute("copy places to 'places.parquet' (format parquet, compression zstd)")
meta = pq.read_metadata("places.parquet").metadata
print([k.decode() for k in meta])
geo = json.loads(meta[b"geo"])
print(geo["version"], geo["columns"]["geom"]["crs"]["id"])

['geo']
1.0.0 {'authority': 'EPSG', 'code': 4326}

Measured, GeoPandas read that file back with the CRS present as full PROJJSON.

The writer records what the column type holds, and after a transform it holds nothing, so the crs key is simply left out. GeoParquet defines a missing crs as OGC:CRS84 — longitude and latitude. Web Mercator metres written that way came back from GeoPandas as OGC:CRS84, a claim that coordinates in the millions are degrees. ST_SetCRS, or a cast to GEOMETRY('EPSG:3857'), puts the CRS back on the type, and the file then records it.

3. State the SRS when the driver needs it

The GDAL drivers take the CRS from the column type when it has one. When it does not — after a transform — a GeoPackage is written with an "Undefined geographic SRS" whose axes are latitude and longitude. Pass srs for those cases, and remember what it does: it labels the output, it does not transform it.

copy places to 'out.geojson' with (format gdal, driver 'GeoJSON', srs 'EPSG:4326');

-- transform first, then label with the CRS the coordinates are now in
copy (select name, st_transform(geom, 'EPSG:4326', 'EPSG:27700', always_xy := true) as geom
      from places)
to 'out.shp' with (format gdal, driver 'ESRI Shapefile', srs 'EPSG:27700');

An earlier version of this page wrote the untransformed places with srs 'EPSG:27700'. The shapefile came back as EPSG:27700 with coordinates from −175.2 to 179.2 — degrees labelled as metres, which puts all 243 towns within a few hundred metres of one another in the Atlantic, south-west of the Isles of Scilly.

4. Check the driver exists before relying on it

select short_name, long_name from st_drivers() where can_create order by 1;

AmigoCloud, CSV, Carto, DGN, DXF, ESRI Shapefile, Elasticsearch, FlatGeobuf, GML, GPKG,
GPSBabel, GPX, GeoJSON, GeoJSONSeq, GeoRSS, Geoconcept, JML, JSONFG, KML, MVT,
MapInfo File, MapML, Memory, NGW, ODS, OGR_GMT, OpenFileGDB, PGDUMP, PMTiles, S57,
SQLite, Selafin, VDV, WAsP, XLSX

Those are the 35 in DuckDB 1.5.5's spatial extension. The bundled GDAL is a smaller build than a system GDAL, and Parquet and Arrow are the notable absences. If the driver you want is missing, write Parquet or GeoPackage from DuckDB and convert with another tool.

5. Make the CSV loss explicit

There are legitimate reasons to write CSV — a spreadsheet, a system that only accepts it. Do it deliberately:

-- coordinates as numbers: readable everywhere, no CRS
copy (select name, round(st_x(geom), 6) as lon, round(st_y(geom), 6) as lat,
             'EPSG:4326' as crs
      from places)
to 'places.csv' (format csv, header true);

-- or WKT, explicitly, with the CRS in a column
copy (select name, st_astext(geom) as wkt, 'EPSG:4326' as crs from places)
to 'places.csv' (format csv, header true);

Six decimal places is about 11 cm at the equator, which is beyond the accuracy of nearly all source data.

6. Verify the export by reading it back

import geopandas as gpd

gdf = gpd.read_parquet("places.parquet")
assert gdf.crs is not None, "the CRS did not survive the export"
assert gdf.geometry.notna().all(), "geometry is missing"
print(f"{len(gdf):,} features, CRS {gdf.crs.to_string()}")

243 features, CRS EPSG:4326

Reading the file in the tool that will consume it is the only check that matters, and it takes three lines.

An export that cannot silently lose geometry

The data used on this page

places is Natural Earth's 1:110m populated places — 243 points in EPSG:4326 — read into DuckDB with ST_Read.

import io
import os
import zipfile

import duckdb
import requests

if not os.path.exists("ne_110m_populated_places.shp"):
    reply = requests.get("https://naciscdn.org/naturalearth/110m/cultural/ne_110m_populated_places.zip",
                         timeout=120, headers={"User-Agent": "spatialworkflow.io guide"})
    zipfile.ZipFile(io.BytesIO(reply.content)).extractall(".")

con = duckdb.connect()
con.install_extension("spatial")
con.load_extension("spatial")
con.execute("create or replace table places as "
            "select NAME as name, geom from st_read('ne_110m_populated_places.shp')")
print(con.execute("select count(*), any_value(typeof(geom)) from places").fetchone())

(243, "GEOMETRY('EPSG:4326')")

Example 1 — an export function that cannot silently lose geometry

import duckdb


SPATIAL_FORMATS = {
    ".parquet": ("native", "format parquet, compression zstd"),
    ".gpkg":    ("gdal", "GPKG"),
    ".geojson": ("gdal", "GeoJSON"),
    ".json":    ("gdal", "GeoJSON"),
    ".shp":     ("gdal", "ESRI Shapefile"),
    ".fgb":     ("gdal", "FlatGeobuf"),
}


def export_spatial(con, relation, path, crs=None, overwrite=True, geom="geom"):
    """Write a spatial format, or refuse and say why.

    crs: required when the geometry column no longer carries one (after
    ST_Transform); refused if it contradicts the one the column declares.
    """
    import os

    declared = con.execute(f"select st_crs({geom}) from {relation} limit 1").fetchone()[0]
    if declared and crs and declared != crs:
        raise ValueError(f"the column says {declared}, you said {crs}: srs labels, "
                         f"it does not transform — use ST_Transform first")
    crs = crs or declared
    if crs is None:
        raise ValueError("the geometry carries no CRS (it was transformed or constructed); "
                         "pass crs= so the file does not default to lon/lat")
    source = f"(select * replace (st_setcrs({geom}, '{crs}') as {geom}) from {relation})"

    ext = os.path.splitext(path)[1].lower()
    if ext in (".csv", ".tsv"):
        raise ValueError(
            f"{ext} cannot hold geometry. Either write coordinates explicitly:\n"
            f"  copy (select *, st_x(geom) lon, st_y(geom) lat from {relation}) "
            f"to '{path}' (format csv)\n"
            f"or choose a spatial format: {', '.join(sorted(SPATIAL_FORMATS))}")

    if ext not in SPATIAL_FORMATS:
        raise ValueError(f"no known writer for {ext!r}")

    kind, spec = SPATIAL_FORMATS[ext]
    if overwrite and os.path.exists(path):
        os.remove(path)

    if kind == "native":
        con.execute(f"copy {source} to '{path}' ({spec})")
    else:
        available = {row[0] for row in con.execute(
            "select short_name from st_drivers() where can_create").fetchall()}
        if spec not in available:
            raise ValueError(
                f"the bundled GDAL has no '{spec}' driver. "
                f"Write .parquet or .gpkg and convert with another tool.")
        con.execute(f"copy {source} to '{path}' "
                    f"with (format gdal, driver '{spec}', srs '{crs}')")

    print(f"wrote {path} ({os.path.getsize(path) / 1e3:,.1f} kB) as {kind}/{spec}, {crs}")
    return path


for path in ("places.parquet", "places.gpkg", "places.geojson", "places.fgb"):
    export_spatial(con, "places", path)

con.execute("create or replace table places_3857 as select name, "
            "st_transform(geom, 'EPSG:4326', 'EPSG:3857', always_xy := true) as geom from places")
try:
    export_spatial(con, "places_3857", "places_3857.parquet")
except ValueError as error:
    print("refused:", error)
for path in ("places_3857.parquet", "places_3857.gpkg"):
    export_spatial(con, "places_3857", path, crs="EPSG:3857")
try:
    export_spatial(con, "places", "places.csv")
except ValueError as error:
    print("refused:", str(error).splitlines()[0])

wrote places.parquet (7.3 kB) as native/format parquet, compression zstd, EPSG:4326
wrote places.gpkg (122.9 kB) as gdal/GPKG, EPSG:4326
wrote places.geojson (34.2 kB) as gdal/GeoJSON, EPSG:4326
wrote places.fgb (35.5 kB) as gdal/FlatGeobuf, EPSG:4326
refused: the geometry carries no CRS (it was transformed or constructed); pass crs= so the file does not default to lon/lat
wrote places_3857.parquet (8.2 kB) as native/format parquet, compression zstd, EPSG:3857
wrote places_3857.gpkg (122.9 kB) as gdal/GPKG, EPSG:3857
refused: .csv cannot hold geometry. Either write coordinates explicitly:

Example 2 — verifying the round trip

def verify_export(path, expected_rows=None, expected_crs=None):
    """Read the file back the way its consumer will."""
    import geopandas as gpd
    import os

    reader = gpd.read_parquet if path.endswith(".parquet") else gpd.read_file
    gdf = reader(path)

    problems = []
    if gdf.geometry.isna().any():
        problems.append(f"{gdf.geometry.isna().sum()} null geometries")
    if gdf.crs is None:
        problems.append("no CRS — the export lost it")
    elif expected_crs and gdf.crs.to_string() != expected_crs:
        found = gdf.crs.to_string()
        found = found if len(found) < 40 else gdf.crs.name      # WKT when there is no code
        problems.append(f"CRS is {found!r}, expected {expected_crs}")
    if expected_rows and len(gdf) != expected_rows:
        problems.append(f"{len(gdf):,} rows, expected {expected_rows:,}")

    for problem in problems:
        print(f"  ! {os.path.basename(path)}: {problem}")
    if not problems:
        print(f"  {os.path.basename(path)}: {len(gdf):,} features, "
              f"CRS {gdf.crs.to_string()}, geometry intact")
    return not problems


for path in ("places.parquet", "places.gpkg", "places.geojson", "places.fgb"):
    verify_export(path, expected_rows=243, expected_crs="EPSG:4326")
for path in ("places_3857.parquet", "places_3857.gpkg"):
    verify_export(path, expected_rows=243, expected_crs="EPSG:3857")

# the same transformed table written without restoring the CRS
con.execute("copy places_3857 to 'unlabelled.parquet' (format parquet)")
con.execute("copy places_3857 to 'unlabelled.gpkg' with (format gdal, driver 'GPKG')")
for path in ("unlabelled.parquet", "unlabelled.gpkg"):
    verify_export(path, expected_rows=243, expected_crs="EPSG:3857")

  places.parquet: 243 features, CRS EPSG:4326, geometry intact
  places.gpkg: 243 features, CRS EPSG:4326, geometry intact
  places.geojson: 243 features, CRS EPSG:4326, geometry intact
  places.fgb: 243 features, CRS EPSG:4326, geometry intact
  places_3857.parquet: 243 features, CRS EPSG:3857, geometry intact
  places_3857.gpkg: 243 features, CRS EPSG:3857, geometry intact
  ! unlabelled.parquet: CRS is 'OGC:CRS84', expected EPSG:3857
  ! unlabelled.gpkg: CRS is 'Undefined geographic SRS', expected EPSG:3857

Example 3 — turning a WKT CSV back into geometry

def recover_from_wkt_csv(con, csv_path, crs="EPSG:4326", wkt_column="geom"):
    """Rescue a CSV that was exported with geometry as text."""
    con.execute("load spatial")
    return con.execute(f"""
        select * exclude {wkt_column},
               st_geomfromtext({wkt_column}) as geom
        from read_csv('{csv_path}')
    """).df(), crs


con.execute("copy places to 'places.csv' (format csv, header true)")
frame, crs = recover_from_wkt_csv(con, "places.csv")
print(len(frame), "rows;", frame.dtypes["geom"], "— CRS supplied from outside:", crs)

243 rows; object — CRS supplied from outside: EPSG:4326

The CRS has to be supplied from outside, because the CSV never carried it. That is the whole cost of the accidental export: a value that has to be remembered rather than read.

Why CSV writes text and the CRS slips away

Why CSV writes WKT rather than failing

CSV has one type: text. Faced with a geometry column, DuckDB does the most useful thing available and writes its textual representation, which is at least reversible.

The alternative — refusing — would break the legitimate case of exporting WKT deliberately. So the behaviour is defensible and it is silent, which is why an export that "worked" produces a file no GIS tool recognises.

Why the Parquet GDAL driver does not exist here

GDAL has a Parquet driver, but it depends on Arrow and is not compiled into the static GDAL bundled with the DuckDB spatial extension. Hence Binder Error: Could not find GDAL driver: Parquet.

This is not a limitation in practice, because DuckDB's native Parquet writer is better for the job: it writes proper GeoParquet metadata, and it is the same writer used for every other Parquet export.

The trap is only that the GDAL syntax is the general answer for spatial formats and Parquet is the exception.

Why the CRS is so easily lost

DuckDB carries the CRS on the geometry column's type, and the annotation does not survive most operations — ST_Transform returns a bare GEOMETRY, and so does any constructed geometry.

Truncated names, GeoJSON's CRS and mixed types

• Shapefile truncates field names to 10 characters and has a 2 GB limit.
• GeoJSON should be EPSG:4326. RFC 7946 allows another CRS only "by prior arrangement", and GDAL will write one; most readers ignore it.
• Mixed geometry types are fine in DuckDB and GeoPackage, and rejected by shapefile.

Internal links

• How to read shapefiles, GeoJSON and GeoParquet in DuckDB — the reading side
• CRS in DuckDB: why ST_Transform moves your data to the wrong place — why the CRS goes missing
• How to move data between DuckDB and GeoPandas — the in-process handover
• GeoParquet and columnar storage explained — what the native writer produces
• WKT, WKB and GeoJSON explained — the representations
• Fixing ST_Read failures on spatial files in DuckDB — the reader's failures
• How to clean CSV coordinates into a GeoDataFrame — recovering a text export
• How to export GeoJSON from GeoPandas — the other library's exports

FAQ

Why is my exported CSV's geometry just text?

Because CSV has no geometry type. DuckDB writes WKT, which is reversible but carries no CRS and is not recognised as spatial by any GIS tool.

How do I write GeoParquet from DuckDB?

copy tbl to 'x.parquet' (format parquet) — the native writer. It stores WKB plus the geo metadata key, and GeoPandas reads the CRS back.

Why does format gdal, driver 'Parquet' fail?

The bundled GDAL build has no Parquet driver. Use the native Parquet writer instead, which writes proper GeoParquet metadata anyway.

How do I write a GeoPackage?

copy tbl to 'x.gpkg' with (format gdal, driver 'GPKG'). Verified: it round-trips through GeoPandas with the geometry and EPSG:4326 intact.

Why did my export lose the CRS?

DuckDB tracks CRS on the column type and drops it in most operations, including ST_Transform, and the writers then fall back to lon/lat rather than failing. Restore it with ST_SetCRS(geom, 'EPSG:…') before writing — that fixes the Parquet and GDAL writers alike.

How do I check the export worked?

Read it back with the library that will consume it and assert the row count, the CRS and that no geometry is null. Three lines, at the point of writing.