从包含多边形/区域和点（纬度、经度）的 shapefile 中，找出每个点属于哪个多边形/区域？在 R 中答案

【问题标题】：From a shapefile with polygons/areas, and points (lat,lon), figure out which polygon/area each point belongs to? In R从包含多边形/区域和点（纬度、经度）的 shapefile 中，找出每个点属于哪个多边形/区域？在 R 中
【发布时间】：2015-11-13 23:13:04
【问题描述】：

在给定一组点和一个 shapefile 的情况下，我正在尝试识别给定点属于哪个多边形（ZCTA...又称邮政编码模拟）。虽然有几个此类问题，但几乎所有问题似乎都将我引向 QGIS。虽然如果需要我会去学习另一个工具，但在 R 中有没有简单的方法来做到这一点？我在 R 环境中经验丰富……在 GIS 领域没有那么多经验。

我使用的 shapefile 位于此处： ftp://ftp.gisdata.mn.gov/pub/gdrs/data/pub/us_mn_state_mngeo/bdry_zip_code_tabulation_areas/shp_bdry_zip_code_tabulation_areas.zip

我的第一次尝试是将 shapefile 加载为 SpatialPolygonsDataFrame，将点加载为 SpatialPointsDataFrame，然后使用“over()”获取匹配的多边形的索引：

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

mn.zip.map <- readShapePoly("zip_code_tabulation_areas.shp")
# The shapefile is the one referenced in the link above

latlon <- data.frame(matrix(0,nrow=2,ncol=1))
latlon$lat <- c(44.730178, 44.784711)
latlon$lon <- c(-93.235381, -93.476415)
latlon[1] <- NULL
coordinates(latlon) = ~lon+lat
indices <- over(latlon, mn.zip.map)

有结果：

> indices
ZCTA5CE10 GEOID10 CLASSFP10 MTFCC10 FUNCSTAT10 ALAND10 AWATER10 INTPTLAT10 INTPTLON10
1      <NA>    <NA>      <NA>    <NA>       <NA>      NA       NA       <NA>       <NA>
2      <NA>    <NA>      <NA>    <NA>       <NA>      NA       NA       <NA>       <NA>
      Shape_Leng Shape_Area
1         NA         NA
2         NA         NA

我希望第一行输出 ZCTA5CE10 == 55124，第二行输出 ZCTA5CE10 == 55379。但是，显然这不会发生。

似乎坐标系没有对齐...但它们都应该是纬度/经度，对吧？

我在这里缺少什么？提前致谢。

【问题讨论】：

标签： r gis shapefile

【解决方案1】：

我认为你必须设置和调整投影：

library(rgdal)
proj4string(mn.zip.map) <- CRS("+proj=utm +zone=15 +datum=NAD83")
mn.zip.map <- spTransform(mn.zip.map, CRS("+proj=longlat"))
proj4string(latlon) <- CRS(proj4string(mn.zip.map))
over(latlon, mn.zip.map)
#   ZCTA5CE10 GEOID10 CLASSFP10 MTFCC10 FUNCSTAT10   ALAND10 AWATER10  INTPTLAT10   INTPTLON10 Shape_Leng Shape_Area
# 1     55124   55124        B5   G6350          S  43572536  1759018 +44.7394617 -093.1938424   27059.59   45295591
# 2     55379   55379        B5   G6350          S 152635134  6181840 +44.7539755 -093.5146083   86609.93  158696544

【讨论】：

谢谢@lukeA！只是为了完整性...有一个需要的库调用... library(rgdal)... 但在那之后它工作得很好。