在 R 中使用 sf 和 ggplot2 绘制河流地图

气象学家

发布于 2022-04-18 20:26:23

2.6K0

发布于 2022-04-18 20:26:23

文章被收录于专栏：气象学家

翻译: Google翻译 作者: 米洛斯·波波维奇 原文链接： https://milospopovic.net/map-rivers-with-sf-and-ggplot2-in-r/

今年的3月22日是“世界水日”，它关注淡水的重要性，并引起人们对无法获得安全用水的22亿人的关注。今年的重点是地下水，这是为我们的泉水、河流、湖泊和湿地提供食物的宝贵来源。

在本教程中，我们将深入探讨重建此地图的具体细节。我们将使用全球河流分类（GloRiC).GloRiC对世界野生动物基金会HydroSHEDS的全球河流网络进行监督分类，以在全球层面创建河流覆盖类型。该数据集包括超过3500万公里的河流和溪流，分为超过8个单独的河流。令人印象深刻😲！

灵感来自这篇(https://blog.benthies.de/blog/mapping-streams-and-rivers-with-ggplot-sf/)很酷的博客文章，我决定写一个简短的教程，根据长期平均流量估计制作一个具有不同宽度的清晰河流地图。让我们开始吧！

Tips：关于HydroSHEDS，请查看即将发布的全球河流数据的最新公告以及新网站：https://www.hydrosheds.org/

获取河流数据

让我们导入几个库，好吗？这次我们将只使用3个库：httr通过GET函数检索数据;tidyverse 和 sf 用于空间分析和数据整理。

windowsFonts(georg = windowsFont('Georgia'))

libs <- c("httr", "tidyverse", "sf")

installed_libs <- libs %in% rownames(installed.packages())
if (any(installed_libs == F)) {
  install.packages(libs[!installed_libs])
}

invisible(lapply(libs, library, character.only = T))

作为热身，我们通过以下方式向GloRiC数据库编写并发送请求此链接（https://ln2.sync.com/dl/89159c590/am6vbeqz-46g6k5w6-syzujnfz-2rgjqrvt/view/default/3443504720013）.该 URL 指向 Sync.com 存储库，其中 shapefile 位于压缩文件夹中。下面，我们下载名为“eu_rivers.zip”的压缩文件夹，设置进度条（如果您不喜欢详细输出，请随时省略后者），然后解压缩文件夹。最后，我们列出包含下载文件名称的所有 shapefile。

# 1. GET RIVERS DATA
#---------

get_data <- function(url, res, filenames, riv_list) {

  url <- "https://m201.syncusercontent1.com/mfs-60:cb0f370cfc9aad110e47af0c14c6d6a0=============================/p/HydroRIVERS_v10_eu_shp.zip?allowdd=0&datakey=ge8BrffNb5ff90cF/3Dr8CGlfSRTS2g+I3poOJEOpd65Nwv0agBsI6yzqEFKLIDb/9iORBIKcrW+jONfWXU+Dm2WNnWn2KcP3R2WQMtiHUpSXuIg9e3gkGUd+nMSCIDRVvkSLQfFifZAeYjNKnNMSWPKn5mJIC/0RZ3kWGVVfDgZRYHLVgMMtxn+47XRsMbvz2WOE8R8J4n0RB7xVI5tlpmFR8Bv3nZB0HYc43cdyYSW5DqyKGji5CkoYTBNzE/jow5oBri2TOZWNxWVV9ET/9/udwd8KuOH++GoRzpgVJZ5WH1LgW7IBaHnDPTHdYENxRe/4JNlXanhauLfqOMNMQ&engine=ln-1.11.18&errurl=GLK4RA2KpWnRCSNSKOrDC8T6XxhjDIuUGGVvUo/zJAGGma48px5C67R8o6en+Qmegmceb/EJkpCEYm4SRBZDNMs+CgNFjFgEYdN4byFeygcAnvSAclpxktxsJc18rczzZtNuL4oyVr7SxjPeDdyHsveNgnibmLp1UWA45i5K12NRlbfK55BDYGOrwLu2RRJUm4AWUNeYpo4f8yqZ8zxQGUuLcKGx68bKP5I7ri8Lf8nkGsxRvT/axtcNmUiaZ706TofwP6mkzzaDZs+xin46q3P//GOuLornz7vN+q/d89qYsgpm5KsqyqEV+DePj2K6/vMDT4BLYWtxLeY1ESZgyg==&header1=Q29udGVudC1UeXBlOiBhcHBsaWNhdGlvbi96aXA&header2=Q29udGVudC1EaXNwb3NpdGlvbjogYXR0YWNobWVudDsgZmlsZW5hbWU9Ikh5ZHJvUklWRVJTX3YxMF9ldV9zaHAuemlwIjtmaWxlbmFtZSo9VVRGLTgnJ0h5ZHJvUklWRVJTX3YxMF9ldV9zaHAuemlwOw&ipaddress=1398165091&linkcachekey=89159c590&linkoid=51000013&mode=101&sharelink_id=3443504720013&timestamp=1649430246943&uagent=0c3c442ba67c66fa3227e09aeb8c763c79676604&signature=1cdeb688c90dc709eb9e652f6ccaf3beaa64197c&cachekey=60:cb0f370cfc9aad110e47af0c14c6d6a0============================="

  res <- GET(url,
             write_disk("eu_rivers.zip"),
             progress())
  unzip("eu_rivers.zip")
  filenames <- list.files("HydroRIVERS_v10_eu_shp", pattern="*.shp", full.names=T)

  riv_list <- lapply(filenames, st_read)

  return(riv_list)
}

我们通过调用前面的函数将欧洲河流 shapefile 读入 R 中，以获取要导入的文件列表。然后，我们将st_read应用于它并检索列表对象。由于我们想要 sf 对象，因此获取列表的第一个组件就足够了，我们的愿望将得到满足。

get_rivers <- function(filenames, list_riv, eu_riv) {

  filenames <- get_data()
  list_riv <- lapply(filenames, st_read)
  eu_riv <- list_riv[[1]] %>% 
  st_cast("MULTILINESTRING")

  return(eu_riv)
}

在下一步中，我们确保我们的线不被视为单个 LINESTRING 对象，而是被视为互连的 MULTILINESTRING 对象。下面是我们的河流对象在表格格式下的外观。

Simple feature collection with 938544 features and 14 fields
Geometry type: MULTILINESTRING
Dimension:     XY
Bounding box:  xmin: -24.47292 ymin: 12.60208 xmax: 69.51667 ymax: 81.78958
Geodetic CRS:  WGS 84
First 10 features:
   HYRIV_ID NEXT_DOWN MAIN_RIV LENGTH_KM DIST_DN_KM DIST_UP_KM CATCH_SKM
1  20000001         0 20000001      0.69        0.0        5.6     13.23
2  20000002         0 20000002      4.09        0.0        7.7     30.58
3  20000003         0 20000003      5.09        0.0        9.2     20.30
4  20000004         0 20000004      3.95        0.0        9.8      8.42
5  20000005         0 20000005      1.77        0.0       12.0      6.01
   UPLAND_SKM ENDORHEIC DIS_AV_CMS ORD_STRA ORD_CLAS ORD_FLOW  HYBAS_L12
1        13.2         0      0.175        1        1        7 2120062410
2        30.6         0      0.401        1        1        7 2120062410
3        20.3         0      0.262        1        1        7 2120062380
4         8.4         0      0.107        1        1        7 2120062380
5        33.7         0      0.420        2        1        7 2120062370
                         geometry
1  MULTILINESTRING ((59.2625 8...
2  MULTILINESTRING ((58.24167 ...
3  MULTILINESTRING ((57.50833 ...
4  MULTILINESTRING ((57.07917 ...
5  MULTILINESTRING ((56.55208 ...

创建河流宽度

欧洲的河流和集水区是一个由数百万条线路组成的错综复杂的网络。区分最突出的河流是值得的。我们可以根据 GloRiC 数据库中的大小类分配不同的宽度。

幸运的是，数据库的创建者已经将河流分类为有序类别。其中之一是ORD_FLOW，这是一种基于其长期平均流量的河流的对数大小类别。有8个这样的类（3-10个），按降序排列。因此，我们只需要根据这些类分配不同的宽度。我们在下面使用 mutate 来创建宽度，并使用case_when来分配宽度。后者在SQL用户中是众所周知的。在此上下文中，如果满足条件，它会根据宽度列分配一个值。

# 2. CREATE RIVER WIDTHS
#---------

get_river_widths <- function(eu_riv) {

  eu_riv <- get_rivers() %>% 
  mutate(width = as.numeric(ORD_FLOW),
         width = case_when(width == 3 ~ 1,
                           width == 4 ~ 0.8,
                           width == 5 ~ 0.6,
                           width == 6 ~ 0.4,
                           width == 7 ~ 0.2,
                           width == 8 ~ 0.2,
                           width == 9 ~ 0.1,
                           width == 10 ~ 0.1,
                           TRUE ~ 0)) %>%
  st_as_sf()

  eu_riv$geometry <- eu_riv$geometry %>%
    s2::s2_rebuild() %>%
    sf::st_as_sfc()

  return(eu_riv)
}

较新版本的 sf 不使用平坦地球模型。相反，包使用 s2 库中的球面几何运算符。在我们的例子中，这会破坏代码，因为某些河流线具有无效的球形几何图形。

一个快速的解决方法是通过sf：：sf_use_s2（FALSE）关闭此功能。理想情况下，我们希望使用无效的球面几何来修复要素，以便 s2 可以对其进行处理。这就是我们在上面的块中应用 s2：：s2_rebuild（）的原因。

边界框

在我们用ggplot2做魔术之前再走几步。我们的目标包括欧洲和中东，因此我们希望确保我们主要占领欧洲。我们通过制作一个边界框来做到这一点。让我们使用 WGS84 坐标定义边界框的参数。在本教程中，我们将使用世界等距圆柱投影来展平地图。因此，我们首先定义此投影，然后转换坐标。

# 3. MAKE BOUNDING BOX
#---------

get_bounding_box <- function(crsLONGLAT, bbox, new_prj, bb) {

  crsLONGLAT <- "+proj=longlat +datum=WGS84 +no_defs"

  bbox <- st_sfc(
  st_polygon(list(cbind(
    c(-10.5, 48.5, 48.5, -10.5, -10.5),
    c(35.000, 35.000, 69.5, 69.5, 35.000)
    ))),
  crs = crsLONGLAT)

  new_prj <- st_transform(bbox, crs = 4087)
  bb <- st_bbox(new_prj)

  return(bb)
}

绘制欧洲河流地图

好了，伙计们，我们准备绘制欧洲河流的地图了。我们首先绘制河流线，并根据定义的宽度根据类和宽度分配特定颜色。

由于我们的目标是将视野缩小到欧洲，因此我们使用coord_sf根据预定义的边界框设置纬度和经度限制。

我们将使用蓝色阴影来绘制我们的河流类。此外，我们将大小限制定义为从 0 到 0.3 的数值范围。我鼓励你玩这个范围，看看你会得到什么。最后，我们使用一系列 alpha 值来使较大的河流在地图上突出显示。

# 4. MAP
#---------

get_river_map <- function(eu_riv, bbox, p) {

  eu_riv <- get_rivers()
  bbox <- get_bounding_box()

  p <- 
    ggplot() +
    geom_sf(data=eu_riv, aes(color=factor(ORD_FLOW), size=width)) +
    coord_sf(crs = 4087,
      xlim = c(bbox["xmin"], bbox["xmax"]), 
      ylim = c(bbox["ymin"], bbox["ymax"])) +
    labs(y="", subtitle="",
         x = "",
         title="Rivers of Europe",
         caption="©2022 Milos Popovic https://milospopovic.net\nSource: ©World Wildlife Fund, Inc. (2006-2013)\n HydroSHEDS database http://www.hydrosheds.org") +
    scale_color_manual(
        name = "",
        values = c('#08306b', '#08519c', '#2171b5', '#4292c6', '#6baed6', '#9ecae1', '#c6dbef', '#deebf7')) +
    scale_size(range=c(0, .3)) +
    scale_alpha_manual(values=c("3" = 1, "4" = 1, "5" = .7, "6" = .6, "7" = .4, "8" = .3, "9" = .2, "10" = .1)) +
    theme_minimal() +
    theme(text = element_text(family = "georg"),
      panel.background = element_blank(), 
      legend.background = element_blank(),
      legend.position = "none",
      panel.border = element_blank(),
      panel.grid.minor = element_blank(),
      panel.grid.major = element_blank(),
      plot.title = element_text(size=40, color="#2171b5", hjust=0.5, vjust=0),
      plot.subtitle = element_text(size=14, color="#ac63a0", hjust=0.5, vjust=0),
      plot.caption = element_text(size=10, color="grey60", hjust=0.5, vjust=10),
      axis.title.x = element_text(size=10, color="grey20", hjust=0.5, vjust=-6),
      legend.text = element_text(size=9, color="grey20"),
      legend.title = element_text(size=10, color="grey20"),
      strip.text = element_text(size=12),
      plot.margin = unit(c(t=1, r=-2, b=-1, l=-2),"lines"),
      axis.title.y = element_blank(),
      axis.ticks = element_blank(),
      axis.text.x = element_blank(),
      axis.text.y = element_blank())

  return(p)
}


p1 <- get_river_map()
ggsave(filename="european_rivers_new.png", width=7, height=8.5, dpi = 600, device='png', p1)

啊，瞧！🤩

好吧，这就是女士们和先生们！在本教程中，您学习了如何导入河流空间文件以及如何在 R 中制作欧洲的炫酷河流地图。随时检查完整代码这里，克隆存储库并根据需要重现、重用和修改代码。

本教程为您从 GloRiC 数据库映射其他河流网络打开了大门。事实上，你可以稍微调整一下我的代码，制作非洲，美洲或亚洲的河流地图。例如这里（https://ln2.sync.com/dl/6c3c8b220/6ej6p3ch-wquydwwn-qnz87mr8-2cggh99e/view/default/3443375760013）是指向非洲形状文件的链接。试试吧，让我知道它是怎么回事！