R中多边形内点（shapefile）的选择和提取

Posted 2023-02-25

【中文标题】R中多边形内点（shapefile）的选择和提取

【英文标题】：Selection and extraction of points (shapefile) within a polygon in R

【发布时间】：2020-06-02 07:51:58

【问题描述】：

我有两个形状文件；一个是点文件（一些世界信息），另一个是 21 个国家的形状文件。

我需要提取一个国家的分数。 我必须在 QGIS 或 ArcGIS 中重复此步骤 21 次。

有没有办法在 R 中做到这一点，如果是在批处理中，那真的很棒，因为我也必须对其他 4 个数据集重复这个。 在此先感谢

【问题讨论】：

也许检查over() 函数的映射overlay and spatial aggregation。

【参考方案1】：

所以作为一个例子如何使用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))

将结果直接写入数据框。

【参考方案2】：

这是一个示例多边形 shapefile 以及如何将其读入 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。它们都按点的顺序返回值，但 extract 还添加了（连续的）“点 id”——这主要用于当你有重叠的多边形时，这样一个点可能会落在两个多边形中。

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")