r - R中的纬度经度坐标到州代码

Question

有没有一种快速的方法可以将纬度和经度坐标转换为R中的州代码？我一直在使用邮政编码包作为查找表，但是当我查询大量纬度/经度值时它太慢了

如果不在 R 中，有什么方法可以使用谷歌地理编码器或任何其他类型的快速查询服务来做到这一点？

谢谢!

Answer 1

这里有两个选项，一种使用sf，一种使用sp包函数。 sf 是用于分析空间数据的更现代的(并且在 2020 年推荐)包，但如果它仍然有用，我将留下我最初的 2012 年答案，展示如何使用 来做到这一点sp相关函数。

<小时/>

方法 1(使用 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         

要获得非常准确的结果，您可以下载包含 GADM 的地理包-从 this page 维护美国的行政边界。然后，加载状态边界数据并像这样使用它们:

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        

<小时/>

方法 2(使用 sp):

这是一个函数，它采用 48 个州内经纬度的 data.frame，并为每个点返回它所在的州。

大部分函数只是简单地准备 sp 中的 over() 函数所需的 SpatialPoints 和 SpatialPolygons 对象code> 包，它负责计算点和多边形的“交集”:

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