geosareal. Some Alternatives to sf Functions And This is a Working

Title

geosareal

geosareal provides fast, geos–backed spatial data operations for R,
with a focus on areal-weighted interpolation as a drop-in alternative to
sf::st_interpolate_aw() and areal::aw_interpolate().

Installation

You can install the development version of geosareal from GitHub with:

# install.packages("pak")
pak::pak("e-kotov/geosareal")

Example

This is a basic example which shows you how to solve a common problem:

## Quick example (using the default `sf` dataset)
library(sf)
library(geosareal)

# data & grid (from sf help)
nc <- st_read(system.file("shape/nc.shp", package = "sf"), quiet = TRUE)
g <- st_make_grid(nc, n = c(10, 5)) # 'to' can be sfc for sf, but geosareal core needs an sf with an ID

# Prepare IDs for geosareal core
nc$src_id <- seq_len(nrow(nc))
g_sf <- st_as_sf(g)
g_sf$tid <- seq_len(nrow(g_sf))

## 1) geosareal core (geos_interpolate_aw) --------------------------------------

# (a) Treat BIR74 as spatially intensive (not mass-preserving)
a1_core <- geos_interpolate_aw(
  .data = hex,
  tid = hex_id,
  source = pop_fr_vec,
  sid = cell_id,
  weight = "sum",
  output = "tibble",
  extensive = "population_density"
)
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 <- geos_interpolate_aw(
  .data = hex,
  tid = hex_id,
  source = pop_fr_vec,
  sid = cell_id,
  weight = "sum",
  output = "tibble",
  extensive = "population_density"
)
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)

To test speed on relatively large data

library(geos)
library(sf)
library(areal)
library(terra)
library(geodata)
library(areal)
library(tidyverse)
library(geosareal)
library(bench)


# download may take > 5s
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))
# plot(st_geometry(fr2_mainland), main = "France NUTS-2 regions (2021)")

# 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)
# plot(pop_fr, main = "GPWv4 Population Density in Mainland France (2020)")

# project the raster to EPSG:3035 (ETRS89 / LAEA Europe)
pop_fr <- terra::project(
  pop_fr,
  "EPSG:3035",
  method = "bilinear"
)


# vectorize the raster
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 1000 m
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)
format(object.size(hex), "Mb")


aw_test <- ar_validate(
  source = pop_fr_vec,
  target = hex,
  varList = "population_density",
  method = "aw",
  verbose = TRUE
)
aw_test


Sys.time()
tictoc::tic()
fr_areal <- aw_interpolate(
  .data = hex,
  tid = hex_id,
  source = pop_fr_vec,
  sid = cell_id,
  weight = "sum",
  output = "tibble",
  extensive = "population_density"
)
tictoc::toc()
Sys.time()


Sys.time()
tictoc::tic()
hex_dst <- geos_interpolate_aw(
  .data = hex,
  tid = hex_id,
  source = pop_fr_vec,
  sid = cell_id,
  weight = "sum",
  output = "tibble",
  extensive = "population_density"
)
tictoc::toc()
Sys.time()


Sys.time()
tictoc::tic()
bm <- bench::mark(
  sf_areal_aw = aw_interpolate(
    .data = hex,
    tid = hex_id,
    source = pop_fr_vec,
    sid = cell_id,
    weight = "sum",
    output = "tibble",
    extensive = "population_density"
  ),
  geosareal_aw = geos_interpolate_aw(
    .data = hex,
    tid = hex_id,
    source = pop_fr_vec,
    sid = cell_id,
    weight = "sum",
    output = "tibble",
    extensive = "population_density"
  ),
  geosareal_aw_matrix = geos_interpolate_aw_matrix(
    .data = hex,
    tid = hex_id,
    source = pop_fr_vec,
    sid = cell_id,
    weight = "sum",
    output = "tibble",
    extensive = "population_density"
  ),
  check = FALSE,
  iterations = 5
)
tictoc::toc()
Sys.time()
bm

# A tibble: 3 × 13
  expression               min   median `itr/sec` mem_alloc `gc/sec` n_itr  n_gc 
  <bch:expr>             <bch:> <bch:>     <dbl>   <bch:byt>   <dbl>  <int> <dbl>
1 sf_areal_aw           10.06s   11.04s    0.0896     369MB    0.842     5    47      
2 geosareal_aw           5.54s    6.26s    0.134      115MB    0.242     5     9      
3 geosareal_aw_matrix    5.46s    6.53s    0.160      136MB    0.351     5    11      
# ℹ 5 more variables: total_time <bch:tm>, result <list>, memory <list>,
#   time <list>, gc <list>