R中州代码的纬度经度坐标
Posted
技术标签:
【中文标题】R中州代码的纬度经度坐标【英文标题】：Latitude Longitude Coordinates to State Code in R 【发布时间】：2012-02-03 19:24:03 【问题描述】：有没有一种快速的方法可以将经纬度坐标转换为 R 中的州代码？我一直在使用 zipcode 包作为查找表，但是当我查询很多纬度/经度值时它太慢了
如果不在 R 中，有没有办法使用谷歌地理编码器或任何其他类型的快速查询服务来做到这一点？
谢谢！
【问题讨论】：
也可以在这里查看我的回答，使用ggmap::revgeocode：***.com/questions/46150851/…
【参考方案1】：
这里有两种选择，一种使用sf，一种使用sp包函数。 sf 是用于分析空间数据的更现代的（并且在 2020 年推荐）包，但如果它仍然有用，我将留下我在 2012 年的原始答案，展示如何使用 sp相关函数。
方法一（使用sf）：
library(sf)
library(spData)
## pointsDF: A data.frame whose first column contains longitudes and
## whose second column contains latitudes.
##
## states: An sf MULTIPOLYGON object with 50 states plus DC.
##
## name_col: Name of a column in `states` that supplies the states'
## names.
lonlat_to_state <- function(pointsDF,
states = spData::us_states,
name_col = "NAME")
## Convert points data.frame to an sf POINTS object
pts <- st_as_sf(pointsDF, coords = 1:2, crs = 4326)
## Transform spatial data to some planar coordinate system
## (e.g. Web Mercator) as required for geometric operations
states <- st_transform(states, crs = 3857)
pts <- st_transform(pts, crs = 3857)
## Find names of state (if any) intersected by each point
state_names <- states[[name_col]]
ii <- as.integer(st_intersects(pts, states))
state_names[ii]
## Test the function with points in Wisconsin, Oregon, and France
testPoints <- data.frame(x = c(-90, -120, 0), y = c(44, 44, 44))
lonlat_to_state(testPoints)
## [1] "Wisconsin" "Oregon" NA
如果您需要更高分辨率的状态边界，请使用sf::st_read() 或其他方式将您自己的矢量数据作为sf 对象读入。一个不错的选择是安装 rnaturalearth 包并使用它从 rnaturalearthhires 加载状态向量层。然后使用我们刚刚定义的lonlat_to_state()函数，如下所示：
library(rnaturalearth)
us_states_ne <- ne_states(country = "United States of America",
returnclass = "sf")
lonlat_to_state(testPoints, states = us_states_ne, name_col = "name")
## [1] "Wisconsin" "Oregon" NA
要获得非常准确的结果，您可以从this page 下载包含GADM 维护的美国行政边界的地理包。然后，加载状态边界数据并像这样使用它们：
USA_gadm <- st_read(dsn = "gadm36_USA.gpkg", layer = "gadm36_USA_1")
lonlat_to_state(testPoints, states = USA_gadm, name_col = "NAME_1")
## [1] "Wisconsin" "Oregon" NA
方法二（使用sp）：
这是一个函数，它采用较低 48 个州内的经纬度数据帧，并为每个点返回它所在的州。
大部分函数只准备sp 包中的over() 函数所需的SpatialPoints 和SpatialPolygons 对象，它真正完成了计算点和多边形的“交点”的繁重工作：
library(sp)
library(maps)
library(maptools)
# The single argument to this function, pointsDF, is a data.frame in which:
# - column 1 contains the longitude in degrees (negative in the US)
# - column 2 contains the latitude in degrees
lonlat_to_state_sp <- function(pointsDF)
# Prepare SpatialPolygons object with one SpatialPolygon
# per state (plus DC, minus HI & AK)
states <- map('state', fill=TRUE, col="transparent", plot=FALSE)
IDs <- sapply(strsplit(states$names, ":"), function(x) x[1])
states_sp <- map2SpatialPolygons(states, IDs=IDs,
proj4string=CRS("+proj=longlat +datum=WGS84"))
# Convert pointsDF to a SpatialPoints object
pointsSP <- SpatialPoints(pointsDF,
proj4string=CRS("+proj=longlat +datum=WGS84"))
# Use 'over' to get _indices_ of the Polygons object containing each point
indices <- over(pointsSP, states_sp)
# Return the state names of the Polygons object containing each point
stateNames <- sapply(states_sp@polygons, function(x) x@ID)
stateNames[indices]
# Test the function using points in Wisconsin and Oregon.
testPoints <- data.frame(x = c(-90, -120), y = c(44, 44))
lonlat_to_state_sp(testPoints)
[1] "wisconsin" "oregon" # IT WORKS
【讨论】：
我必须将 wgs84 更改为 WGS84 才能使此示例正常工作。 @AgustínIndaco 不快，因为在我的代码中，状态的多边形层是由 maps 包提供的，它没有相应的邮政编码边界层。如果你找到这样一个层，你当然可以调整我的代码来使用它。或者，您可能希望查看“反向地理编码”，例如 here。 我发现这个答案产生的结果对于某些应用程序来说缺乏足够的精度。例如，38.83226,-76.98946 编码为马里兰州，而不是哥伦比亚特区。 34.97982,-85.42203 被编码为田纳西州，而不是乔治亚州。如果您像我一样使用 15,000 个点，这种方法会产生很多不正确的结果（我正在使用的数据集中大约有 900 个，我估计）。不过，我不确定有什么更好的解决方案。
通过将“state”更改为“county”，这也适用于县。
@LaissezPasser 感谢您提及这一点。要获得更准确的结果，您可以使用我刚刚在 方法 1 下发布的代码以及该部分底部附近提到的 GADM-maintained 数据集。【参考方案2】：
几行 R 就可以搞定。
library(sp)
library(rgdal)
#lat and long
Lat <- 57.25
Lon <- -9.41
#make a data frame
coords <- as.data.frame(cbind(Lon,Lat))
#and into Spatial
points <- SpatialPoints(coords)
#SpatialPolygonDataFrame - I'm using a shapefile of UK counties
counties <- readOGR(".", "uk_counties")
#assume same proj as shapefile!
proj4string(points) <- proj4string(counties)
#get county polygon point is in
result <- as.character(over(points, counties)$County_Name)
【讨论】：
【参考方案3】：参见 sp 包中的?over。
您需要将状态边界设置为 SpatialPolygonsDataFrame。
【讨论】：
【参考方案4】：示例数据（多边形和点）
library(raster)
pols <- shapefile(system.file("external/lux.shp", package="raster"))
xy <- coordinates(p)
使用 raster::extract
extract(p, xy)
# point.ID poly.ID ID_1 NAME_1 ID_2 NAME_2 AREA
#1 1 1 1 Diekirch 1 Clervaux 312
#2 2 2 1 Diekirch 2 Diekirch 218
#3 3 3 1 Diekirch 3 Redange 259
#4 4 4 1 Diekirch 4 Vianden 76
#5 5 5 1 Diekirch 5 Wiltz 263
#6 6 6 2 Grevenmacher 6 Echternach 188
#7 7 7 2 Grevenmacher 7 Remich 129
#8 8 8 2 Grevenmacher 12 Grevenmacher 210
#9 9 9 3 Luxembourg 8 Capellen 185
#10 10 10 3 Luxembourg 9 Esch-sur-Alzette 251
#11 11 11 3 Luxembourg 10 Luxembourg 237
#12 12 12 3 Luxembourg 11 Mersch 233
【讨论】：
【参考方案5】：使用sf 非常简单：
library(maps)
library(sf)
## Get the states map, turn into sf object
US <- st_as_sf(map("state", plot = FALSE, fill = TRUE))
## Test the function using points in Wisconsin and Oregon
testPoints <- data.frame(x = c(-90, -120), y = c(44, 44))
# Make it a spatial dataframe, using the same coordinate system as the US spatial dataframe
testPoints <- st_as_sf(testPoints, coords = c("x", "y"), crs = st_crs(US))
#.. and perform a spatial join!
st_join(testPoints, US)
ID geometry
1 wisconsin POINT (-90 44)
2 oregon POINT (-120 44)
【讨论】：
以上是关于R中州代码的纬度经度坐标的主要内容，如果未能解决你的问题，请参考以下文章