如何仅向外缓冲一组具有公共边的多边形（从边缘，而不是内部 "frontiers"）

Question

我正在使用 here 找到的爱尔兰共和国各县地图（直接下载）。

我想做的是将每个县的海岸向海延伸1500米，但不移动县与县之间的边界。

编辑：我还需要将沿海县的县界向大海延伸，以满足新的、延伸的“海岸线”（根据 mrhellmann 问题添加）

到目前为止我的代码：

加载包和数据

library(tidyverse)
library(sf)
library(ggplot2)

ireland <- st_read("~Census2011_Admin_Counties_generalised20m.shp")
st_crs(ireland) # units are m

缓冲区

ireland_buf <- ireland %>% 
  st_buffer(1500)

检查边界如何移动

ggplot() + 
  geom_sf(data = ireland_buf, col = "red", fill = "lightgray") + 
  geom_sf(data = ireland, col = "black", fill = NA) + 
  theme_bw()

所以，基本上，红线对应缓冲地图，黑线对应原始地图。我已经在海岸线上获得了预期的效果，但是县边界会根据县的不同向内或向外移动

Answer 1

按问题加载数据

library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
ireland_url <- "http://census.cso.ie/censusasp/saps/boundaries/Census2011_Admin_Counties_generalised20m.zip"
download.file(ireland_url, destfile = "ireland.zip")
f <- unzip("ireland.zip")
unlink("ireland.zip")
ireland <- st_read("Census2011_Admin_Counties_generalised20m.shp")
#> Reading layer `Census2011_Admin_Counties_generalised20m' from data source `/tmp/RtmpbGVDJH/reprexfc55361dad50/Census2011_Admin_Counties_generalised20m.shp' using driver `ESRI Shapefile'
#> Simple feature collection with 34 features and 20 fields
#> geometry type:  MULTIPOLYGON
#> dimension:      XY
#> bbox:           xmin: 17491.14 ymin: 19589.93 xmax: 334558.6 ymax: 466919.3
#> projected CRS:  TM65 / Irish National Grid
unlink(f)

使用 spatstat 说明的问题：

library(spatstat)
#> Loading required package: spatstat.data
#> Loading required package: spatstat.geom
#> spatstat.geom 1.65-8
#> Loading required package: spatstat.core
#> Loading required package: nlme
#> Loading required package: rpart
#> spatstat.core 1.65-9
#> Loading required package: spatstat.linnet
#> 
#> spatstat 2.0-0       (nickname: 'Caution: contains small parts') 
#> For an introduction to spatstat, type 'beginner'
originals <- lapply(st_geometry(ireland), as.owin)
dilated <- lapply(originals, dilation, r = 1500)
bb <- do.call(boundingbox, dilated)
plot(bb, main = "", type = "n")
for(i in seq_along(dilated)){
  plot(dilated[[i]], add = TRUE)
}

单个区域的可能解决方案（你必须自己制作一个额外的循环并根据需要转换回sf格式）

dilated_coast <- dilated[[34]]
for(i in seq_along(originals)[-34]){
  dilated_coast <- setminus.owin(dilated_coast, originals[[i]])
}
plot(dilated_coast, main = "", col = 2)
plot(originals[[34]], add = TRUE, col = 1)

Answer 2

使用基础 shapefile，您可以通过组合、缓冲、空间采样和创建 voronoi 多边形来创建一些其他对象。通过加入和交叉这些你可以得到大部分的方式。

以下示例的问题是它缓冲到北部边界的相邻陆地，一些靠得很近的岛屿可能会造成非常小的问题。

library(sf)
library(ggplot2)
library(dplyr)
library(tmap) # for "MAP_COLORS"


path <- 'Census2011_Admin_Counties_generalised20m.shp'
ireland <- read_sf(path)

base_crs <- st_crs(ireland)

#remove county lines, keeps only an outline of the area
ireland_combined <- st_combine(ireland)  
# add buffer, fortunately the base CRS is in meters
ireland_buffered <- st_buffer(ireland_combined, 1500)   
# polygon(s) of only the buffer
buffer_only <- st_difference(ireland_buffered, ireland_combined) 

# sample points inside the buffer for voronoi polygons
# type can be random, hexagonal, or regular.
# increase size for accuracy, decrease for speed.
sampled_points <- st_sample(buffer_only, size = 5000, type = 'hexagonal')

#Create voronoi polygons from sampled points
voronoi <- sampled_points %>%
  st_union() %>%
  st_voronoi() %>%
  st_cast()

# Crop voronoi polygons to buffer area only
# plot 2 below shows voronoi polygons extend much too far for this case
voronoi_buffer <- 
  st_intersection(
    st_make_valid(voronoi),
    buffer_only
  ) %>%
  st_as_sf()

# Join polygons by proximity to counties in base shapefile
voronoi_joined <- st_join(
  voronoi_buffer,
  ireland,
  join = st_nearest_feature
)

# join voronoi polygons by countyname
vj_summarized <- voronoi_joined %>% 
  group_by(COUNTYNAME) %>%
  summarize()

#mapview was having trouble with this county name
vj_summarized$COUNTYNAME[6] <- 'Dun Laoghaire-Rathdown'
ireland$COUNTYNAME[20] <- 'Dun Laoghaire-Rathdown'

#rbind sea buffer & land, then join with group_by & summarize
buffer_and_land <- rbind(ireland %>% select(COUNTYNAME), vj_summarized) %>%
  group_by(COUNTYNAME) %>%
  summarise()


#plots
p1 <- ggplot() + 
  geom_sf(data = ireland, fill = NA) + 
  geom_sf(data = ireland_buffered, fill = NA, alpha = .4, color = 'green') +
  ggtitle(label = "Ireland counties & 1500m buffer")

p2 <- ggplot() + 
  geom_sf(data = voronoi, alpha = .2) +
  geom_sf(data = voronoi_buffer, alpha = .2, fill = NA, color = 'turquoise') +
  ggtitle(label = 'Voronoi & buffer')

# Counties & buffers joined & colored
p4 <- tm_shape(st_geometry(buffer_and_land)) + 
  tm_polygons(col = "MAP_COLORS") + 
  tm_shape(st_geometry(ireland)) + 
  tm_polygons(border.col = 'black', col = NA, alpha = 0)

p1

绿色区域显示缓冲区和缓冲到边界土地的问题。您可以通过使用整个岛屿的 shapefile 来避免这种情况。

p2

已显示完整的 voronoi。绿松石感兴趣的领域。

p4

Buffer与县相连，以县名填充。县边界以黑色覆盖。

其中一个杂乱区域的特写，增加采样点的数量可能会有所帮助。 ^{由 reprex package (v1.0.0)}

创建于 2021-03-22

如何仅向外缓冲一组具有公共边的多边形（从边缘，而不是内部 "frontiers"）

How can I buffer a set of polygons with common edges only outwards (from the edges, not the inner "frontiers")

r

spatial

sf