2

I'm trying to use st_within() to calculate points within polygons on a fairly complex shapefile. This operation tends to take 4-5 hours which is a bummer for debugging the code. I have access to a supercomputing cluster (120+ cores) and it seems sensible to make use of this tool to address this issue. However, it's unclear to me how I would even go about optimising these geospatial operations (from sf()) for multicore computing.

To get a sense of what I mean, you can try analysis using the UK National Forest Inventory: https://opendata.arcgis.com/datasets/bcd6742a2add4b68962aec073ab44138_0.zip?outSR=%7B%22wkid%22%3A27700%2C%22latestWkid%22%3A27700%7D

if (file.exists("data/National_Forest_Inventory_Woodland_Scotland_2017.shp") == FALSE) {
download.file("https://opendata.arcgis.com/datasets/3cb1abc185a247a48b9d53e4c4a8be87_0.zip?outSR=%7B%22wkid%22%3A27700%2C%22latestWkid%22%3A27700%7D", 
              destfile = "data/National_Forest_Inventory_Woodland_Scotland_2017.zip")
unzip("data/National_Forest_Inventory_Woodland_Scotland_2017.zip", exdir = "data")
}
forestinv <- st_read("data/National_Forest_Inventory_Woodland_Scotland_2017.shp") %>% st_transform(paste0("+init=epsg:",27700))

You can see a concrete example in code here (especially lines 507-510)

Here's the key question-is it possible to parallelise these kinds of geospatial operations in R such as st_within? Even with a simplified set of geometries this process takes ages, but there aren't obvious ways to loop the operation and split across cores.

6
  • 2
    The way to go is to split the complex geometries. PostGIS has a function for that postgis.net/docs/ST_Subdivide.html but I do not know how to do that with R.
    – user30184
    Mar 5 '20 at 21:46
  • I like the sound of that and I guess I could preprocess data in postgis, but sure would be nice if I could keep it all in R Mar 5 '20 at 21:51
  • 1
    I was able to calculate point-in-polygon (random points in 1:50m country polygons) using a single-threaded process against a memory-based geometry table at a rate of 50 microseconds/query. Massive parallel techniques are probably not necessary to get adequate performance, provided you "dice your Godzillas" (esri.com/arcgis-blog/products/arcgis-desktop/analytics/…).
    – Vince
    Mar 5 '20 at 21:59
  • I would move your data to PostGIS or at least make sure you have a spatial index and that it is being used in the R code
    – Ian Turton
    Mar 6 '20 at 9:04
  • 1
    Reverse your argument order and use st_contains instead of st_within. This allows internal indexes to be made on the polygons and should speed up your calculations by 10-20x. Some more details at github.com/r-spatial/sf/issues/1261 . Explicit subdivision, as used in PostGIS, should not be necessary or especially helpful in R, where all geometries are already in memory.
    – dbaston
    Mar 9 '20 at 15:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.