Contents

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

DuckDBSpatial extends DuckDBDataFrame with 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 DuckDBDataFrame uses to serve GEOMETRY columns, and it underpins the MultiAssaySpatialExperiment on-disk format in 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.

1.1 Installation

if (!require("BiocManager"))
    install.packages("BiocManager")
BiocManager::install("DuckDBSpatial")
library(DuckDBSpatial)
library(sf)

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

2 Lazy spatial columns

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 (WKB). 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 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 sf generics to the geometry column without materializing the table.

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
#> DuckDBColumn of length 6
#>               1               2               3               4               5 
#>      LINESTRING      LINESTRING            <NA> MULTILINESTRING MULTILINESTRING 
#>               6 
#> MULTILINESTRING 
#> 8 Levels: POINT LINESTRING POLYGON MULTIPOINT MULTILINESTRING ... UNKNOWN
head(st_area(geom))           # planar area in the geometry's CRS units
#> DuckDBColumn of length 6
#>  1  2  3  4  5  6 
#>  0  0 NA  0  0  0

sf’s st_geometry_type() reports each row’s OGC geometry type, and st_area() 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:

centroids <- st_centroid(geom)  # each geometry's geometric center point
class(centroids)
#> [1] "DuckDBColumn"
#> attr(,"package")
#> [1] "DuckDBDataFrame"
head(st_as_text(centroids))     # WKT: a human-readable text form of a geometry
#> DuckDBColumn of length 6
#>                                             1 
#> POINT (22.639320225002102 27.917960675006306) 
#>                                             2 
#>                                   POINT EMPTY 
#>                                             3 
#>                                          <NA> 
#>                                             4 
#> POINT (22.639320225002102 27.917960675006306) 
#>                                             5 
#>  POINT (25.75049408851471 24.624752955742647) 
#>                                             6 
#>                                   POINT EMPTY

3 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() constructs a single point from its coordinates, and st_sfc() wraps one or more geometries into the sfc list-column type that sf generics expect, here holding just that one point:

query_pt <- st_sfc(st_point(c(30, 10)))
hits <- st_intersects(geom, query_pt)
head(as.vector(hits))
#>     1     2     3     4     5     6 
#>  TRUE FALSE    NA  TRUE  TRUE FALSE

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:

filtered <- st_filter(df, query_pt)
nrow(filtered)
#> [1] 10

4 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 (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() parses WKT text into an sfc geometry:

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
#> [1]  TRUE  TRUE FALSE
layerSubsetByGeometry(pts, poly, coords = c("x", "y"))  # row indices inside poly
#> [1] 1 2
unlink(pts_path)

5 GeoParquet I/O

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

ddb <- readGeoParquet(spatial_path)
nrow(ddb)
#> [1] 24

st_sf() 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:

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)
#> DuckDBDataFrame with 2 rows and 2 columns
#>          id     geometry
#>   <integer>   <geometry>
#> 1         1 01,01,00,...
#> 2         2 01,01,00,...
unlink(path)

6 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 sf object remains preferable for small layers and for interactive geometry editing.

Within BiocDuckDB, 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.

7 Session information

sessionInfo()
#> R version 4.6.1 Patched (2026-06-24 r90190)
#> Platform: x86_64-apple-darwin20
#> Running under: macOS Ventura 13.7.8
#> 
#> Matrix products: default
#> BLAS:   /Library/Frameworks/R.framework/Versions/4.6-x86_64/Resources/lib/libRblas.0.dylib 
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.6-x86_64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
#> 
#> locale:
#> [1] C/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
#> 
#> time zone: America/New_York
#> tzcode source: internal
#> 
#> attached base packages:
#> [1] stats4    stats     graphics  grDevices utils     datasets  methods  
#> [8] base     
#> 
#> other attached packages:
#> [1] DuckDBSpatial_0.99.6    sf_1.1-2                DuckDBDataFrame_0.99.23
#> [4] IRanges_2.47.2          S4Vectors_0.51.6        BiocGenerics_0.59.12   
#> [7] generics_0.1.4          bit64_4.8.4             BiocStyle_2.41.0       
#> 
#> loaded via a namespace (and not attached):
#>  [1] sass_0.4.10           class_7.3-24          SparseArray_1.13.2   
#>  [4] KernSmooth_2.23-27    lattice_0.23-1        digest_0.6.39        
#>  [7] magrittr_2.0.5        evaluate_1.0.5        grid_4.6.1           
#> [10] bookdown_0.47         blob_1.3.0            fastmap_1.2.0        
#> [13] jsonlite_2.0.0        Matrix_1.7-6          e1071_1.7-17         
#> [16] DBI_1.3.0             BiocManager_1.30.27   purrr_1.2.2          
#> [19] jquerylib_0.1.4       abind_1.4-8           duckdb_1.5.5         
#> [22] cli_3.6.6             rlang_1.3.0           units_1.0-1          
#> [25] dbplyr_2.6.0          XVector_0.53.0        withr_3.0.3          
#> [28] cachem_1.1.0          DelayedArray_0.39.6   yaml_2.3.12          
#> [31] otel_0.2.0            S4Arrays_1.13.0       tools_4.6.1          
#> [34] dplyr_1.2.1           assertthat_0.2.1      vctrs_0.7.3          
#> [37] R6_2.6.1              proxy_0.4-29          matrixStats_1.5.0    
#> [40] lifecycle_1.0.5       classInt_0.4-11       bit_4.6.0            
#> [43] arrow_25.0.0          pkgconfig_2.0.3       bslib_0.12.0         
#> [46] pillar_1.11.1         Rcpp_1.1.2            glue_1.8.1           
#> [49] nanoparquet_0.5.1     xfun_0.60             tibble_3.3.1         
#> [52] tidyselect_1.2.1      MatrixGenerics_1.25.0 knitr_1.51           
#> [55] htmltools_0.5.9       rmarkdown_2.31        wk_0.9.5             
#> [58] compiler_4.6.1