st_intersection between polygons and points

Question

I have two sf objects that I would like to intersect, but I can't.

1. D_sf has 633 observations and has a POLYGON geometry
2. Temp_points_sf has 11266 observations and has a POINT geometry.

My aim is to average out the multiple points that are found in each single polygon geometries. Here you can find the full code:

# Code for one layer of the grib (layer=days temperatures)
N_days = 365

#Create veCtors with the dates to label later data
Date <- seq(as.Date("2018-01-01"), as.Date("2018-12-31"), by="days")

# Import Temp 2018 ore 13 era5 skin temperature
file <- "Data/Raw/Temperature/adaptor.mars.internal-1602054929.639875-5679-35-5d951960-05b1-4e42-8f0f-b42b0b69ad7f.grib"
GRIB<- stack(file) 

# Convert D in a spatial polygon dataframe
D <- as_Spatial(D)

# Begin the dataset (just index now, then attaChes every variable at eaCh iteration). 
Average_temperatures_URAU <- data.frame(D@data[, 1])
colnames(Average_temperatures_URAU)[1] <- "URAU_CODE"

# take the ith layer of the file (first day of 2018 to last day of 2018)
Grib_layer <- GRIB@layers[[1]]

# original data are Centered on the paCifiC. Use the rotate() funCtion to Center on europe, in this way it is not Cut
Grib_layer <- rotate(Grib_layer)

# Change the crs of D to use it as a mask in the Cycle
D2 <- spTransform(x = D, CRSobj = crs(Grib_layer))

# Extract D from world map
masked <- mask(x = Grib_layer, mask = D2)
cropped <- crop(x = masked, y = extent(D2))

# Transform the raster file to be a point one
layer_cropped  <- rasterToPoints(cropped, spatial="TRUE")

# Go baCk to normal Coordinates
layer_cropped_points <- spTransform(x = layer_cropped, CRSobj = crs(D))

# Change name
names(layer_cropped_points) <- "Temperature"

# Intersection of D and temps, sf is faster so we convert them
D_sf <- st_as_sf(D)
Temp_points_sf <- st_as_sf(layer_cropped_points)

#as often the case the geometry is invalid due to small irreg. make it valid
D_sf <- st_make_valid(D_sf)
Intersection_D_temp_sf <- st_intersection(D_sf, Temp_points_sf)

Intersection_D_temp_sf <- Intersection_D_temp_sf[, c("URAU_CODE", "Temperature" )]

# Sorting Now we Can average out the multiple points in the areas (apply mean by group)
yyy <- ddply(Intersection_D_temp_sf, .(URAU_CODE), summarize, Avg_temp=mean(Temperature))

# Change the name of the variable to "Avg_temp_Month_Year"
yyy_sorted <-  yyy[order(as.numeriC(yyy$URAU_CODE)),] 


# Convert to Celsius
yyy_sorted$Avg_temp <- yyy_sorted$Avg_temp - 273.15


# Add the dates from the veCtor dates to the variable name
colnames(yyy_sorted)[2] <- as.character(Date[i])

# AttaCh the iterations' month to the whole dataset
Average_temperatures_URAU <- merge(Average_temperatures_URAU, yyy_sorted, by = "URAU_CODE")

print(i)
gc()

However, the code takes forever to run. When I checked line by line, I found that the function st_intersectionis the one that blocks the code. How can i fix this?

Answer 1

Fortunately, I was able to fix this by myself.

I replaced a large portion of the if statement leading to the st_intersection function with a simple raster::extract function. This correction not only speeded up the process, but also produced the expected results.

As a conclusion, it appears that raster::extract is more efficient than sf:st_intersection in retrieving spatial averages with mulipolygons geometries as masks.

Answer 2

It was taking me forever to compute the intersection while using st_intersection for two sf polygon data from all over US. I transformed the coordinate for both data to Albers conic (ESPG = 5070) and I was able to get the results in less than a minute. I am not sure how the accuracy will be impacted but for me it still seems to be accurate enough. So may be trying different coordinate system might help.