This package provides some alternatives to sf functions, which are implemented using duckdb and geoarrow. So far the only implemented function is areal interpolation working similar to the one found in ‘areal’ package, but much faster. I have no idea where this project will go, use at your own risk.
ducksf provides fast, DuckDB Spatial–backed spatial data
operations for R,
with a current focus on areal-weighted interpolation as a drop-in
alternative to
sf::st_interpolate_aw() and areal::aw_interpolate().
Here’s what you can expect for speed gains and memory savings (measured on projected planar data):
As you can see from the figures, ducksf is much faster than than
areal::aw_interpolate() (and sf::st_interpolate_aw(), as it works
the same way, but is functionally more limited), especially for large
datasets. It also uses less memory, which is crucial when working with
large spatial datasets. Even though there is a cost of translating sf
objects to DuckDB tables and back, the overall performance is
significantly improved.
The misterious geosareal package is just a proof-of-concept package
that also implements areal interpolation, but using geos instead of
duckdb. It is not meant to be used in production and is not
maintained, but the code is available here:
https://github.com/e-kotov/r-geos-areal-weighted-interpolation-prototype.
You can install the development version of {ducksf} from
GitHub with:
# install.packages("pak")
pak::pak("e-kotov/ducksf")This is a basic example which shows you how to solve a common problem
(following the sf package’s
example -
our result will be a little different, as we project the data to NAD83 /
Conus Albers (EPSG:5070)):
## Quick example (using the default `sf` dataset)
library(sf)
library(ducksf)
# data & grid (from sf help)
nc <- st_read(system.file("shape/nc.shp", package = "sf"), quiet = TRUE)
# transform to NAD83 / Conus Albers (EPSG:5070)
nc <- nc |> st_transform(5070)
g <- st_make_grid(nc, n = c(10, 5)) # 'to' can be sfc for sf, but ducksf core needs an sf with an ID
# Prepare IDs for ducksf core
nc$src_id <- seq_len(nrow(nc))
g_sf <- st_as_sf(g)
g_sf$tid <- seq_len(nrow(g_sf))
## 1) ducksf core (dst_interpolate_aw) --------------------------------------
# (a) Treat BIR74 as spatially intensive (not mass-preserving)
a1_core <- dst_interpolate_aw(
target_sf = g_sf,
tid = "tid",
source_sf = nc,
sid = "src_id",
weight = "sum", # intensive semantics
intensive = "BIR74",
output = "sf"
)
sum(st_drop_geometry(a1_core)$BIR74, na.rm = TRUE) / sum(nc$BIR74, na.rm = TRUE)
# (b) Treat BIR74 as spatially extensive (mass-preserving)
a2_core <- dst_interpolate_aw(
target_sf = g_sf,
tid = "tid",
source_sf = nc,
sid = "src_id",
weight = "total", # extensive semantics (pycnophylactic)
extensive = "BIR74",
output = "sf"
)
sum(st_drop_geometry(a2_core)$BIR74, na.rm = TRUE) / sum(nc$BIR74, na.rm = TRUE)
# Quick plot (intensive vs extensive)
a_show <- a1_core[, "BIR74"]
names(a_show)[1] <- "intensive"
a_show$extensive <- st_drop_geometry(a2_core)$BIR74
plot(a_show[c("intensive", "extensive")], key.pos = 4)
{ducksf} also provides a way to mask {sf} functions so that they
dispatch to {ducksf} instead of {sf}. This is useful if you want to
use {sf} functions with {ducksf} semantics without changing your
code too much. By default it is of course disabled, so you can use
{sf} functions as usual even if {ducksf} is loaded. The resulting
numbers will be insignificantly different from {sf} results due to
rounding errors, but they will be very close.
## 2) Masked drop-in (sf semantics) -----------------------------------------
# Enable/disable masking so st_interpolate_aw() dispatches to ducksf
sf_use_ducksf(TRUE) # or: sf_use_ducksf(TRUE, allow = "st_interpolate_aw")
sf_use_ducksf(FALSE) # or: sf_use_ducksf(TRUE, allow = "st_interpolate_aw")
a1 <- st_interpolate_aw(nc["BIR74"], g, extensive = FALSE)
sum(a1$BIR74) / sum(nc$BIR74) # not close to 1 (intensive)
a2 <- st_interpolate_aw(nc["BIR74"], g, extensive = TRUE)
sum(a2$BIR74) / sum(nc$BIR74) # ≈ 1 (mass-preserving)
# Compare maps exactly like sf's example
a1$intensive <- a1$BIR74
a1$extensive <- a2$BIR74
plot(a1[c("intensive", "extensive")], key.pos = 4)
# Disable masking when done
sf_use_ducksf(FALSE)This data and code were used for the tests presented at the top of this
README, but with bench::mark() function
install.packages(setdiff(
c(
"duckdb",
"sf",
"areal",
"geoarrow",
"terra",
"geodata",
"giscoR",
"tidyverse",
"tictoc"
),
rownames(installed.packages())
))
library(duckdb)
library(sf)
library(areal)
library(geoarrow)
library(terra)
library(geodata)
library(giscoR)
library(areal)
library(tidyverse)
library(tictoc)
library(ducksf) # remotes::install_github("e-kotov/ducksf")
# also optional if you want to try it yourself, geosareal also does quite well compared to sf-based interpolation
# library(geosareal) # remotes::install_github("e-kotov/r-geos-areal-weighted-interpolation-prototype")
# get data
dir.create("private/data", recursive = TRUE)
pop <- geodata::population(year = 2020, res = 2.5, path = "private/data") # GPWv4 density
nuts2 <- giscoR::gisco_get_nuts(year = "2021", nuts_level = 2, epsg = "4326")
fr2 <- subset(nuts2, CNTR_CODE == "FR")
fr2_mainland <- subset(fr2, !grepl("^FRY|^FRM", NUTS_ID))
# crop pop raster to mainland France
pop_fr <- terra::crop(
pop,
st_bbox(fr2_mainland),
snap = "out"
)
# mask
pop_fr <- terra::mask(pop_fr, fr2_mainland)
# project the raster to EPSG:3035 (ETRS89 / LAEA Europe)
pop_fr <- terra::project(
pop_fr,
"EPSG:3035",
method = "bilinear"
)
# vectorize the raster
# yes, let's imagine there is no terra::extract() or exactextractr::exact_extract()
# and you are bound to the vector world
pop_fr_vec <- as.polygons(pop_fr, aggregate = FALSE) |>
st_as_sf() |>
mutate(cell_id = as.character(row_number()))
# generate a hexagonal grid over mainland France
bbox_fr <- st_bbox(st_transform(fr2_mainland, 3035))
cellsize_hex <- 2500 # 2.5 km hexagons
# also try pushing it down to 1km hexes to push the limits of your computer
hex <- st_make_grid(bbox_fr, cellsize = cellsize_hex, square = FALSE) |>
st_as_sf() |>
st_set_geometry("geometry") |>
mutate(hex_id = as.character(row_number()))
nrow(hex) # 171k, would be just over 1 million for 1km hexes
format(object.size(hex), "Mb") # ~200 mb, would be 1,200 mb for 1km hexes
aw_test <- areal::ar_validate(
source = pop_fr_vec,
target = hex,
varList = "population_density",
method = "aw",
verbose = TRUE
)
aw_test
Sys.time()
tictoc::tic()
hex_areal <- areal::aw_interpolate(
.data = hex,
tid = hex_id,
source = pop_fr_vec,
sid = cell_id,
weight = "total",
output = "tibble",
extensive = "population_density"
)
tictoc::toc()
Sys.time()
Sys.time()
tictoc::tic()
hex_ducksf <- ducksf::dst_interpolate_aw(
target_sf = hex,
tid = "hex_id",
source_sf = pop_fr_vec,
sid = "cell_id",
weight = "total",
output = "tibble",
extensive = "population_density"
)
tictoc::toc()
Sys.time()
Sys.time()
tictoc::tic()
hex_geos <- geosareal::geos_interpolate_aw(
.data = hex,
tid = hex_id,
source = pop_fr_vec,
sid = cell_id,
weight = "total",
output = "tibble",
extensive = "population_density"
)
tictoc::toc()
Sys.time()