我有两个形状文件;一个是点文件(一些世界信息),另一个是21个国家的形状文件。
我需要提取属于一个国家的点,我必须在QGIS或ArcGIS中重复这个步骤21次。我必须在QGIS或ArcGIS中重复这个步骤21次。
是否有任何方法在R中做到这一点,如果是在批量处理,这将是非常伟大的,因为我必须重复这个其他4数据集。先谢谢你
因此,作为一个例子,如何使用 over()
考虑以下模拟数据。
library(raster)
library(sp)
#generate point distribution | with (x,y) coordinates
set.seed(42);pts <- SpatialPoints(data.frame(x = runif(30, 1, 5), y = runif(30, 1, 5)))
#create polygons
df<-data.frame(X = c(rep(1,5),rep(2,5),rep(3,5),rep(4,5),rep(5,5)),
Y = c(rep(seq(1,5,1),5)))
df$cell<-1:nrow(df) #polygon identifier (for later)
#make into spatial object (probably better way to do this)
coordinates(df) <- ~X+Y
rast <- raster(extent(df), ncol = 5, nrow = 5)
grid<-rasterize(df, rast, df$cell, FUN = max)
grid<-rasterToPolygons(grid) #create polygons
#plot to check
plot(grid); points(pts)
#and extract
pointsinpolygon = over(SpatialPolygons(grid@polygons), SpatialPoints(pts))
最后一点你也可以做为
df$pointsinpolygon <- over(SpatialPolygons(grid@polygons), SpatialPoints(pts))
来直接将结果写入数据框。
这里是一个多边形形状文件的例子,以及如何将其读入R中。
library(raster)
# example polygons filename
f <- system.file("external/lux.shp", package="raster")
pol <- shapefile(f)
我没有匹配的点文件,所以我为这个例子生成了一些点。
pts <- spsample(pol, 5, type="regular")
crs(pts) <- crs(pol)
要想知道这些点属于哪个区域,你可以使用 over
或 extract
.
over(pts, pol)
# ID_1 NAME_1 ID_2 NAME_2 AREA
#1 3 Luxembourg 9 Esch-sur-Alzette 251
#2 2 Grevenmacher 7 Remich 129
#3 3 Luxembourg 11 Mersch 233
#4 2 Grevenmacher 6 Echternach 188
#5 1 Diekirch 1 Clervaux 312
extract(pol, pts)
# point.ID poly.ID ID_1 NAME_1 ID_2 NAME_2 AREA
#1 1 10 3 Luxembourg 9 Esch-sur-Alzette 251
#2 2 7 2 Grevenmacher 7 Remich 129
#3 3 12 3 Luxembourg 11 Mersch 233
#4 4 6 2 Grevenmacher 6 Echternach 188
#5 5 1 1 Diekirch 1 Clervaux 312
如果您只是想要交点,比如说落在 "Grenmacher "中的点,就像您可能想对您的国家做的那样,您可以用
g <- pol[pol$NAME_1 == "Grevenmacher", ]
i <- intersect(pts, g)
i
#class : SpatialPointsDataFrame
#features : 2
#extent : 6.27111, 6.27111, 49.53585, 49.7887 (xmin, xmax, ymin, ymax)
#crs : +proj=longlat +datum=WGS84 +no_defs
#variables : 5
#names : ID_1, NAME_1, ID_2, NAME_2, AREA
#min values : 2, Grevenmacher, 6, Echternach, 129
#max values : 2, Grevenmacher, 7, Remich, 188
也许可以这样写回磁盘
shapefile(i, "test.shp")