根据包含点的多边形为点添加标签(NYC civic地理空间数据)

Question

我有纽约市5449棵树的经度和纬度，以及55个不同邻里制表区域(NTA)的shapefile。每个NTA在shapefile中都有一个惟一的NTACode，我需要在long/lat表中追加第三列，告诉我每个树属于哪个NTA (如果有的话)。

我已经在stackoverflow上使用其他多边形点线程取得了一些进展，特别是this one that looks at multiple polygons，但我在尝试使用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是这些(在地理信息系统中可视化的)：http://bit.ly/1A3jEcE

Answer 1

简单地试一下：

ShapeFile <- readShapeSpatial("Shapefile.shp")
points <- data.frame(long=xy$INTPTLON10, 
                         lat =xy$INTPTLAT10,
                         stringsAsFactors=F)
dimnames(points)[[1]] <- seq(1, length(xy$INTPTLON10), 1)
points <- SpatialPoints(points)

df <- over(points, ShapeFile)

我省略了shapefile的转换，因为这不是这里的主要主题。