我有纽约市 5449 棵树的经度和纬度,以及 55 个不同的邻里制表区 (NTA) 的 shapefile。每个 NTA 在 shapefile 中都有一个唯一的 NTACode,我需要在 long/lat 表中附加第三列,告诉我每棵树属于哪个 NTA(如果有)。
我已经在stackoverflow上使用其他多边形点线程取得了一些进展,尤其是这个查看多个多边形的线程,但是在尝试使用gContains时我仍然遇到错误并且不知道如何检查/标记不同多边形的每棵树(我猜是某种 sapply 或 for 循环?)。
下面是我的代码。数据/shapefile 可以在这里找到:http: //bit.ly/1BMJubM
library(rgdal)
library(rgeos)
library(ggplot2)
#import data
setwd("< path here >")
xy <- read.csv("lonlat.csv")
#import shapefile
map <- readOGR(dsn="CPI_Zones-NTA", layer="CPI_Zones-NTA", p4s="+init=epsg:25832")
map <- spTransform(map, CRS("+proj=longlat +datum=WGS84"))
#generate the polygons, though this doesn't seem to be generating all of the NTAs
nPolys <- sapply(map@polygons, function(x)length(x@Polygons))
region <- map[which(nPolys==max(nPolys)),]
plot(region, col="lightgreen")
#setting the region and points
region.df <- fortify(region)
points <- data.frame(long=xy$INTPTLON10,
lat =xy$INTPTLAT10,
id =c(1:5449),
stringsAsFactors=F)
#drawing the points / polygon overlay; currently only the points are appearing
ggplot(region.df, aes(x=long,y=lat,group=group))+
geom_polygon(fill="lightgreen")+
geom_path(colour="grey50")+
geom_point(data=points,aes(x=long,y=lat,group=NULL, color=id), size=1)+
xlim(-74.25, -73.7)+
ylim(40.5, 40.92)+
coord_fixed()
#this should check whether each tree falls into **any** of the NTAs, but I need it to specifically return **which** NTA
sapply(1:5449,function(i)
list(id=points[i,]$id, gContains(region,SpatialPoints(points[i,1:2],proj4string=CRS(proj4string(region))))))
#this is something I tried earlier to see if writing a new column using the over() function could work, but I ended up with a column of NAs
pts = SpatialPoints(xy)
nyc <- readShapeSpatial("< path to shapefile here >")
xy$nrow=over(pts,SpatialPolygons(nyc@polygons), returnlist=TRUE)
我们正在检查的 NTA 是这些(在 GIS 中可视化):http ://bit.ly/1A3jEcE