使用 R 语言计算各省市区县河流密度——基于1:100万公众版基础地理信息数据

今天我们一起来学习下如何使用 R 语言计算各省市区县河流密度，其中河流水系数据来源于1:100万公众版基础地理信息数据库。

在附件中我为大家准备了这个数据：China basic geographic database 2019.gdb，关于这个数据的介绍，大家可以参考这个网页：https://www.webmap.cn/commres.do?method=result100W

首先加载所需的 R 包：

library(tidyverse)
library(sf)
library(ggspatial)

1:100万公众版基础地理信息数据库中的水系数据不仅仅包含了中国范围的，还包含了一些邻国的，因此在后面的数据处理中我们会需要提取中国范围的水系数据，这里就需要用到中国的范围矢量数据，可以从省级行政区划合并得到：

read_sf("2021行政区划/省.shp") %>%
  st_union() -> cn
nngeo::st_remove_holes(cn) -> cn

cn %>%
  write_rds("cn.rds")

由于数据不够精细，合并各省的数据无法完全拼接，会产生一些“小洞”，nngeo::st_remove_holes() 函数可以自动去除多边形内部的“小洞”。

合并好的效果是这样的：

library(leaflet)
library(leafem)
leaflet() %>%
  addTiles("http://map.geoq.cn/ArcGIS/rest/services/ChinaOnlineCommunity/MapServer/tile/{z}/{y}/{x}",
           tileOptions(tileSize = 256, minZoom = 3, maxZoom = 17),
           attribution = "微信公众号 RStata") %>%
  addFeatures(cn, weight = 0.5)

在 China basic geographic database 2019.gdb 数据中，水系数据对应的是 HYDL 图册：

read_rds("cn.rds") -> cn

# 线状水系：HYDL
st_read("China basic geographic database 2019.gdb",
        layer = "HYDL", as_tibble = T) %>%
  st_transform(4326) -> riverline

提取国内的：

riverline %>%
  st_intersection(cn) -> riverline

去除类型为 POINT 的河流：

# 去除 POINT
riverline %>%
  mutate(type = st_geometry_type(.)) %>%
  dplyr::filter(type != "POINT") -> riverline

# 保存
riverline %>%
  write_rds("riverline.rds")

我们可以绘制一幅图来展示下河流水系分布：

# 绘图展示（非常耗时）
library(ggplot2)
# 字体设置
library(showtext)
showtext_auto(enable = TRUE)
font_add("caoshu", regular = "钟齐李洤标准草书符号.ttf")

library(rcartocolor)
read_sf("九段线.geojson") -> jdx
read_sf("海岸线/海岸线.shp") -> hax

riverline %>%
  ggplot() +
  geom_sf(aes(size = SHAPE_Length, color = SHAPE_Length)) +
  geom_sf(data = hax, color = "lightblue", size = 0.1) +
  geom_sf(data = jdx, size = 0.5, color = "white") +
  scale_color_carto_c(palette = "Teal",
                      direction = -1,
                      guide = "none") +
  theme_minimal() +
  theme(panel.background = element_rect(color = NA,
                                        fill = "#000026"),
        plot.background = element_rect(color = "#00003a",
                                       fill = "#00003a"),
        panel.grid.major = element_blank(),
        axis.text.x = element_blank(),
        axis.text.y = element_blank()) +
  scale_size(range = c(0.1, 0.4), guide = "none") +
  coord_sf(crs = "+proj=lcc +lat_1=30 +lat_2=62 +lat_0=0 +lon_0=105 +x_0=0 +y_0=0 +ellps=krass +units=m +no_defs") +
  annotate("text", x = -2000000, y = 1800000,
           label = "中\n国\n河\n流\n水\n系\n图",
           family = "caoshu", hjust = 0,
           color = "white", size = 8) +
  annotate("text", x = -1500000, y = 1800000,
           label = "微\n信\n公\n众\n号\nRStata",
           family = "caoshu", hjust = 0,
           color = "white", size = 6) -> p1
# ggsave(plot = p1, "riverline.png", width = 10,
#        height = 10, device = png)
# knitr::plot_crop("riverline.png")
#
# ggsave(plot = p1, "riverline.pdf", width = 10,
#        height = 10, device = cairo_pdf)
# knitr::plot_crop("riverline.pdf")

下面我们就可以计算各省市区县的河流密度了，这里我们以城市的为例，省级和区县级别的计算方法类似：

read_sf("2021行政区划/市.shp") -> city

通常我们可以直接使用下面的代码提取各城市的河流：

# riverline %>%
#   st_intersection(city) -> riverlinecity

不过因为这个河流数据过大，所以上面的代码很难运行出来结果，所以我们还是分步运行：

lapply(seq(1, nrow(riverline), by = 10000), function(x){
  print(x)
  riverline %>%
    slice(x:(x + 9999)) %>%
    st_intersection(city)
}) %>%
  bind_rows() -> riverlinecity

# riverlinecity %>%
#   write_rds("riverlinecity.rds")

我这里是每次处理 10000 条河流，然后再把处理结果合并起来。

然后就可以分城市计算河流长度并合并了：

read_rds("riverlinecity.rds") -> riverlinecity
# 分城市计算河流长度
riverlinecity %>%
  select(-SHAPE_Length) %>%
  mutate(length = st_length(riverlinecity)) %>%
  st_drop_geometry() %>%
  group_by(省, 省代码, 省类型, 市, 市代码, 市类型) %>%
  summarise(length = sum(length, na.rm = T)) %>%
  mutate(length = as.numeric(length),
         length = length / 1000) %>%
  ungroup() -> riverlinecity2

这里河流密度我们使用 河流长度/区域面积 表示，单位是 km/km^2，所以还需要计算每个城市的面积：

# 计算河流密度
riverlinecity2 %>%
  left_join(city) %>%
  st_sf() %>%
  mutate(totalarea = st_area(.)) -> citydf

citydf %>%
  st_drop_geometry() %>%
  mutate(totalarea = units::set_units(totalarea, km2),
         totalarea = as.numeric(totalarea)) %>%
  mutate(density = length / totalarea) %>%
  rename(河流长度_km = length,
             总面积_km2 = totalarea,
             河流密度_km_per_km2 = density) %>%
  writexl::write_xlsx("中国各城市河流密度.xlsx")

这样我们就计算得到了各个城市的河流密度数据。

最后再绘图展示下：

# 各城市
readxl::read_xlsx("中国各城市河流密度.xlsx") -> cityline
cityline %>%
  left_join(city) %>%
  st_sf() -> cityline

cityline %>%
  ggplot() +
  geom_sf(aes(fill = 河流密度_km_per_km2), size = 0.01) +
  geom_sf(data = hax, color = "lightblue", size = 0.1) +
  geom_sf(data = jdx, size = 0.5, color = "white") +
  scale_fill_carto_c(palette = "Teal",
                     direction = -1,
                     guide = "none") +
  theme_minimal() +
  theme(panel.background = element_rect(color = NA,
                                        fill = "#000026"),
        plot.background = element_rect(color = "#00003a",
                                       fill = "#00003a"),
        panel.grid.major = element_blank(),
        axis.text.x = element_blank(),
        axis.text.y = element_blank()) +
  scale_size(range = c(0.1, 0.4), guide = "none") +
  coord_sf(crs = "+proj=aea +lat_0=0 +lon_0=105 +lat_1=25 +lat_2=47 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs") +
  annotate("text", x = -2000000, y = 1800000,
           label = "各\n市\n河\n流\n密\n度\n图",
           family = "caoshu", hjust = 0,
           color = "white", size = 8) +
  annotate("text", x = -1500000, y = 1800000,
           label = "微\n信\n公\n众\n号\nRStata",
           family = "caoshu", hjust = 0,
           color = "white", size = 6) +
  annotation_scale(
    width_hint = 0.2,
    text_family = "caoshu",
    bar_cols = c("#69889F", "white"),
    text_col = "#69889F"
  ) +
  annotation_north_arrow(
    location = "tr", which_north = "false",
    width = unit(1.6, "cm"),
    height = unit(2, "cm"),
    style = north_arrow_fancy_orienteering(
      text_family = "caoshu",
      fill = c("#69889F", "white"),
      text_col = "#69889F"
    )
  ) -> p3

# ggsave(plot = p3, "cityline.png", width = 10,
#        height = 10, device = png)
# knitr::plot_crop("cityline.png")

省份和区县的数据计算也类似，这里就不再赘述啦！

点击这里跳转到 RStata 短书平台获取附件：使用 R 语言计算各省市区县河流密度——基于1:100万公众版基础地理信息数据