使用apply()函数与简单要素(SF)函数配合使用

Question

我已经编写了一个函数，用于计算质心到其多边形边缘的最大距离，但我不知道如何在简单要素（“sf”）数据框的每个单独多边形上运行它。

library(sf)

distance.func <- function(polygon){
  max(st_distance(st_cast(polygon, "POINT"), st_centroid(polygon)))
}

如果我对单个多边形进行测试，该函数有效。(警告消息与当前问题无关)。

nc <- st_read(system.file("shape/nc.shp", package="sf")) # built in w/package
nc.1row <- nc[c(1),] # Just keep the first polygon

>distance.func(nc.1row)
24309.07 m
Warning messages:
1: In st_cast.sf(polygon, "POINT") :
   repeating attributes for all sub-geometries for which they may not be constant
2: In st_centroid.sfc(st_geometry(x), of_largest_polygon = of_largest_polygon) :
   st_centroid does not give correct centroids for longitude/latitude data

问题是将此函数应用于整个数据框。

nc$distance <- apply(nc, 1, distance.func)
Error in UseMethod("st_cast") :
 no applicable method for 'st_cast' applied to an object of class "list"

我该怎么做才能对类“sf”对象中的每个单独多边形运行此函数（或类似函数）？

- John J.

---

一个好的旧式for循环怎么样？ dist = list(); for (i in seq_along(nc)) dist[[i]] <- distance.func(nc[i,]) - lbusett

你说得对。那个方法完全可行。如果你把它作为答案留下，我会接受的。谢谢！ - John J.

Answer 1

问题在于直接在sf对象上使用类似apply的函数是“有问题的”，因为几何列是一个列表列，这与“apply”结构不良互动。

最简单的解决方法可能是只使用for循环：

library(sf)

nc <- st_read(system.file("shape/nc.shp", package="sf")) %>% 
  st_transform(3857)

distance.func <- function(polygon){
  max(st_distance(st_cast(polygon, "POINT"), st_centroid(polygon)))
}
dist <- list()
for (i in seq_along(nc[[1]])) dist[[i]] <- distance.func(nc[i,]) 

head(unlist(dist))
# [1] 30185.34 27001.39 34708.57 52751.61 57273.54 34598.17

但它的速度相当慢。

为了能够使用类似于apply的函数，您需要将对象的仅几何列传递给函数。像这样的东西会起作用：

library(purrr)

distance.func_lapply <- function(polygon){
  polygon <- st_sfc(polygon)
  max(st_distance(st_cast(polygon, "POINT"), st_centroid(polygon)))
}

dist_lapply <- lapply(st_geometry(nc),  distance.func_lapply)
dist_map    <- purrr::map(st_geometry(nc), distance.func_lapply)

all.equal(dist, dist_lapply)
# [1] TRUE
all.equal(dist, dist_map)
# [1] TRUE

请注意，我必须略微修改距离函数，添加st_sfc调用，否则会出现许多“In st_cast.MULTIPOLYGON(polygon, "POINT") : point from first coordinate only”警告，并且结果不正确（我没有调查这个原因 - 显然，在sfg对象和sfc对象上，st_cast的行为不同）。

就速度而言，无论是lapply解决方案还是map解决方案都比for循环快近一个数量级：

microbenchmark::microbenchmark(
  forloop = {for (i in seq_along(nc[[1]])) dist[[i]] <- distance.func(nc[i,])}, 
  map     = {dist_map <- purrr::map(st_geometry(nc), distance.func_lapply)},  
  lapply  = {dist_lapply <- lapply(st_geometry(nc),  distance.func_lapply)}, times = 10)

Unit: milliseconds
    expr      min       lq     mean   median       uq       max neval
 forloop 904.8827 919.5636 936.2214 920.7451 929.7186 1076.9646    10
     map 122.7597 124.9074 126.1796 126.3326 127.6940  128.7551    10
  lapply 122.9131 125.3699 126.9642 126.8100 129.3791  131.2675    10

- lbusett

---

我是艰难地发现了这个事实。 虽然在R中，for循环速度较慢，但它们仍然在语言中有一席之地。例如，我可以添加一个漂亮的进度条来监视过程的进展。当使用*apply函数时，我发现这很难做到。 - Faustin Gashakamba

Answer 2

有另一种方法可以应用于简单特征，虽然并不比使用for循环更好。您可以先使用lapply创建一个简单特征列表，然后再应用距离函数。

distance.func <- function(polygon){
  max(st_distance(st_cast(polygon, "POINT"), st_centroid(polygon)))
}

distance.func.ls_sf <- function(sf){
  ls_sf <- lapply(1:nrow(sf), function(x, sf) {sf[x,]}, sf)
  dist <- lapply(ls_sf, distance.func)
}

dist_lapply_ls_sf <- distance.func.ls_sf(nc)

all.equal(dist, dist_lapply_ls_sf)
# [1] TRUE

性能几乎和for循环一样差...甚至在4年后（现在是R 4.1.1和sf 1.0-3），它的性能几乎比使用lapply或map和st_geometry(nc)差了近两个数量级（而不是一个）...

microbenchmark::microbenchmark(
  forloop = {for (i in seq_along(nc[[1]])) dist[[i]] <- distance.func(nc[i,])}, 
  map     = {dist_map <- purrr::map(st_geometry(nc), distance.func_lapply)},  
  lapply  = {dist_lapply <- lapply(st_geometry(nc),  distance.func_lapply)},
  ls_sf   = {dist_lapply_ls_sf <- distance.func.ls_sf(nc)},
  times = 10)

Unit: milliseconds
    expr       min        lq      mean    median        uq       max neval
 forloop 7726.9337 7744.7534 7837.6937 7781.2301 7850.7447 8221.2092    10
     map  124.1067  126.2212  135.1502  128.4745  130.2372  182.1479    10
  lapply  122.0224  125.6585  130.6488  127.9388  134.1495  147.9301    10
   ls_sf 7722.1066 7733.8204 7785.8104 7775.5011 7814.3849 7911.3466    10

所以，除非你应用于简单要素对象的函数计算时间比st_distance()更长，否则这是一个不好的解决方案。

如果您需要属性怎么办？

如果您的函数需要sf对象的几何和属性部分，使用mapply是一个很好的选择。以下是一个示例，使用三种方法计算突发性婴儿死亡密度（SID/km²）：

• for
• 在使用lapply之前提取每个要素
• mapply

microbenchmark::microbenchmark(
  forLoop =
    {
      sid.density.for <- vector("list", nrow(nc))
      for (i in seq(nrow(nc))) sid.density.for[[i]] <- nc[i,][["SID74"]] / st_area(nc[i,]) / 1000^2
    },
  list_nc =
    {
      list_nc <- lapply(seq(nrow(nc)), function(x, nc) { nc[x,] }, nc)
      sid.density.lapply <- lapply(list_nc, function(x) { x[["SID74"]] / as.numeric(st_area(x)) / 1000^2 })
    },
  mapply =
    {
      sid.density.func <- function(geometry, attribute) { attribute / st_area(geometry) / 1000^2 }
      sid.density.mapply <- mapply(sid.density.func, st_geometry(nc), nc[["SID74"]], SIMPLIFY = FALSE)
    },
  times = 10)

Unit: milliseconds
    expr       min        lq       mean     median        uq       max neval
 forLoop 4511.7203 4515.5997 4557.73503 4542.75200 4560.5508 4707.2877    10
 list_nc 4356.3801 4400.5640 4455.35743 4440.38775 4475.2213 4717.5218    10
  mapply   17.4783   17.6885   18.20704   17.99295   18.3078   20.1121    10

- Jean-François Bourdon