I'm trying to implement this query that speeds up st_intersections by using faster functions, namely ST_Contains and ST_Within, to do ST_Intersection only where absolutely necessary, from PostGIS to R
SELECT CASE
WHEN ST_Contains(a.geom, b.geom)
THEN b.geom
WHEN ST_Within(a.geom, b.geom)
THEN a.geom
ELSE ST_Intersection(a.geom, b.geom)
END ...
FROM a JOIN b
ON ST_Intersects(a.geom, b.geom)
what I have so far:
library(sf)
library(dplyr)
library(raster)
# anp_fed = raster::getData(country = "MEX", level = 1) %>% st_as_sf() %>% filter(NAME_1 == "Veracruz")
# veg = read_sf("usv250ks4gw.shp")
# downloadable from http://www.conabio.gob.mx/informacion/gis/maps/geo/usv250ks4gw.zip
contienen <- st_contains(anp_fed, veg, sparse = T)
contenidas <- st_within(anp_fed, veg, sparse = T)
contienen_list <- c()
for(i in 1:dim(contienen)[1]) {
contienen_list <- c(contienen[[i]], contienen_list)
}
contienen_list <- sort(contienen_list)
contenidas_list <- c()
t_contenidas <- t(contenidas)
for(j in 1:dim(t_contenidas)[1]) {
contenidas_list <- c(t_contenidas[[j]], contenidas_list)
}
contenidas_list <- sort(contenidas_list)
all_in_one <- rbind(st_join(anp_fed[contenidas_list,], veg),
st_join(veg[contienen_list,], anp_fed) %>% as.data.frame() %>%
dplyr::select(names(anp_fed)[-length(names(anp_fed))], names(veg)) %>% st_as_sf(),
st_intersection(anp_fed[ifelse(length(contenidas_list) == 0, 1:dim(anp_fed)[1], -contenidas_list),],
veg[-contienen_list,])
)
# system.time(source("postgis_to_sf.R")) = 491 s
# system.time(st_intersection(anp_fed, veg)) = 666 s
I'm trying to do this without calling rpostgis because I want to do it directly on files in my system without having to put them into Postgresql; I think there's a lot of space for improvement, maybe there's a more straightforward way of translating sgbp to a vector to filter sf objects, and also the re-ordering of columns of the second st_join() call, transforming to data.frame; I'm trying to make a more concise code; note also the evaluation for empty sf objects needed (it should be done in every object, if it is to become a function)
Execution with st_intersection() takes 666 s, while the optimized version takes 491 s, almost 3 mins less.