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

Question

我正在尝试使用位于 here:

的 census.gov shapefile 通过 lat/lon 坐标查找人口普查区信息

这与@lukeA 为 ZCTA shapefile 提出的代码一起工作非常成功。查看原始问题和解决方案 :

但是，当我将 shapefile 更改为人口普查区时，查找功能突然中断，我似乎无法弄清楚原因。 shapefile 上的 documentation 表示它使用的是 NAD83 投影，与 ZCTA shapefile 相同。该代码本身没有 return 任何错误...但输出都是 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 和地理编码的一些细微差别。

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

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

Looking up polygons in a shapefile that point belong in... why does it work with some shapefiles, not with others?

gis

r

shapefile