我有一个包含两个几何列的sf数据集。它看起来是这样的:
> trip_geo
dstid sourceid dest_geom source_geom
1 1 1 MULTIPOLYGON (((-2.607985 5... MULTIPOLYGON (((-2.607985 5...
2 1 2 MULTIPOLYGON (((-2.607985 5... MULTIPOLYGON (((-2.57022 51...
3 1 3 MULTIPOLYGON (((-2.607985 5... MULTIPOLYGON (((-2.593213 5...
4 1 4 MULTIPOLYGON (((-2.607985 5... MULTIPOLYGON (((-2.608686 5...
5 1 5 MULTIPOLYGON (((-2.607985 5... MULTIPOLYGON (((-2.512852 5...
激活的几何体是dest_geom
。
每一行对应于邻居之间的一次旅行。对于每一次旅行,我想知道哪些邻居阻碍了旅行。也就是说,如果您在每一行的source_geom和dest_geom之间绘制一条直线,哪些几何图形会接触到这条直线?我想要获得行的所有感人的几何图形,然后合并它们。
我有另一个数据集,其中的几何图形对应于每个id:
> id_geo
dstid geometry
1 1 MULTIPOLYGON (((-2.607985 5...
2 2 MULTIPOLYGON (((-2.57022 51...
3 3 MULTIPOLYGON (((-2.593213 5...
4 4 MULTIPOLYGON (((-2.608686 5...
5 5 MULTIPOLYGON (((-2.512852 5...
我想第一步是在每次旅行的source_geom和dest_geom的质心之间定义一条线。然后,创建一个sf,其中geometry列包含与线接触的多边形的列表(我不知道在sf中是否可以在一列中包含多个几何图形)。然后,合并包含在同一行/列表中的几何图形。
我不认为我处理这个问题的方式是正确的,因为据我所知,一个人不能在sf的两个几何上执行操作,比如定义一条线。此外,我不知道如何将列表集成到dataframe/sf。
如果你能提出一种更现实的方法来解决这个问题,我将不胜感激。
发布于 2020-03-16 01:50:46
@Spacedman的响应完美地在感兴趣的多边形之间画出了线。下面是(1)绘制线条和(2)选择与每条线相交的多边形并合并选择的代码。
#Draw a line between each source and destination
library(sf)
trip_geo$join_geom = st_sfc(sapply(1:nrow(trip_geo),function(j)
+ st_cast(st_union(c(st_centroid(trip_geo$dest_geom[j]),
+ st_centroid(trip_geo$source_geom[j]))),"LINESTRING")), crs=4326)
#For each trip, select the polygons that intersect the line, then union them together
df=data.frame()
for (i in 1:nrow(trip_geo)){
id_geo$touch=apply(st_intersects(id_geo, trip_geo$join_geom[i]), 1, any)
touch=subset(id_geo, touch==T)
touch=st_union(touch)
df[i,1]=touch #append the unioned polygons to an empty dataframe
}
#Add the unioned polygons to the original dataset
sf=st_sf(df, crs=4326)
bound=cbind(trip_geo, sf)
计算需要一些时间(线条绘制大约需要1分钟,循环大约需要1分钟,其中包含3,500个多边形),我确信我的代码可以改进。
发布于 2020-03-08 13:14:39
在sf数据框中有两个几何图形没有什么问题,但当它与采用隐式几何图形的函数一起使用时,其中只有一个可以是几何图形,例如获取集合几何图形质心的st_centroid(foo)
。
对于其他几何图形,您可以处理该命名列,例如st_centroid(foo$source_geom)
。
对于具有两个多边形几何图形的数据框nc
,您可以通过以下方法计算两个质心之间的线:首先将点合并为多点,然后将它们投射到直线上。例如,对于第一行:
> st_cast(st_union(c(st_centroid(nc$source_geom[1]), st_centroid(nc$dest_geom[1]))),"LINESTRING")
Geometry set for 1 feature
geometry type: LINESTRING
dimension: XY
bbox: xmin: -81.49826 ymin: 36.42228 xmax: -77.41056 ymax: 36.4314
epsg (SRID): NA
proj4string: NA
LINESTRING (-81.49826 36.4314, -77.41056 36.42228)
您必须逐行执行此操作,否则最终将对整个几何向量进行操作。
完整的示例。执行library(spdep)
和example(poly2nb)
以获取nc.sids
。
首先将其减少到两列和5个随机行:
> nc = nc.sids[,c("NAME","FIPS")]
> nc = nc[sample(1:nrow(nc.sids),5),]
> nc
Simple feature collection with 5 features and 2 fields
geometry type: MULTIPOLYGON
dimension: XY
bbox: xmin: -83.73952 ymin: 34.36442 xmax: -78.16968 ymax: 36.54607
epsg (SRID): NA
proj4string: NA
NAME FIPS geometry
82 Cumberland 37051 MULTIPOLYGON (((-78.49929 3...
96 Bladen 37017 MULTIPOLYGON (((-78.2615 34...
13 Granville 37077 MULTIPOLYGON (((-78.74912 3...
78 Macon 37113 MULTIPOLYGON (((-83.10629 3...
14 Person 37145 MULTIPOLYGON (((-78.8068 36...
假设特征1转到特征2,2到3,等等。创建新的几何列:
> nc$dest_geom = nc$geometry[c(2,3,4,5,1)]
现在制作连接质心的线:
> nc$join_geom = st_sfc(sapply(1:nrow(nc),function(i)st_cast(st_union(c(st_centroid(nc$geometry[i]), st_centroid(nc$dest_geom[i]))),"LINESTRING")))
绘图:
> plot(nc$geom)
> plot(nc$join_geom,add=TRUE,lty=2)
https://stackoverflow.com/questions/60579832
复制