在 shapefile 中查找该点所属的多边形......为什么它适用于某些 shapefile，而不适用于其他 shapefile？答案

【问题标题】：Looking up polygons in a shapefile that point belong in... why does it work with some shapefiles, not with others?在 shapefile 中查找该点所属的多边形......为什么它适用于某些 shapefile，而不适用于其他 shapefile？
【发布时间】：2015-11-23 20:29:51
【问题描述】：

我正在尝试使用位于here 的 census.gov shapefile 按纬度/经度坐标查找人口普查区信息：

@lukeA 为 ZCTA shapefile 提出的代码已经非常成功。见原问题及解决here：

但是，当我将 shapefile 更改为人口普查区时，查找功能突然中断，我似乎无法弄清楚原因。 shapefile 上的documentation 表示它使用的是 NAD83 投影，与 ZCTA shapefile 相同。代码本身并没有返回任何错误......但输出都是NA。

代码如下：

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

mn.zip.map <- readShapePoly("tl_2015_27_tract.shp")
proj4string(mn.zip.map) <- CRS("+proj=utm +zone=15 +datum=NAD83")
mn.zip.map <- spTransform(mn.zip.map, CRS("+proj=longlat"))
latlong <- data.frame(matrix(0,nrow=3,ncol=1))
latlong$lat <- c(44.730178, 44.784711, 44.943008)
latlong$long <- c(-93.235381, -93.476415, -93.303201)
coordinates(latlong) = ~long+lat
proj4string(latlong) <- CRS(proj4string(mn.zip.map))
over(latlong, mn.zip.map)

有输出：

STATEFP COUNTYFP TRACTCE GEOID NAME NAMELSAD MTFCC FUNCSTAT ALAND AWATER INTPTLAT INTPTLON
1    <NA>     <NA>    <NA>  <NA> <NA>     <NA>  <NA>     <NA>    NA     NA     <NA>     <NA>
2    <NA>     <NA>    <NA>  <NA> <NA>     <NA>  <NA>     <NA>    NA     NA     <NA>     <NA>
3    <NA>     <NA>    <NA>  <NA> <NA>     <NA>  <NA>     <NA>    NA     NA     <NA>     <NA>

同样，我显然错过了 GIS 和地理编码的一些细微差别。

【问题讨论】：

标签： r gis shapefile

【解决方案1】：

您的预测有问题。使用rgdal::readOGR() 而不是maptools::readShapePoly()，这会自动处理投影信息：

library("rgdal")
mn.zip.map <- readOGR(".", "tl_2015_27_tract")
proj4string(mn.zip.map)
# [1] "+proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0"
latlong <- data.frame(lat=c(44.730178, 44.784711, 44.943008),
                      long=c(-93.235381, -93.476415, -93.303201))
coordinates(latlong) <- ~long+lat
proj4string(latlong) <- CRS(proj4string(mn.zip.map))
over(latlong, mn.zip.map)
#   STATEFP COUNTYFP TRACTCE       GEOID   NAME            NAMELSAD MTFCC
# 1      27      037  060812 27037060812 608.12 Census Tract 608.12 G5020
# 2      27      139  080301 27139080301 803.01 Census Tract 803.01 G5020
# 3      27      053  108000 27053108000   1080   Census Tract 1080 G5020
#   FUNCSTAT    ALAND  AWATER    INTPTLAT     INTPTLON
# 1        S  3654009  299711 +44.7246681 -093.2316029
# 2        S 29771331 2237086 +44.7886836 -093.4454977
# 3        S   678381       0 +44.9444245 -093.3014752

plot(mn.zip.map, axes=TRUE)
plot(latlong, add=TRUE, col="red")

【讨论】：

谢谢！这很好用。从现在开始，我肯定会使用 readOGR() 而不是 readShapePoly()。我对投影问题是什么有点好奇......但主要是很高兴有一个解决方案。再次感谢。
+proj=utm +zone=15 对 shapefile 错误，您在 (-97.49, 0.00040) 附近获得坐标；请参阅上面正确的 proj4string。