r-duckspatial

duckspatial: Memory-Efficient Spatial Analysis

Overview

duckspatial bridges R's sf ecosystem with DuckDB's spatial extension for memory-efficient analysis of large spatial datasets. Data stays in DuckDB until explicitly collected, enabling operations on datasets larger than RAM.

Core principle: Lazy evaluation keeps data outside R memory. Operations execute in DuckDB's query engine with optimization before results return.

References

Read references/API.md before writing code.

• references/API.md - Complete function reference and spatial operations

When to Use

• Dataset too large for sf to load into memory
• Spatial joins with millions of features
• Areal-weighted interpolation across large geographies
• Repeated spatial operations benefit from DuckDB's query optimization
• Working with spatial files larger than 1GB

When NOT to Use

• Small datasets (<100K features) - sf is simpler
• Complex spatial topology operations not in DuckDB spatial extension
• Need GEOS-specific functionality (sf provides more geometry operations)
• Interactive spatial visualization - materialize with collect() first

Quick Reference

Operation	Function	Key Parameters
Load file	ddbs_open_dataset(path)	Returns lazy duckspatial_df
Convert from sf	as_duckspatial_df(sf_obj)	Registers sf as lazy table
Materialize	ddbs_collect() or st_as_sf()	Brings data into R memory
Spatial join	ddbs_join(x, y, join, predicate)	join: left/inner/right/full
Spatial filter	ddbs_filter(x, y, predicate)	Subset by spatial relation
Crop to extent	ddbs_crop(x, y)	Clip geometries to y extent
Areal interpolation	ddbs_interpolate_aw(x, y, values)	Weight by overlap area
Buffer	ddbs_buffer(x, distance)	Distance in CRS units
Union/dissolve	ddbs_union(x)	Merge into single geometry
Transform CRS	ddbs_transform(x, crs)	EPSG code or proj4string
Validate	ddbs_make_valid(x)	Fix invalid geometries

See references/API.md for complete function reference.

Core Pattern

library(duckspatial)

# Load large dataset as lazy table (not in R memory)
buildings <- ddbs_open_dataset("buildings.gpkg")

# Or convert existing sf object
neighborhoods <- as_duckspatial_df(sf_neighborhoods)

# Chain operations - all happen in DuckDB
result <- buildings |>
  ddbs_make_valid() |>                    # Fix geometries
  ddbs_transform(crs = 32633) |>          # Reproject to UTM
  ddbs_join(neighborhoods,
            join = "left",
            predicate = "within") |>       # Spatial join
  ddbs_collect()                           # NOW bring into R memory

# Result is sf object, ready for plotting/further analysis
plot(result["neighborhood_name"])

Common Workflows

Spatial Join (Point-in-Polygon)

# Points to polygons - find which polygon each point falls in
points <- ddbs_open_dataset("locations.geojson")
zones <- as_duckspatial_df(zone_boundaries)

joined <- ddbs_join(points, zones,
                    join = "left",
                    predicate = "within") |>
  ddbs_collect()

Areal-Weighted Interpolation

# Transfer population from census tracts to neighborhoods
census <- ddbs_open_dataset("census_tracts.gpkg")
neighborhoods <- as_duckspatial_df(neighborhoods_sf)

interpolated <- ddbs_interpolate_aw(
  census,
  neighborhoods,
  values = "population"  # Extensive variable (sums)
) |>
  ddbs_collect()

Large File Processing

# Process multi-GB shapefile without loading entirely
parcels <- ddbs_open_dataset("parcels_nationwide.shp")

city_parcels <- parcels |>
  ddbs_filter(study_area, predicate = "intersects") |>
  ddbs_make_valid() |>
  ddbs_simplify(tolerance = 0.5) |>
  ddbs_collect()

Integration with sf and dplyr

# duckspatial_df supports dplyr verbs
filtered <- buildings |>
  filter(height > 50) |>           # dplyr filter (attribute)
  mutate(area = ddbs_area(.)) |>   # Add geometry measurement
  ddbs_filter(downtown,            # Spatial filter
              predicate = "within")

# Support standard sf methods
crs_info <- st_crs(buildings)
bbox <- st_bbox(buildings)
geom_col <- st_geometry(buildings)

Common Mistakes

Mistake	Fix
Forgetting ddbs_collect()	Data stays in DuckDB - call collect() to get sf
Not validating before union	Always ddbs_make_valid() before ddbs_union()
Wrong CRS for measurements	Transform to projected CRS before distance/area
Mixing sf and duckspatial ops	Convert with as_duckspatial_df() or collect()
Large dataset in collect()	Filter/aggregate in DuckDB first, then collect

Performance Characteristics

When duckspatial is faster:

• Spatial joins with >100K features per layer
• Operations on files >500MB
• Repeated queries (DuckDB caches and optimizes)
• Chained spatial operations (single optimized query plan)

When sf is faster:

• Small datasets (<100K features)
• Operations requiring immediate results
• Complex GEOS operations not in DuckDB spatial extension

Memory pattern:

# BAD: Loads entire file into R memory
data <- st_read("huge_file.gpkg") |>
  st_filter(bbox)

# GOOD: Filters in DuckDB, only loads result
data <- ddbs_open_dataset("huge_file.gpkg") |>
  ddbs_intersects_extent(bbox) |>
  ddbs_collect()

Gotchas

1. First call overhead: Initial ddbs_open_dataset() takes seconds for connection setup
2. Lazy evaluation: Operations don't execute until collect() - errors appear late
3. Geometry validity: Invalid geometries cause cryptic DuckDB errors - validate early
4. CRS handling: DuckDB spatial extension uses EPSG codes - proj4 strings may not work
5. Non-persistent: In-memory DuckDB discards data after session unless using file-backed connection

Spatial Predicates Reference

Predicate	Meaning	Use Case
intersects	Geometries share any space	Most permissive, default
within	x completely inside y	Points in polygons
contains	x completely contains y	Find polygons containing pt
covers	x covers y (boundary can touch)	Relaxed containment
touches	Share boundary but not interior	Adjacent polygons
overlaps	Share space but neither contains other	Partial overlap
crosses	Geometries cross (lines)	Line intersections
equals	Spatially identical	Exact match
disjoint	No spatial interaction	Exclusion filtering

Real-World Impact

Before (sf): 45 minutes to spatial join 2M buildings to 500 neighborhoods, 32GB RAM After (duckspatial): 3 minutes, 4GB RAM

Pattern: Operations scale sub-linearly with data size due to DuckDB's columnar storage and parallel execution.