DuckDB Spatial Analytical SQL for GIS

Writing Python UDFs Over Geometry Columns

A geometry user-defined function is a small amount of code with a large amount of surrounding cost, and most of that cost is decided by four declarations made at registration time rather than by the function body — this walkthrough, part of the Python spatial UDFs and predicate pushdown reference, covers those declarations, the three shapes a geometry UDF can take, and the failures that are specific to passing geometry across the boundary rather than to UDFs in general.

Root-Cause Analysis: what makes a geometry UDF expensive

The function body is rarely the problem. Five other things are, and they compound.

  • The boundary crossing is per row. A scalar registration is invoked once per row, and each invocation acquires the interpreter lock, marshals the argument, and marshals the result. On a million-row table that is a million round trips regardless of how trivial the body is.
  • Geometry has to be serialised twice. The engine holds a native GEOMETRY; Python wants bytes or an object. Passing geometry in means a serialisation, and returning geometry means a deserialisation, so a blob-to-blob function pays the conversion at both ends of every call.
  • The surrounding expression stops being vectorised. An expression containing a UDF is evaluated row at a time in its entirety, so the ST_Area sitting next to it in the same SELECT list stops being a SIMD kernel too.
  • The optimizer cannot see through it. A predicate involving a UDF cannot be pushed into a scan, so rows the function would have rejected are read and decoded first. This is usually a larger cost than the function itself.
  • Its memory is outside the budget. memory_limit governs the engine. Anything the Python side allocates — a captured object, an accumulating list, a library’s internal cache — is on top of it and invisible to it.

The distinguishing question is whether the function is in a projection or in a predicate. In a projection it costs its own per-row overhead. In a predicate it costs that plus every row the query could have avoided reading.

The four registration choices and what each one decides Return type, argument types, null handling and purity — each changes how often the function is called and what it costs per call. REGISTRATION CHOICE WHAT IT DECIDES GETTING IT WRONG COSTS declared return type storage and cast-free use a silent conversion per row declared argument types what the engine hands you an unnecessary conversion per row null handling whether nulls reach the function one invocation per null row purity / side effects whether the optimizer may fold calls repeated evaluation of the same input

Four small declarations, each of which is a per-row cost when it is wrong.

Deterministic Configuration 

import duckdb
from shapely import from_wkb, to_wkb

con = duckdb.connect("gis.duckdb")
con.execute("INSTALL spatial; LOAD spatial")

# The engine's budget does not cover the Python side, so leave headroom for
# whatever the function and its libraries allocate on top.
con.execute("SET memory_limit = '4GB'")   # engine only
con.execute("SET threads = 8")            # does not speed up the UDF itself

Optimized Execution Pattern

Three shapes for a geometry UDF, by boundary cost A scalar blob-to-blob function is simplest and dearest; a scalar blob-to-double saves the return serialisation; a vectorised registration crosses the boundary once per batch. scalar: WKB → WKB simplest to write one crossing per row most expensive shape and it serialises the return value on every row too scalar: WKB → DOUBLE a measurement rather than a new geometry same crossing, cheaper return no serialisation on the way back out vectorised: array → array one crossing per batch rather than per row removes the largest cost still opaque to the optimizer, so still no pushdown

The blob-to-double shape is worth reaching for whenever the answer is a measurement rather than a new shape, because it halves the serialisation: the geometry comes in as bytes and a plain number goes out. That is a meaningful saving on a function whose body is short, where marshalling dominates.

# ANTI-PATTERN: returns a Shapely object and takes the native type, so the
# engine converts at both ends and the declared types do not match reality.
def bad_centroid_x(geom):
    return from_wkb(geom).centroid          # a Python object, not a value

con.create_function("bad_centroid_x", bad_centroid_x)
# PATTERN: bytes in, a plain scalar out, nulls short-circuited, types declared.
def centroid_x(wkb: bytes) -> float:
    if wkb is None:
        return None
    return from_wkb(wkb).centroid.x

con.create_function(
    "centroid_x", centroid_x,
    [duckdb.typing.BLOB], duckdb.typing.DOUBLE,   # declared, so no per-row cast
    null_handling="special",                       # engine short-circuits nulls
)

con.execute("SELECT parcel_id, centroid_x(ST_AsWKB(geom)) FROM parcels LIMIT 5").fetchall()

When the function should not be a UDF at all

The example above is deliberately one that should never be written, because ST_X(ST_Centroid(geom)) already exists and is vectorised. That is the common case: a UDF is written for something SQL already does, because the author reached for the language they were thinking in rather than the one the data was in. Before registering anything, check the function catalogue — SELECT function_name FROM duckdb_functions() WHERE function_name ILIKE 'st_%' is one query and it regularly ends the discussion.

Where a genuine gap exists, the second question is whether the operation is per-geometry. A function that needs to see several geometries cannot be a scalar UDF at all, and trying to make it one usually means passing a serialised collection in as an argument, which is far more expensive than the join it is standing in for.

Five ways a geometry UDF fails, and how each presents Unparsed blobs error immediately; returning objects or WKT costs conversion; unhandled nulls error partway through; a captured object inflates memory outside the engine budget. FAILURE HOW IT PRESENTS SEVERITY blob not parsed an exception on the first row friendly — found at once returning a Shapely object conversion error, or stringified depends on the declared type returning WKT instead of WKB works; re-parsed every row a performance failure null input unhandled an exception partway through a run wasteful — work already done large object in the closure memory outside the engine budget invisible to memory_limit

Only the first row fails early. The rest fail late, or not at all.

Diagnostic Queries & Plan Validation 

-- The UDF appears as its own node. What matters is how many rows reach it:
-- compare that against the table size to see whether anything filtered first.
EXPLAIN ANALYZE
SELECT count(*) FROM parcels WHERE centroid_x(ST_AsWKB(geom)) > 2.3;

Diagnostic — does the query respond to threads? Run it at two threads and at eight. A vectorised query improves substantially; one dominated by a UDF barely moves, because the interpreter lock serialises the Python half. It is the cheapest available test for a per-row bottleneck and needs nothing installed.

for t in (2, 8):
    con.execute(f"SET threads = {t}")
    # time the same statement; a flat result means the UDF, not the engine, is the cost

Geometry Validation & Fallback Routing

A UDF is a poor place to discover invalid geometry, because the exception arrives partway through a run with an arbitrary amount of work already done and no indication of which row caused it. Validate before the function rather than inside it, and have the function fail safe rather than fail hard.

def safe_metric(wkb: bytes) -> float:
    """Return NULL rather than raising, so one bad row does not lose the run.
    The count of NULLs afterwards is the diagnostic."""
    if wkb is None:
        return None
    try:
        return _expensive_library_call(from_wkb(wkb))
    except Exception:
        return None
-- Then the count of NULLs tells you how many rows the function refused,
-- which is a number you can act on rather than a stack trace you cannot.
SELECT count(*) AS rows_total,
       count(*) FILTER (WHERE metric IS NULL) AS rows_refused
FROM (SELECT safe_metric(ST_AsWKB(geom)) AS metric FROM parcels WHERE ST_IsValid(geom));

Testing a UDF outside the engine

A UDF that is only ever exercised through SQL is difficult to debug, because a failure arrives as a query error with the row that caused it already discarded. Keeping the function a plain Python function, registered separately, means it can be tested directly against a handful of WKB values — including the degenerate ones — before it ever meets a table.

# The function is ordinary Python; registration is a separate line. That
# separation is what makes it testable against the inputs that break it.
from shapely import to_wkb, Point, Polygon

CASES = {
    "point":        to_wkb(Point(2.35, 48.85)),
    "empty":        to_wkb(Polygon()),
    "null":         None,
}
for name, wkb in CASES.items():
    print(name, safe_metric(wkb))       # expect a value or None, never an exception

Running that list on every change costs nothing and catches the two failures that are otherwise found halfway through a production run: the unhandled null, and the empty geometry whose library call raises rather than returning a sensible zero.

Frequently Asked Questions

Should geometry arrive as a blob or as the native type?

As a blob. The Python side cannot do anything with a native GEOMETRY handle, so the engine would convert it anyway — declaring the argument as BLOB and passing ST_AsWKB(geom) makes that conversion explicit and stops a second one being added implicitly.

Why declare the return type?

Because an undeclared return type makes the engine infer one and add a cast on every row when the inference is wrong. Declaring it is one argument at registration and removes a per-row cost that is otherwise invisible — it does not appear as a separate plan node, it is simply time inside the function’s own node.

Does null handling matter for performance?

On a column with many nulls, yes. By default the function is invoked for null inputs and has to return early, which is a full boundary crossing to produce a null. Declaring that the engine should short-circuit nulls skips the invocation entirely, which on a sparse column is a large fraction of the calls.

Can a UDF see more than one row?

Not a scalar one — it receives a single value. A vectorised registration, where the driver offers one, receives a whole array and returns one, which removes most of the boundary cost. What neither offers is arbitrary access to other rows, which is what a join or a window function is for.

How do I stop a UDF exhausting memory?

Keep it pure and keep it small. Anything captured in the closure or accumulated across calls sits outside memory_limit entirely, so the engine will not spill it and will not warn about it. If the function genuinely needs a large lookup structure, pass it in as a registered table and join against it instead of capturing it.

Up: Python Spatial UDFs and Predicate Pushdown