R语言st_intersection处理多边形与点相交运行过慢如何解决

我有两个sf对象执行相交操作时耗时极长：

• D_sf：633条POLYGON几何观测
• Temp_points_sf：11266条POINT几何观测
  目标是计算每个多边形内部所有对应点的温度平均值，原代码如下：

# 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()

逐行调试后发现st_intersection函数是性能瓶颈，请问该如何优化解决该问题？

核心问题 #

你完全不需要用st_intersection做面和点的匹配，st_intersection会计算几何相交后的新几何，对于点面归属场景属于过度计算，耗时自然很高。

第一层优化：替换点面匹配逻辑 #

直接用空间连接匹配点所属的多边形，性能提升至少10倍以上：

# 原st_intersection代码替换为以下内容
# 先给点数据匹配所属的多边形URAU_CODE
point_with_poly <- st_join(Temp_points_sf, D_sf["URAU_CODE"], join = st_within)
# 过滤掉不在任何多边形内的点
point_with_poly <- point_with_poly[!is.na(point_with_poly$URAU_CODE), ]
# 分组求平均
yyy <- aggregate(Temperature ~ URAU_CODE, data = point_with_poly, FUN = mean)

如果追求更快的速度，可以用st_intersects先获取匹配索引，直接计算均值，跳过属性连接的步骤：

# 获取每个多边形包含的点的索引
point_index <- st_intersects(D_sf, Temp_points_sf)
# 直接计算每个多边形的平均温度
yyy <- data.frame(
  URAU_CODE = D_sf$URAU_CODE,
  Avg_temp = sapply(point_index, function(idx) mean(Temp_points_sf$Temperature[idx], na.rm = TRUE))
)

第二层优化：跳过转点步骤，直接用栅格提取 #

你当前的逻辑是把栅格转成点再做匹配，完全可以直接用栅格的extract函数计算每个多边形的均值，省去转点和空间匹配两步，性能提升更高：

# 不需要执行rasterToPoints、转sf、空间匹配的相关代码
# 直接在裁剪后的栅格上提取每个多边形的均值
temp_mean <- extract(cropped, D2, fun = mean, na.rm = TRUE)
# 直接和结果表拼接即可
yyy_sorted <- data.frame(
  URAU_CODE = D$URAU_CODE,
  Avg_temp = temp_mean - 273.15
)

另外可以提前把D的sf对象、坐标转换这些操作放到循环外，不要每次迭代都重复执行，进一步减少冗余计算。

内容的提问来源于stack exchange，提问作者nflore