不正确的NA回报R中

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在密歇根州的地位。

当使用“世界”地图作为我的“映射”变量的形式参数时，代码返回“USA”。我一直在努力寻找解决方法，并没有运气。我玩过SpatialPointDataFrames，各种投影，并查看状态多边形对象本身。我在Windows 7系统上使用R版本3.3.1。我认为有问题的数据可能正在落在边界线上。在这种情况下，我认为会出现“不适用”。我使用的代码如下。

代码中使用：

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表示不落入多边形的点。我们可以放大地图上看到这种情况

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()提供的多边形。为什么会出现这种情况呢？正如你所说，其他地图在边界的美国一侧定位了这一点？我不认为这是一个投影问题，因为我们对多边形和点使用相同的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") 

得到我们这个图中，我们看到的一点是美国的边界内：

因此，我建议的解决方案，就是使用这些更高的分辨率，更可靠边界多边形，特别是如果你有兴趣准确地解决接近边界的点。