根據哪個多邊形包含它（紐約市公民地理空間數據）標記一個點

Question

我有紐約市5449棵樹的經度和緯度，以及55個不同鄰居製表區（NTA）的shapefile。每個NTA在shapefile中都有一個唯一的NTACode，我需要在long/lat表中追加第三列，告訴我每棵樹屬於哪個NTA（如果有的話）。

我已經取得了一些進展，已經使用其他在多邊形中的點在多邊形線程上，尤其是this one that looks at multiple polygons，但我仍然在嘗試使用gContains時出現錯誤，不知道如何檢查/標記每棵樹對於不同的多邊形（我猜一些sapply或for循環？）。

以下是我的代碼。數據/ shape文件可以在這裏找到：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) 

我們檢查的NTAS是這些的（在GIS可視化）：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因爲這不是這裏的主要主題。