How to speed up the intersection of two large shapefiles in R with sf (or any other way)

Question

I have two large-ish shapefiles.

Shapefile 1:

One can be downloaded from here:

https://data.statistik.gv.at/data/OGDEXT_DSR_1_STATISTIK_AUSTRIA_20111031.zip

• It is ~20 MB in size
• Seems to be super detailed
• It has three features (classes of settlement) and two columns

When I print the object it takes forever to load the output which looks roughly like this:

Simple feature collection with 3 features and 2 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 112294.6 ymin: 275328.4 xmax: 685584 ymax: 570580.4
Projected CRS: SOURCECRS

Shapefile 2:

• can be downloaded from here: https://data.statistik.gv.at/data/OGDEXT_GEM_1_STATISTIK_AUSTRIA_20220101.zip
• Is roughly 61 MB in size
• Has 2115 Features (2115 administrative areas in Austria) and two columns
• When I print it, it prints the output much faster (why?! when its bigger), which looks like this:

Simple feature collection with 2115 features and 2 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 112518.2 ymin: 275472 xmax: 685444.5 ymax: 570431.1
Projected CRS: SOURCECRS

Now I would like to find out how much area of each polygon is intersecting with each class in the first shapefile (the settlement layer). A normal st_intersection() takes way too long. It did not finish and I let it run for the entire day.

How do I speed this up, perhaps by chunking up some layers?

Answer 1

What you can also do is to use the parallel processing with smaller intersections. For example, iterating or using each core to do the process for one unit from your region layer.

This is the process:

1. Read layers.
2. Validate geometries (most of the time this open layers come with some problems that can ruin your operations).
3. Set parallel parameters.
4. Run the code (I did crop and then inetersection, but you can change it to use directly the intersection).
5. Merge results of parallel process.
6. Export.

Here is the code:

# PARALLLEL INTERSECTION OF LARGE LAYERS
library(sf)
library(mapview)
library(s2)
library(tictoc)
library(parallel)
library(doParallel)

# a download: https://data.statistik.gv.at/data/OGDEXT_DSR_1_STATISTIK_AUSTRIA_20111031.zip
# b download: https://data.statistik.gv.at/data/OGDEXT_GEM_1_STATISTIK_AUSTRIA_20220101.zip

# read data
a <- read_sf('./data/20111031/STATISTIK_AUSTRIA_DSR_20111031.shp')
b <- read_sf('./data/20220101/STATISTIK_AUSTRIA_GEM_20220101.shp')

# Check layers
# mapview(a)+b

tic()
# Correct topology (this is just a recommendation)
# This can take a while but it worth
a <- st_make_valid(st_cast(a, "POLYGON")) # transform into POLYGONS
b <- st_make_valid(b)
toc()

tic()
# Process in parallel (intersection by "village" (name) in b)
# Proceso de cálculo de para cada TIPO (MUNI, EEO, ESF, EEO_ESF)
cl = parallel::makeCluster(detectCores(), type="FORK")
doParallel::registerDoParallel(cl, detectCores())

# Union in parallel
names <- unique(b$name) #different munis
pols = foreach(i=1:length(names)) %dopar% {
    x <- b %>% dplyr::filter(name == names[[i]]) %>% st_make_valid()
    crop <- st_crop(a, x) %>% st_make_valid()

    # remove invalid polygons
    # this allows the process to continue. Sometimes topology is not perfectly 
    # made and tinny pols can ruin your process
    notvalid <- which(s2_is_valid_detail(crop)==FALSE)
    if(length(notvalid) > 0){crop <- crop[-notvalid,]}
    
    # intersection between x and eeoval and esfval
    x2 <- crop %>% st_intersection(x) %>% st_make_valid()
    
    # "dissolve" to get
    x3 <- x2 %>% group_by(NAME, id, name) %>% 
        summarize()
    
    # return
    x3
}

parallel::stopCluster(cl)
toc()

cast <- lapply(pols, function(x) st_cast(x, "MULTIPOLYGON"))

# merge all
# Merge polygons
ab <- sf::st_as_sf(data.table::rbindlist(cast)) # superfast
ab <- ab %>% rename(NAME_A = NAME, 
              NAME_B = name,
              ID = id)
mapview(ab)

# export
st_write(ab, "ab.gpkg", append = FALSE)

I have a AMD 7 3700X and it performed the process in +-400 seconds (including validation). The issue here is that during the process around 30 Gb of ram are in use. If you have less than 32 Gb, you will need to set some subsets and run the code in some other way.

Final note: the code is writen for linux parallel process (but you can tweak it for windows easily).

Here is the result:

Answer 2

By just changing the geometry to polygons before the intersection you'll gain a lot of time in the processing.

You can use st_cast for both shapefiles

st_cast(mynameshape, "POLYGON")