Р: Как заменить значения в геометрии многоугольника без цикла for?

Происходит что-то такое, что, возможно, кто-то другой сможет объяснить — возможно, весь столбец geometry перезаписывается при каждой итерации с использованием цикла for. Подход Map работает на порядки быстрее с более крупным объектом sf:

library(sf)

# load test data
myshape <- st_read(system.file("shape/nc.shp", package = "sf"))[sample(100, 1e4, 1),]

# make some random data
new_values <- runif (1e4)
myshape$geometry2 <- myshape$geometry3 <- myshape$geometry

system.time({
  for (row in 1:nrow(myshape)) {
    # replace
    myshape$geometry2[row][[1]][[1]][[1]][,2] <- new_values[row]
  }
})
#>    user  system elapsed 
#>  105.19    0.39  105.64

system.time({
  f <- function(g, x) {
    g[[1]][[1]][,2] <- x
    g
  }
  
  myshape$geometry3[] <- Map(f, myshape$geometry, new_values)
})
#>    user  system elapsed 
#>    0.07    0.00    0.06

identical(myshape$geometry2, myshape$geometry3)
#> [1] TRUE

Это относительно хорошо масштабируется. Повторное выражение основано на данных nc.shp, более крупный пример см. ниже. Рабочий процесс:

1. Создайте new_values ​​df со значениями z и соответствующим идентификатором объекта, например. nc$NAME
2. st_cast() myshape для многоугольника, назначайте pid отдельным многоугольникам, st_cast() для точек, объединяйте значения XYZ и создавайте геометрию, перестраивайте мультиполигоны

library(sf)
library(sfheaders)
library(dplyr)
library(ggplot2)

# Load test data
myshape <- st_read(system.file("shape/nc.shp", package = "sf"))

# Make some random z coordinate data with associated feature ID e.g. myshape$NAME
set.seed(1)
new_values <- data.frame(NAME = myshape$NAME,
                         z = runif (nrow(myshape)))

# Create output with z values assigned
sf_poly <- myshape %>%
  # Split multipolygons and assign unique id to each polygon for 
  st_cast("POLYGON") %>%
  mutate(pid = 1:n(), .by = NAME) %>%
  # Cast to point, retrieve lon/lat values, left_join() z values
  st_cast("POINT") %>%
  mutate(x = st_coordinates(.)[,1],
         y = st_coordinates(.)[,2]) %>%
  left_join(., new_values, by = "NAME") %>%
  # Drop original geometry, create new XYZ geometry, assign CRS
  st_drop_geometry() %>%
  sf_pt(., keep = TRUE) %>%
  rename(geometryz = geometry) %>%
  st_set_crs(st_crs(myshape)) %>%
  # Rebuild multipolygons in two steps, otherwise non-contiguous features will be joined
  summarise(geometryz = st_combine(geometryz), .by = c(NAME, pid)) %>%
  st_cast("POLYGON") %>%
  summarise(geometryz = st_combine(geometryz), .by = NAME) %>%
  st_cast("MULTIPOLYGON")

Результат из большего набора данных на ПК с низкими характеристиками:
Используя эти данные с 60 324 полигонами:

# Set extent
sfc <- st_sfc(st_polygon(list(rbind(c(0,0), c(20,0), c(20,20), c(0,0)))))

# Make polygon sf
sf_grid <- st_make_grid(sfc, cellsize = .088, square = FALSE) %>%
  st_as_sf() %>%
  mutate(NAME = paste0("poly", sprintf("%02d", 1:n()))) %>%
  st_set_crs(st_crs(myshape)) %>%
  rename(geometry = x)

Это заняло ~4,63 минуты. Не совсем репрезентативно, поскольку sf_grid не имеет функций мультиполигонов, но, тем не менее, обеспечивает приличное сокращение времени.