如何根据 shapefile 将纬度/经度数据分组到不同的组中？

Question

原始问题：我有一个数据集，其中每一行的纬度和经度都在纽约的范围内。现在我需要将每一行分组到纽约的一个邮政编码中。我有来自https://gis.ny.gov/gisdata/inventories/details.cfm?DSID=934 的所有可用边界的 shapefile。

添加经纬度样本数据http://pastebin.com/mXntxhK2

Answer 1

over / %over% 工作得很好，正如@RobertH 建议的那样。

library(sp)
library(raster)
library(rgdal)
library(dplyr)

# get the shapefile without wasting bandwidth
URL <- "http://gis.ny.gov/gisdata/data/ds_934/zip_codes_shp.zip"
fil <- "nyzips.zip"
if (!file.exists(fil)) download.file(URL, fil)
shp <- grep("shp$", unzip(fil), value=TRUE)
ny <- readOGR(shp[2], ogrListLayers(shp[2])[1], stringsAsFactors=FALSE)

# you didn't give us data so we have to create some by random sampling
# within the bounding box
ny_area <- as(extent(bbox(ny)), "SpatialPolygons")
set.seed(1492) # reproducible
pts <- spsample(ny_area, 3000, "random")
proj4string(pts) <- proj4string(ny)

# this does the lon/lat to zip mapping
zip_where <- pts %over% ny

# since we fabricated data, not all will be in a zip code since
# ny isn't a rectangle, so we remove the "bad" data
zip_where <- zip_where[complete.cases(zip_where),]

arrange(count(zip_where, POSTAL), desc(n))

## Source: local data frame [602 x 2]
## 
##    POSTAL     n
##     (chr) (int)
## 1   12847    16
## 2   12980    14
## 3   13367    14
## 4   13625    10
## 5   12843     9
## 6   12986     9
## 7   12134     8
## 8   12852     7
## 9   13324     7
## 10  13331     7
## ..    ...   ...

由于您提供了坐标样本，以下是如何读取它们并将它们转换为您的 NY shapefile 的投影，以便您进行聚合：

pts <- read.csv("http://pastebin.com/raw.php?i=mXntxhK2", na.strings="null")
pts <- pts[complete.cases(pts),]
coordinates(pts) <- ~longitude+latitude
proj4string(pts) <- CRS("+proj=longlat +datum=WGS84")
pts <- spTransform(pts, proj4string(ny))

# this does the lon/lat to zip mapping
zip_where <- pts %over% ny

# but since we fabricated data, not all will be in a zip code since
# ny isn't a rectangle, so we remove the "bad" data
zip_where <- zip_where[complete.cases(zip_where),]

arrange(count(zip_where, POSTAL), desc(n))
## Source: local data frame [158 x 2]
## 
##    POSTAL     n
##     (chr) (int)
## 1   11238    28
## 2   11208    25
## 3   11230    20
## 4   10027    19
## 5   11229    17
## 6   11219    16
## 7   11385    16
## 8   11206    15
## 9   11211    15
## 10  11214    14
## ..    ...   ...

Answer 2

这是一种方法：

library(raster)
library(rgeos)
# example data
filename <- system.file("external/lux.shp", package="raster")
zip <- shapefile(filename)
set.seed(0)
xy <- coordinates(spsample(zip, 10, 'random'))
plot(zip, col='gray')
points(xy, pch=20, col='red', cex=2)
# 
extract(zip, xy)

你也可以使用 sp::over