0

我有纽约市 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

4

1 回答 1

0

尝试简单：

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 的转换，因为这不是这里的主要主题。

于 2015-02-24T14:53:10.630 回答