R（或Python）有什么方法可以加快地理空间分析？答案

【问题标题】：R (or Python) any way to speed up geospatial analysis?R（或Python）有什么方法可以加快地理空间分析？
【发布时间】：2018-01-27 19:08:26
【问题描述】：

我正在解决一个问题，我正在确定彼此数据点之间的距离内的数据点的性质。基本上，对于每一行数据，我都会尝试确定一个地理范围内数据点的“邻域”，然后找出这个“邻域”的特征。

问题是这是 O^2 问题，因为我目前嵌套了 for 循环，这意味着我正在运行 nrow^2 计算（我有 70k 行，所以 4.9B！计算......哎哟）

所以我拥有的 R（伪）代码是

for (i in 1:n.geopoints) {
   g1<-df[i,]
   for (j in 1:n.geopoints) {
      g2 <- df[j,]
      if (gdist(lat.1 = g1$lat, lon.1=g1$lon, lat.2 = g2$lat, lon.2 = g2$lon, units = "m") <= 1000) {
         [[[DO SOME STUFFF]]]
      }
   }
}

如何以更直接的方式实现这一点？有我可以依靠的功能吗？我可以矢量化吗？

我在 R 中有这个，但如果有更好的函数可用，可以很容易地将它放到 Python 中。

谢谢

【问题讨论】：

您希望它以多快的速度执行？您的目标是为每个点找到小半径（邻域）内的所有点吗？
是的。我不需要它快速运行，但我在本地钻机上运行它的第 20 小时.. 很想在 5 小时内完成它。我知道嵌套的 for 循环是一个容易实现的目标，但我一直在努力寻找一种方法来消除它。
尝试查看列表以外的数据结构。 k-d trees，例如。点云库肯定有适合您问题的东西（并且执行速度非常快）。

标签： python r for-loop optimization geospatial

【解决方案1】：

您可能需要考虑不同的方法。我会使用像 QGIS 这样的 GIS 工具来分割您的数据。就像你说的，你不需要数据的完整笛卡尔连接，只需要本地集群。看看一些聚类问题。

GIS Stackexchange 上的这个问题解决了一个具有 800k 数据点的类似类型问题。 https://gis.stackexchange.com/questions/211106/clustering-points-polygons-based-on-proximity-within-specifed-distance-using-q

【讨论】：

谢谢你——这为我指明了正确的方向。我使用了固定距离窃听器，然后在 PostGIS 中使用空间查询来总结与每条缓冲线相交的所有点的特征。与笛卡尔关节的 48 小时以上相比，处理能力只需要 30 分钟左右。

【解决方案2】：

这是使用data.table 的一种方法，以及我为this question 所做的重写的haversine 公式，以便它可以在data.table 操作中工作

这个想法是在每个点上对每个点进行data.table 连接，但在连接内计算每对点之间的距离，并删除那些超出阈值的点。这是受到@Jaap 优秀的answer here的启发

设置

半正弦公式是

## Haversine formula
dt.haversine <- function(lat_from, lon_from, lat_to, lon_to, r = 6378137){
  radians <- pi/180
  lat_to <- lat_to * radians
  lat_from <- lat_from * radians
  lon_to <- lon_to * radians
  lon_from <- lon_from * radians
  dLat <- (lat_to - lat_from)
  dLon <- (lon_to - lon_from)
  a <- (sin(dLat/2)^2) + (cos(lat_from) * cos(lat_to)) * (sin(dLon/2)^2)
  return(2 * atan2(sqrt(a), sqrt(1 - a)) * r)
}

我在此示例中使用的数据来自我的 googleway 包，它是墨尔本 City Loop 电车上的电车站

library(googleway)

## Tram stops data
head(tram_stops)
#   stop_id                                     stop_name stop_lat stop_lon
# 1   17880           10-Albert St/Nicholson St (Fitzroy) -37.8090 144.9731
# 2   17892    10-Albert St/Nicholson St (East Melbourne) -37.8094 144.9729
# 3   17893 11-Victoria Pde/Nicholson St (East Melbourne) -37.8083 144.9731
# 4   18010    9-La Trobe St/Victoria St (Melbourne City) -37.8076 144.9709
# 5   18011  8-Exhibition St/La Trobe St (Melbourne City) -37.8081 144.9690
# 6   18030    6-Swanston St/La Trobe St (Melbourne City) -37.8095 144.9641

计算

现在我们有了数据和距离公式，我们可以构造data.table连接

library(data.table)

## set the tram stop data as a data.table
dt1 <- as.data.table(tram_stops)

## add a column that will be used to do the join on
dt1[, joinKey := 1]

## find the dinstance between each point to every other point
## by joining the data to itself
dt2 <- dt1[
  dt1
  , {
    idx = dt.haversine(stop_lat, stop_lon, i.stop_lat, i.stop_lon) < 500 ## in metres
    .(stop_id = stop_id[idx],
      near_stop_id = i.stop_id)
  }
  , on = "joinKey"
  , by = .EACHI
]

结果

dt2 现在包含两列 stop_id，它们之间的距离在 500 米以内（包括与其自身相同的停靠点，因此可以将其移除）

dt2 <- dt2[stop_id != near_stop_id]

情节

由于我们使用googleway，让我们绘制一些结果（为此，您需要 Google Maps API 密钥，或使用其他映射库，例如传单）

mapKey <- "your_api_key"

## Just pick one to look at
myStop <- 18048
dt_stops <- dt3[stop_id == myStop ]

## get the lat/lons of each stop_id
dt_stops <- dt_stops[
  dt1      ## dt1 contains the lat/lons of all the stops
  , on = c(near_stop_id = "stop_id")
  , nomatch = 0
]

google_map(key = mapKey) %>%
  add_circles(data = dt1[stop_id == myStop], lat = "stop_lat", lon = "stop_lon", radius = 500) %>%
  add_markers(dt_stops, lat = "stop_lat", lon = "stop_lon")

注意事项

data.table 连接应该非常高效，但显然我在这里使用的数据只有 51 行；你必须让我知道这种方法对你的数据的扩展程度

【讨论】：