【问题标题】:Fastest way to determine COUNTRY from millions of GPS coordinates [R]从数百万个 GPS 坐标确定 COUNTRY 的最快方法 [R]
【发布时间】:2019-03-25 10:38:49
【问题描述】:

我有数百万个 GPS 坐标,想快速添加坐标国家的一列。

我目前的方法有效,但速度极慢:

library(data.table)

#REPRODUCE DATA
data <- data.table(latitude=sample(seq(47,52,by=0.001), 1000000, replace = TRUE),
                   longitude=sample(seq(8,23,by=0.001), 1000000, replace = TRUE))

#REQUIRED PACKAGES
if (!require("sp")) install.packages("sp")
if (!require("rworldmap")) install.packages("rworldmap")
if (!require("sf")) install.packages("sf")
library(sp)
library(rworldmap)
library(sf)

#CURRENT SLOW FUNCTION
coords2country = function(points,latcol,loncol){  
  countriesSP <- getMap(resolution='low')
  pointsSP <- st_as_sf(points,coords=c(loncol,latcol),crs=4326)
  pointsSP<- as(pointsSP,"Spatial")
  # use 'over' to get indices of the Polygons object containing each point 
  indices = over(pointsSP, countriesSP)
  # return the ADMIN names of each country
  indices$ADMIN  
  #indices$ISO3 # returns the ISO3 code 
  #indices$continent   # returns the continent (6 continent model)
  #indices$REGION   # returns the continent (7 continent model)
}

#SLOW!
> system.time(data[,country:=coords2country(data,"latitude","longitude"),])
   user  system elapsed 
121.293   7.849 130.226 

有没有更快/更好的方法来做到这一点?谢谢!

【问题讨论】:

  • 我想你见过this post。这篇文章的另一个选择是使用geonames 包。我想知道这是否适合你。
  • this question。我测试了map.where()。此功能运行速度非常快。
  • @jazzurro,感谢您的帮助! map.where() 大约快 10 倍。干杯!
  • 如果您的大部分点都在美国、加拿大、俄罗斯、澳大利亚...您可以手动找到仅包含给定国家/地区的最大矩形,并对那些进行预过滤,即除非提议的软件包已经进行了这种优化,否则可以大大加快速度。

标签: r coordinates geospatial sp sf


【解决方案1】:

有两个类似的问题。它们在我上面的 cmets 中。问题是询问如何从坐标中获取国家名称。在这里,OP 询问哪种方法可以更快地完成任务。

根据帖子,我们有三个选项。

  1. 在本题中使用自定义函数;
  2. 使用geonames 包;或
  3. map 包中使用map.where()

第二个选项需要一些设置。所以我刚刚测试了map.where()。以下是结果。正如 OP 所说,此功能的运行速度要快得多。

library(maps)
set.seed(111)
data <- data.table(latitude=sample(seq(47,52,by=0.001), 1000000, replace = TRUE),
                   longitude=sample(seq(8,23,by=0.001), 1000000, replace = TRUE))

system.time(data[, country := map.where(x = longitude, y = latitude)])

#   user  system elapsed 
#   7.20    0.05    7.29 

【讨论】:

  • 这确实很快,但有没有办法只获得主要国家而不是澳大利亚::塔兹马尼亚?
  • @HermanToothrot 对于附加问题:您可以使用strsplit(country, split = ":") 生成country:region-name 的列表,或者 - 或者 - 在特殊字符之后“替换”全部,即sub(pattern = ":.*", replacement = "")
猜你喜欢
  • 2018-10-31
  • 1970-01-01
  • 2023-03-09
  • 1970-01-01
  • 2023-03-22
  • 1970-01-01
  • 2016-08-09
  • 1970-01-01
  • 1970-01-01
相关资源
最近更新 更多