在R中将纬度和经度坐标转换为国家名称

Question

在R中将纬度和经度坐标转换为国家名称

27

我有一组纬度和经度坐标，想要确定它们所在的国家。

我修改了一个关于将经纬度转换为美国州的答案（来源于这个问题），并编写了一个可行的函数。但是我遇到了一个问题，就是 worldHires 地图（来自 mapdata 包）极其过时，包含了许多已经不存在的国家，例如南斯拉夫和苏联。

我该如何修改此函数，以使用更现代的包，例如 rworldmap？到目前为止，我只让自己感到沮丧...

library(sp)
library(maps)
library(rgeos)
library(maptools)

# The single argument to this function, points, is a data.frame in which:
#   - column 1 contains the longitude in degrees
#   - column 2 contains the latitude in degrees
coords2country = function(points)
{
    # prepare a SpatialPolygons object with one poly per country
    countries = map('worldHires', fill=TRUE, col="transparent", plot=FALSE)
    names = sapply(strsplit(countries$names, ":"), function(x) x[1])

    # clean up polygons that are out of bounds
    filter = countries$x < -180 & !is.na(countries$x)
    countries$x[filter] = -180

    filter = countries$x > 180 & !is.na(countries$x)
    countries$x[filter] = 180

    countriesSP = map2SpatialPolygons(countries, IDs=ids, proj4string=CRS("+proj=longlat +datum=wgs84"))


    # convert our list of points to a SpatialPoints object
    pointsSP = SpatialPoints(points, proj4string=CRS("+proj=longlat +datum=wgs84"))


    # use 'over' to get indices of the Polygons object containing each point 
    indices = over(pointsSP, countriesSP)


    # Return the state names of the Polygons object containing each point
    myNames = sapply(countriesSP@polygons, function(x) x@ID)
    myNames[indices]
}

##
## this works... but it has obsolete countries in it
## 

# set up some points to test
points = data.frame(lon=c(0, 5, 10, 15, 20), lat=c(51.5, 50, 48.5, 47, 44.5))

# plot them on a map
map("worldHires", xlim=c(-10, 30), ylim=c(30, 60))
points(points$lon, points$lat, col="red")

# get a list of country names
coords2country(points)
# returns [1] "UK"         "Belgium"    "Germany"    "Austria"    "Yugoslavia"
# number 5 should probably be in Serbia...

- Iain S

地图包现已更新，新增了许多国家。不再包括苏联、南斯拉夫等国家。 - ZacharyST

2个回答

11

您可以使用我的geonames包从http://geonames.org/服务中查找：

> GNcountryCode(51.5,0)
$languages
[1] "en-GB,cy-GB,gd"

$distance
[1] "0"

$countryName
[1] "United Kingdom of Great Britain and Northern Ireland"

$countryCode
[1] "GB"

> GNcountryCode(44.5,20)
$languages
[1] "sr,hu,bs,rom"

$distance
[1] "0"

$countryName
[1] "Serbia"

$countryCode
[1] "RS"

从r-forge获取，因为我不确定是否已发布到CRAN：

https://r-forge.r-project.org/projects/geonames/

是的，它依赖于外部服务，但至少它知道共产主义发生了什么... :)

- Spacedman

似乎API要求每个调用中都有一个内置的用户名，但GNcountryCode函数不允许使用用户名。您会调整API调用吗？@Spacedman - Amit Kohli

对于大数据集来说，这可能不是理想的（本地）解决方案。 - undefined

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Andy · Accepted Answer

感谢您细心构思的问题。只需要做几行更改就能使用包含最新国家信息的rworldmap，如下所示。我不是CRS的专家，但我认为我对proj4string所做的更改没有任何影响。其他人可能会对此发表评论。

这对我有用，并给出了：

> coords2country(points)
[1] United Kingdom     Belgium            Germany            Austria           
[5] Republic of Serbia

祝一切顺利，安迪

library(sp)
library(rworldmap)

# The single argument to this function, points, is a data.frame in which:
#   - column 1 contains the longitude in degrees
#   - column 2 contains the latitude in degrees
coords2country = function(points)
{  
  countriesSP <- getMap(resolution='low')
  #countriesSP <- getMap(resolution='high') #you could use high res map from rworldxtra if you were concerned about detail

  # convert our list of points to a SpatialPoints object

  # pointsSP = SpatialPoints(points, proj4string=CRS(" +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"))

  #setting CRS directly to that from rworldmap
  pointsSP = SpatialPoints(points, proj4string=CRS(proj4string(countriesSP)))  


  # 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)
}