如何根據shapefile將緯度/經度數據分組到不同的組？

Question

原問題：我有一個數據集，每行有紐約範圍內的經度和緯度。現在我需要將每一行分組到紐約的一個郵編中。我有形狀文件與https://gis.ny.gov/gisdata/inventories/details.cfm?DSID=934可用的所有邊界。

添加樣本數據的緯度，經度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 
## .. ... ... 

由於您提供的座標的樣本，這裏是如何在閱讀它們，並將它們轉換爲您的紐約shape文件的投影，所以你可以做的聚合：用於快速回答

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 ::在