问题描述
许多用户询问如何检查多边形是否相交,但是,这些问题的当前答案不适用于我的情况。
我有两个对象。一个是名为“farms”的“Large Spatialpolygons”。这个对象有几个多边形(2011 年的总数),每个多边形表示不同农场的限制(参见打印屏幕)。
我的第二个对象是一个名为“slope_RJ_100m”的“Large Spatialpolygons DataFrame”。该对象将一个大区域划分为几个正方形多边形,每个正方形的面积为 10.000m^2(总共 310000 个多边形)。
对于“slope_RJ_100m”中的每个正方形(多边形),我想知道它们是否与“农场”中的任何多边形相交。换句话说,我想知道“slope_RJ_100m”中的每个特定方块是否都有一个农场(即使只是农场的一部分)。 我期望结果是 310000 行和两个变量,一个表示slope_RJ_100m 中的多边形,另一个用TRUE 或FALSE 表示该多边形是否有农场。
我试过了:
inters = gIntersection(slope_RJ_100m,farms)
这会产生大约 1500 个多边形的输出。 我不知道如何使用它来知道我的 310000 个多边形中的哪个有一个农场。
inters = raster::intersect(slope_RJ_100m,farms)
输出有 29144 个多边形。与前一种情况一样,不确定如何使用它来知道广场是否有农场。
和
inters = st_intersects(slope_RJ_100m,farms)
UseMethod("st_intersects") 中的错误: 没有适用于“c('SpatialpolygonsDataFrame','Spatialpolygons','Spatial','SpatialpolygonsNULL','SpatialVector','SPNULL') 类的对象的“st_intersects”方法。
我确定我的问题是微不足道的,我事先道歉。
解决方法
这是一个带有 terra
的示例(要运行它,您需要 terra 1.1-17)
library(terra)
# polygons
p1 <- vect("POLYGON ((0 0,8 0,8 9,0 9,0 0))")
p2 <- vect("POLYGON ((5 6,15 6,15 15,5 15,5 6))")
p3 <- vect("POLYGON ((8 2,9 2,9 3,8 3,8 2))")
p4 <- vect("POLYGON ((2 6,3 6,3 8,2 8,2 6))")
p5 <- vect("POLYGON ((2 12,3 12,3 13,2 13,2 12))")
p6 <- vect("POLYGON ((10 4,12 4,12 7,11 7,11 6,10 6,10 4))")
p <- rbind(p1,p2,p3,p4,p5,p6)
plot(p,col=rainbow(6,alpha=.5))
lines(p,lwd=2)
text(p)
relate(rbind(p1,p2),rbind(p3,p6),"intersects")
# [,1] [,2] [,3] [,4]
#[1,] TRUE TRUE FALSE FALSE
#[2,] FALSE FALSE FALSE TRUE
使用您的 SpatialPolygons*
,您应该能够执行以下操作:
s <- vect(slope_RJ_100m)
f <- vect(farms)
x <- relate(s,f,"intersects")
同样,要使用 st_intersects
,您需要使用 sf
对象,而不是 Spatial*
对象。类似的东西
library(sf)
ss <- st_as_sf(slope_RJ_100m)
ff <- st_as_sf(farms)
inters <- st_intersects(ss,ff)
使用上面的示例数据
s1 < st_as_sf(rbind(p1,p2))
s2 <- st_as_sf(rbind(p3,p6))
st_intersects(s1,s2)
#Sparse geometry binary predicate list of length 2,where the predicate was `intersects'
# 1: 1,2
# 2: 4