用具有最大共享边界的相邻多边形替换多边形
问题描述 投票：0回答：1
我有一个代表土地覆盖类型的多边形形状文件（其中包括“道路”类别），以及一个代表道路的线条形状文件。我在道路线周围创建了一个缓冲区，我想合并两个形状文件以将道路集成到土地覆盖类型的形状文件中。
我还希望删除土地覆盖形状文件中剩余的“道路”多边形，并将其替换为最接近并与道路多边形共享最多边界的土地覆盖类型的多边形。
这是我每个步骤的代码：
第 1 步：在道路线周围创建缓冲区
land_cover_types <- terra::vect("C:/R_scripts/land_cover_types.shp")
## terra::plot(land_cover_types)
roads <- terra::vect("C:/R_scripts/roads.shp")
## terra::plot(roads)
roads <- terra::buffer(roads, 20, quadsegs = 10, capstyle = "round", joinstyle = "round", mitrelimit = NA, singlesided = FALSE)
roads <- terra::aggregate(roads, by = "habitat_ty", dissolve = TRUE)
roads <- roads["habitat_ty"]
## terra::plot(roads)
第 2 步：合并两个 shapefile。 我的代码有问题。我期望将单个几何图形添加到土地覆盖形状文件中，并在其他土地覆盖类型下方添加“新道路”类别。
land_cover_types <- terra::erase(land_cover_types, roads)
land_cover_types <- terra::union(land_cover_types, roads)
## terra::plot(land_cover_types)
## terra::plot(roads, add = T, col = "red")
第3步：删除旧的道路多边形。我不知道如何做这部分。
例如，应从图像中删除道路多边形（带有蓝线和黑色星形），并替换为具有共享最多边界的黑色星形的土地覆盖多边形。土地覆盖多边形应包含道路多边形，类似于 ArcGIS 中的Eliminate 功能。
1个回答
0
投票
投票
事实证明这比预期的要长一些，但它仍然是一个部分答案，只是为了演示一个概念。但想法很简单 - 相邻多边形的交点是一组线（或点），通过为每个旧道路多边形选择最长的此类特征，我们可以选择目标多边形。
小型玩具数据集使调试和遵循变得更容易，
st_set_precision()lwgeom::st_snap_to_grid()library(sf)
#> Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is TRUE
library(dplyr, warn.conflicts = FALSE)
library(ggplot2)
# example data ------------------------------------------------------------
example_precision <- 100
old_road <-
st_as_sfc("LINESTRING(1 1, 2 2, 3 1, 4 2)") |>
st_buffer(.3) |>
st_set_precision(example_precision) |>
st_sf(geometry = _, class = "old")
other_cover_ty <-
st_make_grid(st_buffer(old_road, .2), n = c(3,1)) |>
st_set_precision(example_precision) |>
st_sf(geometry = _, class = LETTERS[1:3])
cover_ty <-
bind_rows(old_road, other_cover_ty) |>
st_difference() |>
st_set_precision(example_precision)
tibble(cover_ty)
#> # A tibble: 4 × 2
#> class geometry
#> <chr> <GEOMETRY>
#> 1 old POLYGON ((1.79 2.21, 1.8 2.22, 1.81 2.23, 1.82 2.24, 1.84 2.25, 1.85 2.…
#> 2 A POLYGON ((1.83 0.5, 0.5 0.5, 0.5 2.5, 1.83 2.5, 1.83 2.245, 1.82 2.24, …
#> 3 B MULTIPOLYGON (((3.17 1.59, 3 1.42, 2.21 2.21, 2.2 2.22, 2.19 2.23, 2.18…
#> 4 C POLYGON ((4.5 2.5, 4.5 0.5, 3.17 0.5, 3.17 0.755, 3.18 0.76, 3.19 0.77,…
new_road <-
st_as_sfc("LINESTRING(1 1.5, 4 1.5)") |>
st_buffer(.3) |>
lwgeom::st_snap_to_grid(.01) |>
st_sf(geometry = _, class = "new")
tibble(new_road)
#> # A tibble: 1 × 2
#> class geometry
#> <chr> <POLYGON>
#> 1 new ((4 1.8, 4.02 1.8, 4.03 1.8, 4.05 1.8, 4.06 1.79, 4.08 1.79, 4.09 1.79,…
# we can combine 2 sf objects and remove overlaps with bind_rows + st_difference
# (order matters(!); and it's probably not the most performant option)
combine_remove_overlaps <- function(x, y, add_id = TRUE){
result <-
# we want y on top as first shapes remain in case of an overlap
bind_rows(y, x) |>
st_difference() |>
st_cast("MULTIPOLYGON") |>
st_cast("POLYGON")
if (add_id) tibble::rowid_to_column(result) else result
}
cover_plot <- function(x){
ggplot(x) +
# add alpha to check for overlapping polygons
geom_sf(aes(fill = class), alpha = .7) +
geom_sf_label(data = ~st_point_on_surface(.x), aes(label = rowid), alpha = .7) +
scale_fill_viridis_d(drop = FALSE) +
theme_minimal()
}
cover_ty_all_road <-
combine_remove_overlaps(cover_ty, new_road) |>
# for consistent fill values when plotting
mutate(class = as.factor(class))
#> Warning in st_cast.sf(st_cast(st_difference(bind_rows(y, x)), "MULTIPOLYGON"),
#> : repeating attributes for all sub-geometries for which they may not be
#> constant
此时我们有一个小示例，其中新的道路缓冲区（类new）被添加到原始覆盖类型多边形（A，B，C，old），重叠区域被删除：
cover_ty_all_road
#> Simple feature collection with 9 features and 2 fields
#> Geometry type: POLYGON
#> Dimension: XY
#> Bounding box: xmin: 0.5 ymin: 0.5 xmax: 4.5 ymax: 2.5
#> CRS: NA
#> rowid class geometry
#> 1 1 new POLYGON ((4 1.8, 4.02 1.8, ...
#> 2 2 old POLYGON ((1.8 2.22, 1.81 2....
#> 3 3 old POLYGON ((3.79 2.21, 3.8 2....
#> 4 4 old POLYGON ((3.21 0.79, 3.2 0....
#> 5 5 old POLYGON ((1.21 0.79, 1.2 0....
#> 6 6 A POLYGON ((0.5 2.5, 1.83 2.5...
#> 7 7 B POLYGON ((2.2 2.22, 2.19 2....
#> 8 8 B POLYGON ((2.38 1.2, 2.79 0....
#> 9 9 C POLYGON ((3.17 0.755, 3.18 ...
cover_plot(cover_ty_all_road)
#> Warning: st_point_on_surface assumes attributes are constant over geometries
目标是通过使用最长的共享边来选择目标，将剩余的old（2,3,4,5）合并到相邻的A或B或C之一。
# Build a recoding table based on shared edge length
# class_col: classification column
# from: source class (old road)
# ignore: polygon class(es) to be ignored (new road)
recode_by_edge_lenght <- function(x, class_col, from, ignore, id_col = "rowid"){
`%notin%` <- Negate(`%in%`)
old_idx <- which(x[[class_col]] == from) # [1] 2 3 4 5
other_idx <- which(x[[class_col]] %notin% c(ignore, from)) # [1] 6 7 8 9
# rowid_recoded class_recoded
# "rowid.1" "class"
resulting_names <- c(paste0(id_col, ".1"), class_col)
names(resulting_names) <- c(paste0(id_col, "_recoded"), paste0(class_col, "_recoded"))
# intersections of touching polygons are shared edges or points,
# by measuring intersection lengths and picking longest by "from" column we
# can find polygons we should merge those to
st_intersection(
x[old_idx, id_col],
x[other_idx,]
) |>
# example with lengths before slicing by rowid would look like this:
# rowid rowid.1 class geometry shared_edge_len
# 2 2 6 A MULTILINESTRING ((1.8 2.22,... 0.6334343
# 5 5 6 A MULTILINESTRING ((1.21 0.79... 1.6159754
# 2.1 2 7 B MULTILINESTRING ((1.83 2.24... 1.0136534
# 4 4 8 B MULTILINESTRING ((3.17 0.75... 1.0136534
# 3 3 9 C MULTILINESTRING ((3.79 2.21... 1.6159754
# 4.1 4 9 C MULTILINESTRING ((3.21 0.79... 0.6334343
slice_max(order_by = st_length(geometry), by = all_of(id_col)) |>
# we only need a recoding table for non-spatial join, drop geometry
st_drop_geometry() |>
rename(all_of(resulting_names))
}
old_road_recode <- recode_by_edge_lenght(cover_ty_all_road, class_col = "class", from = "old", ignore = "new")
#> Warning: attribute variables are assumed to be spatially constant throughout
#> all geometries
old_road_recode
#> rowid rowid_recoded class_recoded
#> 2.1 2 7 B
#> 5 5 6 A
#> 4 4 8 B
#> 3 3 9 C
# add recoded rowid and class values through join, update with coalesce, union
union_by_recode <- function(x, recode_df, class_col){
id_col <- names(recode_df)[1]
# set one-to-one as a safeguard, throws error if there are multiple matches
left_join(x, recode_df, by = id_col, relationship = "one-to-one") |>
# coalesce across all (x_recoded, x) column pairs, result in x
mutate(across(ends_with("_recoded"), \(x) coalesce(x, get(sub("_recoded", "", cur_column()))), .names = '{sub("_recoded", "", .col)}')) |>
# mutate(across(ends_with("_recoded"), \(x) coalesce(x, get(str_split_i(cur_column(), "_", 1))), .names = '{str_split_i(.col, "_", 1)}')) |>
# rowid class rowid_recoded class_recoded geometry
# 1 1 new NA <NA> POLYGON ((4 1.8, 4.02 1.8, ...
# 2 7 B 7 B POLYGON ((1.8 2.22, 1.81 2....
# 3 9 C 9 C POLYGON ((3.79 2.21, 3.8 2....
# 4 8 B 8 B POLYGON ((3.21 0.79, 3.2 0....
# 5 6 A 6 A POLYGON ((1.21 0.79, 1.2 0....
# 6 6 A NA <NA> POLYGON ((0.5 2.5, 1.83 2.5...
# 7 7 B NA <NA> POLYGON ((2.2 2.22, 2.19 2....
# 8 8 B NA <NA> POLYGON ((2.38 1.2, 2.79 0....
# 9 9 C NA <NA> POLYGON ((3.17 0.75, 3.18 0...
summarise(across(geometry, st_union), .by = all_of(c(id_col, class_col)))
}
(cover_ty_upd <- union_by_recode(cover_ty_all_road, old_road_recode, class_col = "class"))
生成的sf对象：
#> Simple feature collection with 5 features and 2 fields
#> Geometry type: POLYGON
#> Dimension: XY
#> Bounding box: xmin: 0.5 ymin: 0.5 xmax: 4.5 ymax: 2.5
#> CRS: NA
#> rowid class geometry
#> 1 1 new POLYGON ((4 1.8, 4.02 1.8, ...
#> 2 7 B POLYGON ((1.81 2.23, 1.82 2...
#> 3 9 C POLYGON ((3.18 0.76, 3.19 0...
#> 4 8 B POLYGON ((3.17 0.5, 1.83 0....
#> 5 6 A POLYGON ((0.84 1.25, 0.85 1...
cover_plot(cover_ty_upd)
#> Warning: st_point_on_surface assumes attributes are constant over geometries
但是...这里的问题是它并不总是适用于真实数据，因此可能需要一些预处理、拓扑检查和修复和/或只是摆弄精度。
让我们像以前一样加载形状和处理。我试图不要过多地改变输入数据，但恐怕有必要从 Web-Mercator 进行转换，以便测量（缓冲区半径和后面的区域）才有意义。如果需要，请随意将其转换回来 EPSG:3857 和/或选择任何其他更合适的投影而不是 UTM。
library(mapview)
shp_url <- "https://www.dropbox.com/scl/fo/6245s22xgso3thyjn2qp6/AFL0cHneRZ5Fd_nBjrDnB74?rlkey=7co2m5ev6vw1eps2lwgcb8o9r&e=1&st=hmvn93o4&dl=1"
shp <- curl::curl_download(shp_url, "tmp.shp.zip")
st_layers(shp)
#> Driver: ESRI Shapefile
#> Available layers:
#> layer_name geometry_type features fields crs_name
#> 1 roads Line String 107 1 WGS 84 / Pseudo-Mercator
#> 2 land_cover_types Polygon 212 1 WGS 84 / Pseudo-Mercator
# load and transform from Web-Mercator to UTM zone 18N (for proper metric buffer);
# create buffer, union by habitat_ty
road_buff <-
read_sf(shp, layer = "roads") |>
st_transform(32618) |>
st_buffer(dist = 20, nQuadSegs = 10) |>
summarise(across(geometry, st_union), .by = habitat_ty)
road_buff
#> Simple feature collection with 1 feature and 1 field
#> Geometry type: MULTIPOLYGON
#> Dimension: XY
#> Bounding box: xmin: 654117.3 ymin: 5017527 xmax: 657038.1 ymax: 5021843
#> Projected CRS: WGS 84 / UTM zone 18N
#> # A tibble: 1 × 2
#> habitat_ty geometry
#> <chr> <MULTIPOLYGON [m]>
#> 1 New road (((654990.8 5017573, 654992.1 5017570, 654993 5017567, 654993.3 50…
# load, transform, combine with road_buff
land_cover <-
read_sf(shp, layer = "land_cover_types") |>
st_cast("POLYGON") |>
st_transform(32618) |>
combine_remove_overlaps(road_buff)
#> Warning in st_cast.sf(read_sf(shp, layer = "land_cover_types"), "POLYGON"):
#> repeating attributes for all sub-geometries for which they may not be constant
#> Warning in st_cast.sf(st_cast(st_difference(bind_rows(y, x)), "MULTIPOLYGON"),
#> : repeating attributes for all sub-geometries for which they may not be
#> constant
land_cover
#> Simple feature collection with 380 features and 2 fields
#> Geometry type: POLYGON
#> Dimension: XY
#> Bounding box: xmin: 654117.3 ymin: 5017527 xmax: 657048.2 ymax: 5021843
#> Projected CRS: WGS 84 / UTM zone 18N
#> # A tibble: 380 × 3
#> rowid habitat_ty geometry
#> <int> <chr> <POLYGON [m]>
#> 1 1 New road ((654990.8 5017573, 654992.1 5017570, 65499…
#> 2 2 New road ((655131.5 5017814, 655129.9 5017812, 65512…
#> 3 3 New road ((655198.2 5018359, 655196.8 5018356, 65519…
#> 4 4 New road ((654239.7 5021776, 654242.4 5021777, 65424…
#> 5 5 New road ((655018.6 5021793, 655020.9 5021791, 65502…
#> 6 6 Deciduous and mixed forest ((656081.5 5018400, 656080.2 5018397, 65607…
#> 7 7 Deciduous and mixed forest ((655852.4 5018191, 655849.6 5018190, 65584…
#> 8 8 Deciduous and mixed forest ((656186.4 5018372, 656185.8 5018363, 65618…
#> 9 9 Deciduous and mixed forest ((656012.1 5018227, 656017.6 5018231, 65602…
#> 10 10 Deciduous and mixed forest ((656176.6 5018579, 656185.7 5018584, 65618…
#> # ℹ 370 more rows
# check areas by cover types
land_cover |>
summarise(across(geometry, st_union), n_poly = n(), total_area = st_area(geometry), .by = habitat_ty) |>
st_drop_geometry()
#> # A tibble: 8 × 3
#> habitat_ty n_poly total_area
#> * <chr> <int> [m^2]
#> 1 New road 5 1459240.
#> 2 Deciduous and mixed forest 242 7544179.
#> 3 Water 4 75772.
#> 4 Agriculture 24 1042579.
#> 5 Urban 47 248348.
#> 6 Wetland 26 1432519.
#> 7 Road 26 14096.
#> 8 Coniferous forest 6 66056.
mapview(land_cover, zcol = "habitat_ty")
在尝试将Road多边形与其他多边形合并之前：
(old_road_recode <- recode_by_edge_lenght(land_cover, class_col = "habitat_ty", from = "Road", ignore = "New road"))
#> Warning: attribute variables are assumed to be spatially constant throughout
#> all geometries
#> # A tibble: 25 × 3
#> rowid rowid_recoded habitat_ty_recoded
#> * <int> <int> <chr>
#> 1 104 17 Deciduous and mixed forest
#> 2 93 38 Deciduous and mixed forest
#> 3 94 38 Deciduous and mixed forest
#> 4 95 38 Deciduous and mixed forest
#> 5 114 43 Deciduous and mixed forest
#> 6 117 282 Deciduous and mixed forest
#> 7 98 59 Deciduous and mixed forest
#> 8 99 59 Deciduous and mixed forest
#> 9 100 59 Deciduous and mixed forest
#> 10 101 59 Deciduous and mixed forest
#> # ℹ 15 more rows
# distribution of target types
old_road_recode |>
count(habitat_ty_recoded)
#> # A tibble: 2 × 2
#> habitat_ty_recoded n
#> <chr> <int>
#> 1 Coniferous forest 1
#> 2 Deciduous and mixed forest 24
land_cover_upd <-
union_by_recode(land_cover, recode_df = old_road_recode, class_col = "habitat_ty")
mapview(land_cover_upd, zcol = "habitat_ty")
合并后：
land_cover_upd |>
summarise(across(geometry, st_union), n_poly = n(), total_area = st_area(geometry), .by = habitat_ty) |>
st_drop_geometry()
#> # A tibble: 8 × 3
#> habitat_ty n_poly total_area
#> * <chr> <int> [m^2]
#> 1 New road 5 1459240.
#> 2 Deciduous and mixed forest 242 7558060.
#> 3 Water 4 75772.
#> 4 Agriculture 24 1042579.
#> 5 Urban 47 248348.
#> 6 Wetland 26 1432519.
#> 7 Road 1 26.2
#> 8 Coniferous forest 6 66245.
mapview(land_cover_upd, zcol = "habitat_ty", alpha = .2, alpha.regions = .2) +
mapview(filter(land_cover_upd, habitat_ty == "Road"), layer.name = "remaining Road", col = "red", col.regions = "red")
最新问题
- 如何将\Illuminate\Http\UploadedFile存储到公共磁盘上？
- 删除上一步中的字段
- 反转 FloatVector（1/每个元素）
- 在 Apache 服务器上为 React Router 配置 htaccess 文件
- 编译 C 项目时包含 pdcurses 时出错
- 以编程方式将小部件添加到桌面
- 为什么我在 ASP.NET Core 中向 Facebook 请求长期令牌时收到“验证码格式无效”？
- SpringBoot“调用了预验证入口点。拒绝访问”
- 将 YouTube 视频安全地嵌入网页的 div 中
- otel 配置期间 aws s3 导出器错误
- SQL Server XML 数据 - 如何插入（如果为空）或更新（如果不为空）
- 鸭嘴兽水滴文本更改
- NumericUpDown 与 UpdateEditText()
- 存储位置 - OAuth 2.0 中的访问令牌和刷新令牌
- 使用 zIndex 进行元素标高
- 如何让jquery select2动态禁用一个选项？
- 使用带有 webpack 模板的 electro-forge 的 puppeteer-extra-plugin-stealth 时找不到模块“is-plain-object”错误
- 在角度 @for 列表中对变化进行动画处理，而在每个项目更改时列表不会“跳转”
- 为什么 String.intern() 似乎对相同的字符串返回不同的实例？
- 通过原始日期减去两个日期之间的周数来查找日期
© www.soinside.com 2019 - 2024. All rights reserved.