作为我昨天发布的一个问题的后续,R 中是否有包/函数可以为我提供经度和纬度定义的点所在的国家(和大陆)?类似于在 MatLab 中所做的工作。
数据帧看起来像这样...
Point_Name Longitude Latitude
University of Arkansas 36.067832 -94.173655
Lehigh University 40.601458 -75.360063
Harvard University 42.379393 -71.115897
我想输出上面添加了国家和大陆列的数据帧。作为一个额外的奖励,在美国有一个美国各州的专栏(在美国以外有“其他”专栏)?
要获取大洲,您可以使用此答案中的worldmap
修改coords2country
函数的最后一行,以创建一个Coords2continental
功能,如下所示。选择您想要的是6大陆模式还是7大陆模式。我将考虑将此代码放入worlddmap
。
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
coords2continent = function(points)
{
countriesSP <- getMap(resolution='low')
#countriesSP <- getMap(resolution='high') #you could use high res map from rworldxtra if you were concerned about detail
# converting points to a SpatialPoints object
# 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)
#indices$continent # returns the continent (6 continent model)
indices$REGION # returns the continent (7 continent model)
#indices$ADMIN #returns country name
#indices$ISO3 # returns the ISO3 code
}
这是一个测试。
points = data.frame(lon=c(0, 90, -45, -100, 130), lat=c(52, 40, -10, 45, -30 ))
coords2continent(points)
#[1] Europe Asia South America North America Australia
coords2country(points)
#[1] United Kingdom China Brazil United States of America Australia
因此,这是在Google上使用反向地理编码API的替代方法。此代码部分基于此引用。
在 df
上方调用数据帧,
reverseGeoCode <- function(latlng) {
require("XML")
require("httr")
latlng <- as.numeric(latlng)
latlngStr <- gsub(' ','%20', paste(round(latlng,2), collapse=","))
url <- "http://maps.google.com"
path <- "/maps/api/geocode/xml"
query <- list(sensor="false",latlng=latlngStr)
response <- GET(url, path=path, query=query)
if (response$status !=200) {
print(paste("HTTP Error:",response$status),quote=F)
return(c(NA,NA))
}
xml <- xmlInternalTreeParse(content(response,type="text"))
status <- xmlValue(getNodeSet(xml,"//status")[[1]])
if (status != "OK"){
print(paste("Query Failed:",status),quote=F)
return(c(NA,NA))
}
xPath <- '//result[1]/address_component[type="country"]/long_name[1]'
country <- xmlValue(getNodeSet(xml,xPath)[[1]])
xPath <- '//result[1]/address_component[type="administrative_area_level_1"]/long_name[1]'
state <- xmlValue(getNodeSet(xml,xPath)[[1]])
return(c(state=state,country=country))
}
st.cntry <- t(apply(df,1,function(x)reverseGeoCode(x[2:3])))
result <- cbind(df,st.cntry)
result
# Point_Name Longitude Latitude state country
# 1 University of Arkansas 36.06783 -94.17365 Arkansas United States
# 2 Lehigh University 40.60146 -75.36006 Pennsylvania United States
# 3 Harvard University 42.37939 -71.11590 Massachusetts United States
在API定义中,“administrative_area_level_1”是国家以下的最高行政区域。在美国,这些是州。在其他国家,定义有所不同(例如,可能是省份)。
顺便说一句,我相当肯定你把纬度和经度弄反了。