r - 将纬度/经度坐标转换为 R 中的位置时，NA 返回不正确

Question

我正在尝试使用以下链接中找到的 R 代码的修改版本:

Latitude Longitude Coordinates to State Code in R

为了测试代码，我创建了以下形式参数:

mapping = "state"
pointsDF = data.frame(x = c(-88.04607, -83.03579), y = c(42.06907, 42.32983))
latlong2state(pointsDF, mapping)

代码返回以下内容:

[1] "Illinois" NA

第一个坐标集返回正确答案，即“伊利诺伊州”。但是，当我将第二个坐标集(即 -83.03579、42.32983)输入在线转换器时，我得到以下结果:

Downtown, Detroit, MI, USA

( http://www.latlong.net/Show-Latitude-Longitude.html )

再次运行代码，但将第二个坐标从 42.32983 更改为 43.33 会将点置于密歇根州。

当使用“世界” map 作为“映射”变量的正式参数时，代码返回“美国”。我几天来一直在努力解决这个问题，但没有运气。我尝试过 SpatialPointDataFrames、各种投影，并研究了状态多边形对象本身。我在 Windows 7 系统上使用 R 版本 3.3.1。我认为有问题的数据点可能落在边界线上。在这种情况下，我认为“NA”是可以预料的。我使用的代码如下。

使用的代码:

library(sp)

library(maps)  
library(maptools)  
library(rgdal)

latlong2state = function(pointsDF, mapping) {

        local.map = map(database = mapping, fill = TRUE, col = "transparent", plot = FALSE) 
        IDs = sapply(strsplit(local.map$names, ":"), function(x) x[1])
        maps_sp = map2SpatialPolygons(map = local.map, ID = IDs, 
                                      proj4string = CRS("+proj=longlat +datum=WGS84"))                        
        pointsSP = SpatialPoints(pointsDF, 
                                 proj4string = CRS("+proj=longlat +datum=WGS84"))
        indices = over(x = pointsSP, y = maps_sp)
        mapNames = sapply(maps_sp@polygons, function(x) {x@ID})
        mapNames[indices]
}

我学习 R 才两个月，到目前为止我很喜欢这门语言。这是我第一次找不到答案。我非常感谢就此事提供的帮助!!!

Answer 1

首先，问题不是由位于边界上的点引起的。事实上，对于边界上的点，over() 不会返回 NA，而是“如果一个点落在多个多边形中，则记录最后一个多边形。”

NA表示不落在多边形内的点。我们可以放大您的 map 来查看情况

plot(local.map,  xlim = c(-83.2, -82.8), ylim=c(42.2,42.6), type="l")
polygon(local.map, col="grey60")
points(local.map)
points(pointsDF[2,], col="red")

根据 maps::map() 提供的多边形，该点位于加拿大毗邻的美国之外。正如你所说，当其他 map 将这一点定位在边境的美国一侧时，为什么会出现这种情况？我不认为这是投影问题，因为我们对多边形和点使用相同的 WGS84 地理坐标。因此，maps::map() 提供的多边形本身可能是错误的。

我们可以通过与其他来源的多边形进行比较来检查这一点。我从 http://www2.census.gov/geo/tiger/GENZ2015/shp/cb_2015_us_state_500k.zip 下载了美国人口普查部门最高分辨率的州边界。然后，

shp.path <- "C:/Users/xxx/Downloads/cb_2015_us_state_500k/cb_2015_us_state_500k.shp"
states <- readOGR(path.expand(shp.path), "cb_2015_us_state_500k")
plot(states,  xlim = c(-83.2, -82.8), ylim=c(42.2,42.6))
points(pointsDF[2,], col="red")

给我们提供了这张 map ，我们在其中看到该点位于美国边界之内:

因此，我推荐的解决方案是使用这些分辨率更高、更可靠的边界多边形，特别是如果您有兴趣准确解析靠近边界的点。