为什么apply似乎放弃了CRS？

Question

标题确实总结了我的问题。Apply 似乎会删除 CRS，但其他函数不会。计算点向量上的地理函数的最佳方法是什么？

library(tidyverse)
library(sf)

# Generate 1000 lat longs, save as sf, and set crs
df1 <- data.frame(lat = runif(1000, 30, 33.4), long = runif(1000, -95, -82)) %>%
  st_as_sf(coords = c("long", "lat"),
           crs = 4326)

# Single point, with identical crs
df2 <- data.frame(lat = 32, long = -96) %>%
  st_as_sf(coords = c("long", "lat"),
           crs = 4326)

apply(df1, 1, function(x) st_distance(x, df2))

这给出了错误：Error in st_distance(x, df2) : st_crs(x) == st_crs(y) is not TRUE

但这些都工作得很好：

st_distance(df1[1,], df2)

final.df <- NULL
for(i in 1:nrow(df1)){
  ith.distance <- st_distance(df1[i,], df2)
  
  final.df <- rbind(final.df, ith.distance)
}

for 循环当然不是最有效的方法。...是吗？？

Answer 1

要点：

1. for循环apply在性能上没有真正的差异，这是一个多年前修复的旧条件（它曾经是正确的）。在大多数情况下，使用其中一种而不是另一种的决定应该基于数据争论和熟练程度的原因，而不是“速度”。通常，当循环性能出现问题时for，它与循环内部所做的事情有关，而与循环for本身无关。

2. 从长远来看，迭代调用的final.df <- rbind(final.df, ith.distance)性能会很差。每次迭代时，它都会生成所有数据的完整副本；当您有几行时，这很好，但是当您进入 100 和 1000 行时，每次循环都需要多一点时间。不要这样做，最好将结果附加到 alist并对 进行一次调用rbind。稍后我将演示这一点。

3. 通常循环较小的帧会更有效；在这种情况下，这将df2代替df1. （也许这只是您生成的示例数据的问题，但我是为了以防万一而这么说的。）

我建议迭代行号并以这种方式进行距离计算。

### iterating over the LARGER frame
system.time(res1 <- sapply(seq_len(nrow(df1)), function(rn) as.numeric(sf::st_distance(df1[rn,], df2))))
#    user  system elapsed 
#   4.782   0.000   4.783 

### iterating over the SMALLER frame
system.time(res2 <- sapply(seq_len(nrow(df2)), function(rn) as.numeric(sf::st_distance(df1, df2[rn,]))))
#    user  system elapsed 
#   0.008   0.000   0.009 

identical(res1, c(res2))
# [1] TRUE
head(res1)
# [1]  974915.8  184703.5 1215753.1  260598.1 1063713.6  745808.6
head(res2)
#           [,1]
# [1,]  974915.8
# [2,]  184703.5
# [3,] 1215753.1
# [4,]  260598.1
# [5,] 1063713.6
# [6,]  745808.6

res1是一个向量，res2是一个矩阵（1 列），并且值相同。后者运行得更快，这是在可用时使用 R 向量化计算的结果（它们在这里）。

作为记录，这里是一个for循环，但避免了 per-pass rbind。

lst1 <- list()
system.time({
  for (rn in 1:nrow(df1)) {
    dist <- sf::st_distance(df1[rn,], df2)
    lst1 <- c(lst1, list(dist))
  }
  lst1 <- do.call(rbind, lst1)
})
#    user  system elapsed 
#   4.717   0.000   4.717 

lst2 <- list()
system.time({
  for (rn in 1:nrow(df2)) {
    dist <- sf::st_distance(df1, df2[rn,])
    lst2 <- c(lst2, list(dist))
  }
  lst2 <- do.call(rbind, lst2)
})
#    user  system elapsed 
#   0.010   0.000   0.011 

identical(lst1, lst2)
# [1] TRUE
identical(res2, lst2)
# [1] TRUE

至于为什么apply不行：我们看看内部提供了什么：

apply(df1, 1, function(x) {browser();1;})
Browse[1]> debug at #1: [1] 1
Browse[2]> x
$geometry
POINT (-85.94798 30.29108)

apply如果你看一下( )的开头src/library/base/R/apply.R#29，它所做的第一件事就是将输入转换为 a matrix：

    ## Ensure that X is an array object
    dl <- length(dim(X))
    if(!dl) stop("dim(X) must have a positive length")
    if(is.object(X))
    X <- if(dl == 2L) as.matrix(X) else as.array(X)
###                   ^^^^^^^^^^^^  this is a problem for sf-objects

除此之外，这会从"sf"-class 对象中剥离 CRS。

（郑重声明，在对具有混合列和非字符串列的as.matrix(X)帧执行逐行操作时，对 的调用也是一个问题；如果在 之前没有充分子集化，则所有值都可能会转换为字符串。）characterXapply