将 Lat/Long 坐标转换为 R 中的位置时不正确的 NA return

Incorrect NA return when converting Lat/Long Coordinates to location in R

我正在尝试使用在以下 link 中找到的 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

第一个坐标集returns一个正确答案,即"Illinois"。但是,当我将第二个坐标集(即 -83.03579、42.32983)输入在线转换器时,我得到以下信息:

Downtown, Detroit, MI, USA

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

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

当使用 "world" 映射作为 "mapping" 变量的正式参数时,代码 returns "USA"。几天来我一直在努力解决这个问题并且没有运气。我玩过 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 语言只用了两个月,到目前为止我很喜欢这门语言。这是我第一次找不到答案。如果能就此事提供帮助,我将不胜感激!!!

首先,问题不是因为点在边界上。事实上,对于边界上的点,over() 不会 return NA,而是 "if a point falls in multiple polygons, the last polygon is recorded."

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")

获取这张地图,我们在其中看到该点位于美国边界内:

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