将 apply() 与简单特征 (SF) 函数一起使用
Use apply() with a simple features (SF) function
我已经编写了一个函数来计算质心与其多边形边缘之间的最大距离，但我无法弄清楚如何 运行 它在简单特征的每个单独多边形上（ "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"
对于 class "sf" 对象中的每个单独的多边形，我可以为 运行 这个函数（或类似函数）做些什么？
这里的问题是直接在 sf 对象上使用类似应用的函数是 "problematic" 因为几何列是一个列表列，它与 "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
，但是很慢。
为了能够使用类似应用的函数，您只需要将对象的几何列传递给函数。这样的事情会起作用：
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 对象上的行为不同）。
在速度方面，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
还有另一种方法可以应用于简单的特征，尽管并不比使用 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
我已经编写了一个函数来计算质心与其多边形边缘之间的最大距离，但我无法弄清楚如何 运行 它在简单特征的每个单独多边形上（ "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"
对于 class "sf" 对象中的每个单独的多边形，我可以为 运行 这个函数（或类似函数）做些什么？
这里的问题是直接在 sf 对象上使用类似应用的函数是 "problematic" 因为几何列是一个列表列，它与 "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
，但是很慢。
为了能够使用类似应用的函数，您只需要将对象的几何列传递给函数。这样的事情会起作用：
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 对象上的行为不同）。
在速度方面，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
还有另一种方法可以应用于简单的特征，尽管并不比使用 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