代码之家 › 专栏 › 技术社区 › dozyaustin

R(或Python)有没有加快地理空间分析的方法?

geospatial optimization for-loop r python

dozyaustin · 技术社区 · 7 年前

问题是这是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中。

2 回复 | 直到 7 年前

SymbolixAU Adam Erickson 7 年前

这里有一种方法 data.table ,以及我为其重新编写的哈弗森公式 this question 这样它就能在里面工作 数据桌子 操作

在每个点上连接到每个点,但在连接内计算每对点之间的距离,并删除阈值之外的点。这是受到@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 这是墨尔本城市环线有轨电车上的车站

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

计算

参加

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 <- dt2[stop_id != near_stop_id]

情节

googleway

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")

笔记

join应该相当有效,但显然我在这里使用的数据只有51行;你必须让我知道这个方法对你的数据的扩展性有多好

WombatPM 7 年前

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