8

我编写了一个函数来计算质心与其多边形边缘之间的最大距离，但我不知道如何在简单特征（“sf）data.frame 的每个单独多边形上运行它。

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

问题是将此函数应用于整个data.frame。

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”类的对象中的每个单独的多边形运行这个函数（或类似的函数）？

4

2 回答 2

12

这里的问题是直接在sf对象上使用类似应用的函数是“有问题的”，因为几何列是一个列表列，它不能与“应用”构造很好地交互。

最简单的解决方法可能是只使用一个 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

，但速度很慢。

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

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”警告，结果不正确（我没有调查原因 - 显然 st_cast 在对象上的行为与在sfg对象上sfc的行为不同）。

在速度方面，thelapply和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
于 2018-01-30T15:48:27.270 回答
0

还有另一种应用简单功能的方法，尽管它并不比使用 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 对象的几何和属性部分，那么 usingmapply是一个不错的方法。以下是使用三种方法计算婴儿猝死密度 (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
于 2021-10-29T15:40:03.260 回答