今天我们一起来看一下如何使用 R 语言判断每个城市边界是否和河流交汇。
首先加载所需的 R 包：
library(tidyverse)
library(sf)
|
读取二级河流矢量数据：
read_sf("china/二级河流.shp") %>%
st_transform(4326) -> rivers
plot(rivers[1])
|
读取市级行政区划数据：
read_sf("2021行政区划/市.shp") %>%
select(-contains("类型")) %>%
st_make_valid() -> city
|
提取城市边界：
st_cast(city, "MULTILINESTRING") -> cityline
|
然后就可以计算提取两者相交的部分了：
在当前坐标系下，会采用球面几何运算方法，效率很低：
转换成 aea 坐标系计算速度更快（平面几何的算法）：
mycrs <- "+proj=aea +lat_0=0 +lon_0=105 +lat_1=25 +lat_2=47 +x_0=0 +y_0=0 +ellps=krass +units=m +no_defs"
cityline %>%
st_transform(mycrs) %>%
st_intersection(
st_transform(rivers, mycrs)
) -> cityres1
cityres1
|
#> Simple feature collection with 1473 features and 17 fields
#> Geometry type: GEOMETRY
#> Dimension: XY
#> Bounding box: xmin: -2454645 ymin: 2245482 xmax: 2207103 ymax: 5921142
#> Projected CRS: +proj=aea +lat_0=0 +lon_0=105 +lat_1=25 +lat_2=47 +x_0=0 +y_0=0 +ellps=krass +units=m +no_defs
#> # A tibble: 1,473 × 18
#> 省 省代码 市 市代码 OBJECTID FNODE_ TNODE_ LPOLY_ RPOLY_ LENGTH
#> * <chr> <dbl> <chr> <dbl> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 内蒙古自治区 150000 巴彦淖尔市… 150800 35 136 136 143 148 0.03
#> 2 内蒙古自治区 150000 鄂尔多斯市… 150600 35 136 136 143 148 0.03
#> 3 内蒙古自治区 150000 巴彦淖尔市… 150800 39 140 140 143 153 0.028
#> 4 内蒙古自治区 150000 鄂尔多斯市… 150600 39 140 140 143 153 0.028
#> 5 内蒙古自治区 150000 巴彦淖尔市… 150800 47 148 148 161 143 0.038
#> 6 内蒙古自治区 150000 鄂尔多斯市… 150600 47 148 148 161 143 0.038
#> 7 内蒙古自治区 150000 巴彦淖尔市… 150800 48 149 149 143 156 0.075
#> 8 内蒙古自治区 150000 鄂尔多斯市… 150600 48 149 149 143 156 0.075
#> 9 内蒙古自治区 150000 巴彦淖尔市… 150800 53 153 153 143 163 0.029
#> 10 内蒙古自治区 150000 鄂尔多斯市… 150600 53 153 153 143 163 0.029
#> # ℹ 1,463 more rows
#> # ℹ 8 more variables: HYD2_4M_ <dbl>, HYD2_4M_ID <dbl>, GBCODE <int>,
#> # NAME <chr>, LEVEL_RIVE <int>, LEVEL_LAKE <int>, Shape_Leng <dbl>,
#> # geometry <MULTIPOINT [m]>
|
绘图观察结果：
library(mapview)
mapview(cityres1) +
mapview(rivers)
|
该计算结果包含的城市列表：
cityres1 %>%
select(1:4) %>%
st_drop_geometry() %>%
distinct() -> cityres1
|
考虑到由于数据不够精细的原因，城市边界数据和河流数据可能存在实际相交但是数据计算结果不相交的问题，可以给河流创建个半径例如为 1km 的缓冲区：
rivers %>%
st_buffer(dist = units::set_units(1, km)) -> riversbuffer
|
绘图观察：
再次计算：
cityline %>%
st_transform(mycrs) %>%
st_intersection(
st_transform(riversbuffer, mycrs)
) -> cityres2
cityres2 %>%
select(1:4) %>%
st_drop_geometry() %>%
distinct() -> cityres2
|
可以看出这个时候的结果多了两个城市：
cityres2 %>%
anti_join(cityres1) %>%
left_join(city) %>%
st_sf() -> temp
|
绘图观察这两个城市是否符合要求：
mapview(riversbuffer) +
mapview(temp)
|
可以看出这两个城市不算是以河流为边界，还是舍去，因此最终我们还是选择 cityres1 的结果：
cityres1 %>%
writexl::write_xlsx("边界与二级河流交汇的城市列表.xlsx")
|
点击这里跳转到 RStata 短书平台获取附件：使用 R 语言判断每个城市边界是否和河流交汇
评论