将 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")
获取这张地图,我们在其中看到该点位于美国边界内:
因此,我推荐的解决方案是使用这些分辨率更高、更可靠的边界多边形,特别是如果您有兴趣准确地解析靠近边界的点。
我正在尝试使用在以下 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")
获取这张地图,我们在其中看到该点位于美国边界内:
因此,我推荐的解决方案是使用这些分辨率更高、更可靠的边界多边形,特别是如果您有兴趣准确地解析靠近边界的点。