基于空间设置更新多边形内点的属性

问：

假设我有如下所示的多边形和点：

library(sf)
library(sp)

####Create minimum working example data####
n_rows <- 10  # Number of rows of rectangles
n_cols <- 10  # Number of columns of rectangles
rect_width <- 1  # Width of each rectangle
rect_height <- 1  # Height of each rectangle

#Initialize a list to store rectangle polygons
rectangles <- list()

#Create rectangles and add them to the list with gaps
#Create rectangles and add them to the list
for (i in 1:n_rows) {
  for (j in 1:n_cols) {
    #Define the coordinates for the four corners of each rectangle
    x_coords <- c((j - 1) * rect_width, j * rect_width, j * rect_width, (j - 1) * rect_width, (j - 1) * rect_width)
    y_coords <- c((i - 1) * rect_height, (i - 1) * rect_height, i * rect_height, i * rect_height, (i - 1) * rect_height)
    
    #Combine the coordinates into a matrix
    rect_coords <- cbind(x_coords, y_coords)
    
    #Create a polygon from the coordinates and add it to the list
    rectangles[[length(rectangles) + 1]] <- st_polygon(list(rect_coords))
  }
}

#Create an sf object representing the rectangles with gaps
polygons_sf <- st_sf(geometry = st_sfc(rectangles))

polygons_sf$ID <- 1:100  # Add an ID column to represent each rectangle

#Generate random points within each rectangle
set.seed(123)  # for reproducibility
points <- st_sample(polygons_sf, size = 100)
points_sf <- st_sf(geometry = st_sfc(points))

#Generate a vector of values 
Depth <- as.integer(runif(100, min = 1, max = 5))

#Add the new column to random_points
points_sf$Z <- Depth

对于每个多边形中的点，我想做一些计算来改变点的深度。对于计算，我需要收集多边形周围多边形内点的平均深度半径内的所有点。例如，让我们对多边形 #45 执行此操作：

# Function to calculate average Z value for points within a polygon
calculate_average_Z <- function(polygon, points) {
  points_in_polygon <- st_intersection(points, polygon)  
  if (nrow(points_in_polygon) == 0) {
    return(NA)  # No points in the polygon
  }
  avg_Z <- mean(points_in_polygon$Z)
  return(avg_Z)
}

#Do select points in a buffer around polygon #45
polygon <- polygons_sf[45, ]
    
avg_Z <- calculate_average_Z(polygon, points_sf)
    
#if (is.na(avg_Z)) {
 # next 
#}
    
buffer_distance <- avg_Z  
buffered_polygon <- st_buffer(polygon, dist = buffer_distance)
    
# Select points within the buffered area
points_in_buffered_area <- st_filter(points_sf, buffered_polygon)

现在我想根据多边形的属性和缓冲区内点的 X/Y/Z 来计算一些东西，这将导致多边形内点的深度发生变化：

##For the real calculations I need X and Y as columns
#point_coordinates <- st_coordinates(points_in_buffered_area)
#points_in_buffered_area$X <- point_coordinates[, "X"]
#points_in_buffered_area$Y <- point_coordinates[, "Y"]

#while loop for changing depth of the points     
    while (TRUE) {
      Depth_t <- 4
      ##here there is in reality a very complicated function
      #that uses the properies of the polygon and the X/Y/Z of the points as input
      #to calculate the depth needed at this location. 
      Depth <- avg_Z
      
    if (Depth>= Depth_t) {
      break  # Exit the loop if the condition is met
    }
    
    avg_Z <- avg_Z + 1
    }

现在，我想更新原始 points_sf 数据帧中多边形内点的 Z 值：

new_Z <- data.frame(X = points_in_buffered_area$X, Y = points_in_buffered_area$Y, new_Z = avg_Z, geom=points_in_buffered_area$geom)
# Convert the data frame to an sf object
new_Z_values <- st_as_sf(new_Z, crs = st_crs(points_sf))

# Select points within original polygon
new_Z_to_use <- st_intersection(new_Z_values, polygon)
new_Z_to_use <- new_Z_to_use[, -c(1,2,4)]
#print(new_Z_to_use)
# Join new Z values to points_sf and update Z values
points_sf <- st_join(points_sf, new_Z_to_use, by = c("geometry"= "geometry"))
points_sf$Z <- ifelse(!is.na(points_sf$new_Z), points_sf$new_Z, points_sf$Z)

#delete not needed column
points_sf <- points_sf[, -2]

当我手动执行此操作时，这非常有效，因为现在多边形 #45 内的点现在已经更新了 Z 值。但是，我希望这是自动化的，并尝试了以下方法：

# Create a vector of polygon indices
polygon_indices <- 1:nrow(polygons_sf)

# Initialize a list to store results
results_list <- vector("list", length = length(polygon_indices))

#function to loop through all polygons:
process_polygons <- function(indices) {
  points_sf<-points_sf
  for (i in indices) {
      polygon <- polygons_sf[i, ]
      avg_Z <- calculate_average_Z(polygon, points_sf)
    
      if (is.na(avg_Z)) {
         next
      }
    
      buffer_distance <- avg_Z  
    
      buffered_polygon <- st_buffer(polygon, dist = buffer_distance)
    
      # Select points within the buffered area
      points_in_buffered_area <- st_filter(points_sf, buffered_polygon)
      point_coordinates <- st_coordinates(points_in_buffered_area)
      points_in_buffered_area$X <- point_coordinates[, "X"]
      points_in_buffered_area$Y <- point_coordinates[, "Y"]
    
      while (TRUE) {
        Depth_t <- 4
        ##here there is in reality a very complicated function
        #that uses the properies of the polygon and the X/Y/Z of the points as input
        #to calculate the depth needed at this location. 
        Depth <- avg_Z
   
      if (Depth>= Depth_t) {
        break  # Exit the loop if the condition is met
      }
    
      avg_Z <- avg_Z + 1
      }
    
      new_Z <- data.frame(X = points_in_buffered_area$X, Y = points_in_buffered_area$Y, new_Z = avg_Z, geom=points_in_buffered_area$geom)
    
      # Convert the data frame to an sf object
      new_Z_values <- st_as_sf(new_Z, crs = st_crs(points_sf))
    
      # Select points within original polygon
      new_Z_to_use <- st_intersection(new_Z_values, polygon)
      new_Z_to_use <- new_Z_to_use[, -c(1,2,4)]
      #print(new_Z_to_use)
      # Join new Z values to points_sf and update Z values
      points_sf <- st_join(points_sf, new_Z_to_use, by = c("geometry"= "geometry"))
      points_sf$Z <- ifelse(!is.na(points_sf$new_Z), points_sf$new_Z, points_sf$Z)
      points_sf <- points_sf[, -2]
      # Store the result in the list
      result <- c(avg_Z)
     #print(result)
     return(result)
  }
}  

results_list <- lapply(polygon_indices, process_polygons)

我觉得这应该有效，但它不会使用更新的 Z 值更新原始points_sf数据帧。每次缓冲新多边形时，它都应该考虑前一个多边形的新 Z 值，但现在不会发生这种情况。知道如何（有效地）解决这个问题吗？

对不起，这篇文章很长，但我试图让它尽可能全面和可行。