---
title: "Introduction to DuckDBSpatial"
author:
- name: Patrick Aboyoun
  email: aboyounp@gene.com
  affiliation: Genentech, Inc.
date: "Compiled: `r BiocStyle::doc_date()`; Modified: 6 July 2026"
package: DuckDBSpatial
vignette: >
  %\VignetteIndexEntry{1. Introduction to DuckDBSpatial}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
output:
  BiocStyle::html_document:
    number_sections: true
    toc: true
    toc_depth: 2
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>",
                      error = FALSE, warning = FALSE, message = FALSE)

# Probe once whether DuckDB's spatial extension is available in this build
# environment (it is fetched on first use, or loaded from a provisioned cache).
# The live spatial examples below are gated on this so the vignette still builds
# where the extension cannot be obtained (e.g. an offline developer machine);
# CI / a networked build renders them for real.
.spatial_ok <- tryCatch({
    con <- DuckDBDataFrame::acquireDuckDBConn()
    DBI::dbGetQuery(con,
        "SELECT ST_Area(ST_GeomFromText('POLYGON((0 0,0 1,1 1,1 0,0 0))')) AS a")
    TRUE
}, error = function(e) FALSE)
```

```{r spatial-unavailable-note, echo=FALSE, results='asis', eval=!.spatial_ok}
cat("> **Note:** DuckDB's `spatial` extension could not be loaded in this build ",
    "environment, so the live spatial examples below are shown but not evaluated. ",
    "See the README for how to make the extension available (it is fetched ",
    "automatically where the DuckDB extension repository is reachable).\n", sep = "")
```

# Introduction

Spatial datasets are increasingly large: a spatial-transcriptomics experiment can
carry millions of transcript centroids and cell-boundary polygons, and imaging or
geographic layers run to hundreds of millions of features. The
[sf](https://r-spatial.github.io/sf/) package is the standard R interface for such
data, `st_area()`, `st_intersects()`, `st_centroid()`, but an ordinary `sf`
object holds every geometry in memory.

`r Biocpkg("DuckDBSpatial")` extends `r Biocpkg("DuckDBDataFrame")` with
`r CRANpkg("sf")`-compatible spatial methods that run on **DuckDB-backed columns
and tables**, powered by DuckDB's native spatial extension. The geometries stay
**on disk** in columnar Parquet, and spatial operations are recorded as lazy SQL
that DuckDB pushes down, so you can measure, filter, and transform geometries
**without loading the layer into memory**. It reads and writes **GeoParquet 1.0**,
so the same files interoperate with GeoPandas, GDAL, QGIS, and DuckDB itself.

Within the **BiocDuckDB** suite it is the spatial layer: it is what
`r Biocpkg("DuckDBDataFrame")` uses to serve `GEOMETRY` columns, and it underpins
the `MultiAssaySpatialExperiment` on-disk format in `r Biocpkg("BiocDuckDB")`.

This vignette is a practical introduction. For how the `sf` generics map onto
DuckDB spatial SQL and how GeoParquet I/O works, see *Design and extension of
DuckDBSpatial*.

## Installation

```{r install, eval=FALSE}
if (!require("BiocManager"))
    install.packages("BiocManager")
BiocManager::install("DuckDBSpatial")
```

```{r load}
library(DuckDBSpatial)
library(sf)
```

The DuckDB spatial extension loads automatically the first time a spatial
operation runs (it is fetched once and cached).

# Lazy spatial columns

`r CRANpkg("sf")` represents a set of geometries (points, lines, polygons, ...)
as a `geometry` column: a list where each entry is one feature's shape, stored
internally as [Well-Known Binary](https://libgeos.org/specifications/wkb/)
(WKB). `r Biocpkg("DuckDBSpatial")`'s spatial columns hold the same WKB, just
kept on disk in DuckDB rather than materialized in memory, so the familiar `sf`
generics (see the
[sf reference index](https://r-spatial.github.io/sf/reference/index.html) for
the full list) work the same way, they just push the computation down to SQL
instead of running in R.

The package bundles a small example layer as partitioned Parquet. Open it as a
lazy `DuckDBDataFrame` and apply `r CRANpkg("sf")` generics to the geometry column
without materializing the table.

```{r columns, eval=.spatial_ok}
spatial_path <- system.file("extdata", "spatial", package = "DuckDBSpatial")
df <- DuckDBDataFrame(spatial_path)
df <- df[which(!is.na(df$type)), ]

geom <- df[["geometry"]]
head(st_geometry_type(geom))  # one of POINT, LINESTRING, POLYGON, ... per row
head(st_area(geom))           # planar area in the geometry's CRS units
```

`r CRANpkg("sf")`'s
[`st_geometry_type()`](https://r-spatial.github.io/sf/reference/st_geometry_type.html)
reports each row's OGC geometry type, and
[`st_area()`](https://r-spatial.github.io/sf/reference/geos_measures.html)
its planar area (this bundled example has no coordinate reference system set,
so the values are in raw coordinate units rather than a real-world unit like
square meters).

Measurements (`st_area()`, `st_length()`, `st_perimeter()`) and geometry
transforms (`st_centroid()`, `st_buffer()`, `st_convex_hull()`) all return
**lazy `DuckDBColumn` objects**, nothing is computed until the values are
pulled:

```{r transforms, eval=.spatial_ok}
centroids <- st_centroid(geom)  # each geometry's geometric center point
class(centroids)
head(st_as_text(centroids))     # WKT: a human-readable text form of a geometry
```

# Spatial predicates and filtering

Spatial predicates, true/false tests of how two geometries relate, such as
`st_intersects()`, `st_within()`, and `st_contains()`, build SQL against a query
geometry and stay lazy. The query geometry itself is built the same way any
`sf` geometry is: [`st_point()`](https://r-spatial.github.io/sf/reference/st.html)
constructs a single point from its coordinates, and
[`st_sfc()`](https://r-spatial.github.io/sf/reference/sfc.html) wraps one or
more geometries into the `sfc` list-column type that `sf` generics expect,
here holding just that one point:

```{r predicates, eval=.spatial_ok}
query_pt <- st_sfc(st_point(c(30, 10)))
hits <- st_intersects(geom, query_pt)
head(as.vector(hits))
```

The table-level `st_filter()` is the convenient way to keep only the rows whose
geometry satisfies a predicate against a query, evaluated in DuckDB:

```{r filter, eval=.spatial_ok}
filtered <- st_filter(df, query_pt)
nrow(filtered)
```

# Layer-level engines for coordinate columns

Point data often lives as plain `x`/`y` columns (for example transcript
centroids) rather than a geometry column. The `layer*` helpers run spatial
queries directly on coordinate columns, so no geometry column is needed. The
query polygon below is written as
[WKT](https://en.wikipedia.org/wiki/Well-known_text_representation_of_geometry)
(Well-Known Text), a standard textual encoding for geometries: a comma-separated
list of `x y` vertex coordinates, closed by repeating the first vertex last.
[`st_as_sfc()`](https://r-spatial.github.io/sf/reference/st_as_sfc.html) parses
WKT text into an `sfc` geometry:

```{r layers, eval=.spatial_ok}
pts_path <- tempfile(fileext = ".csv")
write.csv(data.frame(x = c(1, 5, 30), y = c(1, 5, 10)), pts_path, row.names = FALSE)
pts <- DuckDBDataFrame(pts_path, datacols = c("x", "y"))

poly <- st_as_sfc("POLYGON((0 0, 6 0, 6 6, 0 6, 0 0))")  # a 6x6 square, from WKT
layerSpatialOverlaps(pts, poly, coords = c("x", "y"))   # which points fall in poly
layerSubsetByGeometry(pts, poly, coords = c("x", "y"))  # row indices inside poly
unlink(pts_path)
```

# GeoParquet I/O

`readGeoParquet()` opens a GeoParquet file as a lazy `DuckDBDataFrame` with a
native `GEOMETRY` column; `writeGeoParquet()` writes an `r CRANpkg("sf")` object
with GeoParquet 1.0 metadata (requires the suggested `r CRANpkg("nanoparquet")`
package). The result is readable by any GeoParquet-aware tool.

```{r geoparquet-read, eval=.spatial_ok}
ddb <- readGeoParquet(spatial_path)
nrow(ddb)
```

[`st_sf()`](https://r-spatial.github.io/sf/reference/sf.html) combines a data
frame of feature attributes (here just an `id` column) with an `sfc` geometry
column into a full `sf` object, the standard way to pair non-spatial attributes
with geometries:

```{r geoparquet-write, eval=.spatial_ok && requireNamespace("nanoparquet", quietly=TRUE)}
pts_sf <- st_sf(id = 1:2,
                geometry = st_sfc(st_point(0:1), st_point(2:3)))
path <- tempfile(fileext = ".parquet")
writeGeoParquet(pts_sf, path)
readGeoParquet(path)
unlink(path)
```

# When to use DuckDBSpatial

A good fit when the spatial layer is **larger than memory**, when the workload is
**filter-heavy** (selecting features by region or predicate before the expensive
step), or when the data lives on disk as **GeoParquet** shared with other spatial
tooling. An in-memory `r CRANpkg("sf")` object remains preferable for small layers
and for interactive geometry editing.

Within `r Biocpkg("BiocDuckDB")`, `r Biocpkg("DuckDBSpatial")` is what makes a
`MultiAssaySpatialExperiment` hold its cell boundaries and landmarks on disk as
lazy `GEOMETRY` columns. For the SQL translation and GeoParquet details, see
*Design and extension of DuckDBSpatial*.

# Session information

```{r sessioninfo}
sessionInfo()
```