如何做一个“完整的"与 R 包 sf 结合 [英] How to do a "full" union with the R package sf

问题描述

我尝试使用 sf::st_union 在三个多边形之间进行联合.在下图中显示了来自 ArcGIS叠加、联合、全部"的结果，我希望通过使用 R 中的 sf 包获得与输出"中的五个不同多边形类似的结果.

I try to do a union between three polygons using sf::st_union. In the figure below showing the result from ArcGIS "Overlay, Union, All" I wish to obtain a similar result as the five different polygons in 'OUTPUT' by using the sf package in R.

library(sf)
a1 <- st_polygon(list(rbind(c(0, 10), c(45, 10), c(45, 90), c(0, 90), c(0, 10))))
a2 <- st_polygon(list(rbind(c(45, 10), c(90,10), c(90, 90), c(45, 90), c(45, 10))))
b <- st_polygon(list(rbind(c(15, 5), c(75, 5), c(75, 50), c(15, 50), c(15, 5))))
a <- st_sf(c(st_sfc(a1), st_sfc(a2)))
b <- st_sf(st_sfc(b))
a$station <- c(1, 2)
b$type <- "A"
ab_union <- st_union(a, b)

在这个简单的例子中，生成的 sf 对象ab_union"将只包含两个多边形，而不是预期的五个.我可以通过使用sf包中的函数来获得上图五个对象的想要的结果吗?

In this simple example the resulting sf object 'ab_union' will only contain two polygons, not the expected five. Can I get the wanted result with five objects as in the figure above by using functions in the sf package?

推荐答案

我没有找到可以一步完成所有事情的功能，但这是解决您问题的一种方法:

I didn't find a function that make everything in one step, but this is a way to resolve your problem:

library(sf)
#> Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3
library(tidyverse)

a1 <- st_polygon(list(rbind(c(0, 10), c(45, 10), c(45, 90), c(0, 90), c(0, 10))))
a2 <- st_polygon(list(rbind(c(45, 10), c(90,10), c(90, 90), c(45, 90), c(45, 10))))
b1 <- st_polygon(list(rbind(c(15, 5), c(75, 5), c(75, 50), c(15, 50), c(15, 5))))

a <- st_sf(station=c(1, 2), geometry=st_sfc(a1, a2))
b <- st_sf(type="A",   geometry=st_sfc(b1))

st_agr(a) = "constant" #to avoid warnings, but see https://github.com/r-spatial/sf/issues/406
st_agr(b) = "constant"

#Operations
plot(st_geometry(st_union(a,b)))

op1 <- st_difference(a,st_union(b)) #notice the use of st_union()
plot(st_geometry(op1), border="red", add=TRUE)

op2 <- st_difference(b, st_union(a)) #notice the order of b and a and st_union()
plot(st_geometry(op2), border="green", add=TRUE)

op3 <- st_intersection(b, a) #notice the order of b and a
plot(st_geometry(op3), border="blue", add=TRUE)

union <- rbind(op1, op2, op3) #Error because op1 (op2) doesn't have the column "type" ("station")
#> Error in match.names(clabs, names(xi)): names do not match previous names

op11 <- dplyr::mutate(op1, type=NA)
op22 <- dplyr::mutate(op2, station=NA)

union <- rbind(op11, op22, op3)
(as.data.frame(union)) #The row names must be ordered.
#>     station type                       geometry
#> 1         1 <NA> POLYGON ((15 10, 0 10, 0 90...
#> 2         2 <NA> POLYGON ((45 50, 45 90, 90 ...
#> 3        NA    A POLYGON ((75 10, 75 5, 15 5...
#> 11        1    A POLYGON ((15 10, 15 50, 45 ...
#> 1.1       2    A POLYGON ((45 50, 75 50, 75 ...
plot(union)


#Other approach for avoid create the new columns would be:

union2 <- dplyr::bind_rows(op1, op2, op3) #But see discusion here: https://github.com/r-spatial/sf/issues/49
#> Warning in bind_rows_(x, .id): Vectorizing 'sfc_POLYGON' elements may not
#> preserve their attributes

#> Warning in bind_rows_(x, .id): Vectorizing 'sfc_POLYGON' elements may not
#> preserve their attributes

#> Warning in bind_rows_(x, .id): Vectorizing 'sfc_POLYGON' elements may not
#> preserve their attributes

^{由 reprex 包 (v0.2.1) 于 2019 年 4 月 6 日创建}

^{Created on 2019-04-06 by the reprex package (v0.2.1)}

我参考的讨论:
https://github.com/r-spatial/sf/issues/406

https://github.com/r-spatial/sf/issues/49

这篇关于如何做一个“完整的"与 R 包 sf 结合的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！