(10)连接到GIS的桥梁

发表于 更新于 分类于 ， 阅读次数： Waline： 本文字数： 14k 阅读时长 ≈ 49 分钟

本文由SCY翻译自《Geocomputation with R》第十章

R 等具有交互式控制台的语言的一个特点——严格来说是一个读取-求值-打印循环（REPL）—— 是你与它们互动的方式。与其依赖于在屏幕的不同部分上指点和点击，你可以将命令键入控制台，并使用 Enter 键执行它们。使用像RStudio或VS Code这样的交互式开发环境时，一个常见且有效的工作流程是将代码键入源文件的源编辑器中，并使用像Ctrl+Enter这样的快捷方式来控制代码的交互式执行。

前提条件

  • 本章需要安装 QGIS、SAGA和GRASS，并附加以下包：
1
2
3
4
5
6
7
8
9
10
11
12
library(sf)
library(terra)
remotes::install_github("r-tmap/tmap")
#> 
#> ── R CMD build ─────────────────────────────────────────────────────────────────
#> * checking for file ‘/private/var/folders/x6/0mbmh6790tz0rrgmt0bqqs3m0000gn/T/RtmpjyIDsf/remotesc09c4ea71a88/r-tmap-tmap-1a2b4a7/DESCRIPTION’ ... OK
#> * preparing ‘tmap’:
#> * checking DESCRIPTION meta-information ... OK
#> * checking for LF line-endings in source and make files and shell scripts
#> * checking for empty or unneeded directories
#> * looking to see if a ‘data/datalist’ file should be added
#> * building ‘tmap_3.3-4.tar.gz’
1
2
3
4
5
6
remotes::install_github("r-spatial/qgisprocess")
library(qgisprocess)
library(Rsagacmd)
library(rgrass)
library(rstac)
library(gdalcubes)

引言

R 等具有交互式控制台的interpreted语言的一个特点——严格来说是一个读取-求值-打印循环（REPL）—— 是你与它们互动的方式。与其依赖于在屏幕的不同部分上指点和点击，你可以将命令键入控制台，并使用 Enter 键执行它们。使用像RStudio或VS Code这样的交互式开发环境时，一个常见且有效的工作流程是将代码键入源文件的源编辑器中，并使用像Ctrl+Enter这样的快捷方式来控制代码的交互式执行。

命令行界面（CLI）并不是R的独特属性：大多数早期的计算环境都依赖于命令行“shell”，直到计算机鼠标在20世纪90年代的发明和广泛采用后，图形用户界面（GUI）才变得普遍。例如，GRASS 是最长持续开发的开源GIS软件，最初依赖于其命令行界面，直到它获得GUI。大多数流行的GIS软件项目都是由GUI驱动的。你可以通过系统终端和嵌入式CLI与QGIS、SAGA、GRASS和gvSIG 进行交互，但它们的设计鼓励大多数人通过“点和击”与之交互。根据流行QGIS软件的创建者的说法：这样做的无意的后果是，大多数GIS用户错失了由CLI驱动和可编程方法带来的优势。

随着“现代”GIS 软件的出现，大多数人希望点击完成任务。这很好，但命令行为你提供了大量的灵活性和力量。很多时候，你可以在命令行上以GUI的一小部分时间来完成某些操作。

“CLI 与 GUI” 的辩论不必是敌对的：两种工作方式都有优点，取决于一系列因素，包括任务（如绘制新要素适合GUI）、所需的可重复性水平和用户的技能组合。GRASS是一个主要基于CLI但也有突出GUI的GIS软件的很好例子。同样，虽然R主要关注其CLI，但像RStudio这样的IDE为提高可访问性提供了 GUI。软件不能被整齐地分类为基于CLI或基于GUI。然而，交互式命令行界面在以下方面具有几个重要优势：

  • 自动化重复任务
  • 实现透明度和可重复性
  • 通过提供修改现有功能和实现新功能的工具来鼓励软件开发
  • 发展具有高需求的未来防护和高效编程技能
  • 提高触摸打字，这是数字时代的一项关键技能

另一方面，良好的GUI也有优势，包括：

  • “浅”学习曲线意味着地理数据可以在无需数小时学习新语言的情况下进行探索和可视化
  • 支持“数字化”（创建新矢量数据集），包括跟踪、捕捉和拓扑工具。（mapedit R包允许在从R打开的浏览器窗口中快速编辑少量空间特性，但不支持专业的大规模制图数字化。）
  • 通过地面控制点和正射校正实现地理参考
  • 支持立体测图（例如，LiDAR 和运动结构）

专用GIS软件项目的另一个优点是它们通过“GIS桥梁”提供对数百个“地理算法”的访问。这些增强R解决地理数据问题能力的计算配方桥梁是本章的主题。

📌命令行界面是一个通过键入和输入连续命令（命令行）与计算机程序互动的环境。
Linux中的bash和Windows中的PowerShell是众所周知的例子，它们允许用户控制操作系统的几乎任何部分。
像RStudio和VS Code这样的集成开发环境（IDE）提供代码自动补全和其他功能，以改善用户在开发代码时的体验。

R是可重复数据分析工作流与GIS之间建立桥梁的自然选择，因为它起源于一个接口语言。R（及其前身S）的一个关键特性是，它提供了访问其他语言（特别是 FORTRAN和C）中的统计算法的途径，但是来自一种功能强大的高级函数语言，并具有直观的REPL环境，而C和FORTRAN缺乏这样的环境。R继续这一传统，与众多语言接口，特别是 C++。

尽管R不是作为命令行GIS设计的，但其与专用GIS接口的能力赋予了它惊人的地理空间能力。借助GIS桥梁，R可以复制更多样化的工作流，通过从编程环境和一致的CLI控制它们，还能增加额外的可重复性、可扩展性和生产力优势。此外，R 在地理计算的某些领域超越了 GIS，包括交互式/动画地图制作和空间统计建模。

本章重点介绍了三个成熟的开源GIS产品的“桥梁”，在表中总结：

  • QGIS，通过qgisprocess包
  • SAGA，通过Rsagacmd包
  • GRASS，通过rgrass包

还有其他相关的桥梁，包括不再维护的R包RPyGeo，作为专有 GIS软件的接口。
还有一些重大进展，使开源GIS软件能够从QGIS（见 docs.qgis.org）和GRASS（见 grasswiki.osgeo.org）编写和执行R脚本。

Table: Comparison between three open-source GIS. Hybrid refers to the support of vector and raster operations.

GIS First release No. functions Support
QGIS 2002 >1000 hybrid
SAGA 2004 >600 hybrid
GRASS 1982 >500 hybrid

除了上面提到的三种R-GIS桥接方法外,本章还简要介绍了R语言与空间库、空间数据库以及基于云的地球观测数据处理的接口。

qgisprocess：连接到QGIS及更远的桥梁

是最受欢迎的开源GIS。QGIS提供了一个统一的接口来访问QGIS的本地地理算法、GDAL，以及当它们被安装时，来自提供商如GRASS和SAGA。自2020年夏季发布的3.14版本以来，QGIS随附了用于访问大量地理计算功能的qgis_process命令行工具。qgis_process提供了对标准QGIS安装中的300 多个地理算法的访问，以及通过插件对诸如GRASS和SAGA这样的外部提供商的1,000 多个算法的访问。

qgisprocess 包提供了从 R 访问qgis_process的能力。该包要求安装QGIS以及本章中使用的任何其他相关插件，如 GRASS 和SAGA，并且系统可以访问它们。有关安装说明，请参阅qgisprocess的文档。如果您已经安装了Docker，那么使用qgisprocess的快速方法是通过本项目开发的qgis镜像。假设您已经安装了Docker并且有足够的计算资源，您可以使用以下命令运行带有qgisprocess和相关插件的R会话（有关详细信息，请参阅 geocompx/docker 仓库）：

1
docker run -e DISABLE_AUTH=true -p 8786:8787 ghcr.io/geocompx/docker:qgis
1
library(qgisprocess)

此包会自动尝试检测QGIS的安装，并在找不到时发出警告。^[你可以使用qgis_configure()查看检测过程的详细信息。]当配置失败时有几个可能的解决方案：你可以设置options(qgisprocess.path = "path/to/your_qgis_process")，或者设置R_QGISPROCESS_PATH环境变量。上述方法也可用于你有多个 QGIS 安装并想决定使用哪一个时。更多细节，请参考qgisprocess的’入门’小插图。

接下来，我们可以找到在我们的计算机上可用的插件（意味着不同的软件）。

1
qgis_plugins()

这告诉我们，GRASS（grassprovider）和SAGA（processing_saga_nextgen）插件在系统上是可用的，但尚未启用。因为我们在本章后面需要它们，所以让我们启用它们。

1
2
qgis_enable_plugins(c("grassprovider", "processing_saga_nextgen"), 
                    quiet = TRUE)

请注意，除了在您的系统上安装 SAGA 之外，您还需要安装 QGIS Python 插件 Processing Saga NextGen。您可以在 QGIS 内部通过插件管理器进行安装，或者在 Linux 上至少可以借助Python包qgis-plugin-manager以编程方式进行。

qgis_providers()会列出软件的名称和可用地理算法的相应数量。

1
qgis_providers()

输出表确认我们可以通过QGIS接口使用QGIS地理算法（native、qgis、3d）和来自第三方提供商GDAL、SAGA和GRASS 的外部算法。

我们现在已经准备好使用QGIS及其伙伴在R内部进行一些地理计算了！让我们尝试两个示例案例研究。第一个展示如何将具有不同边界的两个多边形数据集统一。第二个专注于从表示为栅格的数字高程模型中获取新信息。

矢量数据

假设你有两个具有不同空间单位（例如，区域、行政单位）的多边形对象。我们的目标是将这两个对象合并为一个，包括所有边界线和相关属性。我们再次使用已经遇到的不一致的多边形。spData包中提供了这两个多边形数据集，我们希望对这两个数据集都使用地理坐标参考系统（CRS）。

1
2
3
data("incongruent", "aggregating_zones", package = "spData")
incongr_wgs = st_transform(incongruent, "EPSG:4326")
aggzone_wgs = st_transform(aggregating_zones, "EPSG:4326")
1
2
3
4
5
6
7
8
#> The legacy packages maptools, rgdal, and rgeos, underpinning the sp package,
#> which was just loaded, will retire in October 2023.
#> Please refer to R-spatial evolution reports for details, especially
#> https://r-spatial.org/r/2023/05/15/evolution4.html.
#> It may be desirable to make the sf package available;
#> package maintainers should consider adding sf to Suggests:.
#> The sp package is now running under evolution status 2
#>      (status 2 uses the sf package in place of rgdal)


Illustration of two areal units: incongruent (black lines) and aggregating zones (red borders).

要找到执行这项工作的算法，我们可以搜索 qgis_algorithms() 函数的输出。这个函数返回一个数据框，包含所有可用的提供者和它们包含的算法。因此，如果你看不到预期的提供者，可能是因为你还需要安装一些外部GIS软件。

1
qgis_algo = qgis_algorithms()

qgis_algo对象包含许多列，但通常我们只对algorithm列感兴趣，该列结合了提供商和算法名称的信息。假设函数的简短描述中包含单词"union"，我们可以运行以下代码来找到感兴趣的算法：

1
grep("union", qgis_algo$algorithm, value = TRUE)

上述列表中的一个算法，"native:union"，听起来很有希望。下一步是找出这个算法是做什么的，以及我们如何使用它。这是qgis_show_help()的作用，它返回算法的简短总结，包括算法的功能、参数和输出。^[我们还可以使用qgis_get_description()、qgis_get_argument_specs()和qgis_get_output_specss()独立提取其中的一些信息。] 这使得输出相当长。以下命令返回一个数据框，每一行代表"native:union"所需的一个参数，每一列包含与之关联的名称、描述、类型、默认值、可用值和可接受值：

1
2
3
alg = "native:union"
union_arguments = qgis_get_argument_specs(alg)
union_arguments

包含在union_arguments$name中的参数有INPUT、OVERLAY、OVERLAY_FIELDS_PREFIX和OUTPUT。union_arguments$acceptable_values包含每个参数可能的输入值的列表。许多函数需要代表矢量层路径的输入；对于这样的参数，qgisprocess函数接受sf对象。当期望的是“指向栅格层的路径”时，可以使用terra和stars包中的对象。这可能非常方便，但是当你仅仅为了将其提交给qgisprocess算法而读取空间数据时，我们建议提供磁盘上空间数据的路径。qgisprocess执行地理算法的第一件事是将存储在你的R会话中的空间数据导出到QGIS所知道的格式，例如.gpkg或.tif文件。这可能会增加算法的运行时间。

qgisprocess的主要功能是qgis_run_algorithm()，它将输入发送到QGIS并返回输出。它接受所使用的算法名称和帮助列表中显示的一组命名参数，并执行预期的计算。在我们的情况下，三个参数看起来很重要——INPUT、OVERLAY和OUTPUT。第一个参数INPUT是我们的主要矢量对象incongr_wgs，而第二个参数OVERLAY是aggzone_wgs。最后一个参数OUTPUT是输出文件名，如果没有提供，qgisprocess会自动选择并在tempdir()中创建。

1
2
3
4
union = qgis_run_algorithm(alg,
  INPUT = incongr_wgs, OVERLAY = aggzone_wgs
)
union

运行上述代码行将会把我们的两个输入对象保存到临时的.gpkg文件中，然后在它们上面运行所选算法，并返回一个临时的.gpkg文件作为输出。qgisprocess包将qgis_run_algorithm()的结果存储为一个列表，在本例中，包含指向输出文件的路径。我们可以用read_sf()将这个文件读回R（例如，union_sf = read_sf(union[[1]])），或者直接用st_as_sf()。

1
union_sf = st_as_sf(union)

请注意，QGIS的并集操作通过使用两个输入层的交集和对称差集来将两个输入层合并为一个层（顺便说一下，这也是在 GRASS 和 SAGA 中进行并集操作时的默认设置）。这与st_union(incongr_wgs, aggzone_wgs)（见练习）不相同！

结果union_sf是一个多边形，特征数量比两个输入对象都要多。然而，注意其中许多多边形较小，并不代表实际区域，而是由于我们两个数据集具有不同细节层次的结果。这些错误的副产品被称为割缝多边形。识别割缝的一种方法是找到相对非常小的区域的多边形，例如，在这里是25000平方米，然后将它们去除。让我们寻找一个合适的算法。

1
grep("clean", qgis_algo$algorithm, value = TRUE)

这次找到的算法，v.clean，不包括在QGIS中，而是在GRASS GIS中。GRASS GIS的v.clean是一个用于清理空间矢量数据拓扑的强大工具。重要的是，我们可以通过qgisprocess来使用它。

与前一步类似，我们应该从查看此算法的帮助开始。

1
qgis_show_help("grass7:v.clean")

我们在这里省略了输出，因为帮助文本相当长，并且包含许多参数。^[还要注意，这些参数与QGIS的参数不同，都是小写的。]这是因为v.clean是一个多功能工具——它可以清理不同类型的几何体并解决不同类型的拓扑问题。

对于这个示例，让我们仅关注几个参数，但我们鼓励你访问该算法的文档以了解更多关于v.clean的功能。

1
2
3
qgis_get_argument_specs("grass7:v.clean") |>
  select(name, description) |>
  slice_head(n = 4)

此算法的主要参数是input——我们的矢量对象。接下来，我们需要选择一个工具——清理方法。^[也可以选择多个工具，然后按顺序执行。]v.clean中存在大约十几个工具，可以用来删除重复的几何体、删除线之间的小角度或删除小面积等。在这个情况下，我们对最后一个工具rmarea感兴趣。其中一些工具，包括rmarea，期望一个附加参数threshold，其行为取决于所选工具。在我们的情况下，rmarea工具将删除所有小于或等于提供的threshold的面积。请注意，无论输入图层的坐标参考系统如何，阈值都必须以平方米为单位指定。

让我们运行这个算法，把输出转化为新的sf对象——clean_sf

1
2
3
4
5
clean = qgis_run_algorithm("grass7:v.clean",
  input = union_sf, 
  tool = "rmarea", threshold = 25000
)
clean_sf = st_as_sf(clean)

结果，图右侧面板，看起来是我们期待的结果——消除了细碎的面。


Sliver polygons colored in red (left panel). Cleaned polygons (right panel).

栅格数据

数字高程模型（DEMs）包含每个栅格单元的高程信息。它们被用于许多目的，包括卫星导航、水流模型、表面分析或可视化。在这里，我们有兴趣从DEM栅格中推导出可用作统计学习预测因子的新信息。例如，各种地形参数可能有助于滑坡的预测。

对于本节，我们将使用dem.tif——蒙贡研究区域的数字高程模型（从Land Process Distributed Active Archive Center下载，也可参见?dem.tif）。它的分辨率约为30乘30米，并使用投影的CRS。

1
2
3
library(qgisprocess)
library(terra)
dem = system.file("raster/dem.tif", package = "spDataLarge")

terra包的terrain()命令已经允许计算一些基本的地形特征，例如坡度、方位、TPI（地形位置指数）、TRI（地形崎岖指数）、粗糙度和流向。然而，GIS程序提供了许多的地形特征，其中一些在某些情境下可能更适合。例如，地形湿润指数（TWI）在研究水文和生物过程方面被发现有用。让我们使用"wetness"作为关键字搜索此索引的算法列表。

1
2
3
qgis_algo = qgis_algorithms()
grep("wetness", qgis_algo$algorithm, value = TRUE) |>
  head(2)

上述代码的输出表明，所需的算法存在于SAGA GIS软件中。^[TWI也可以使用r.topidx GRASS GIS功能来计算。]虽然SAGA是混合GIS，但其主要关注点一直是栅格处理，特别是数字高程模型（土壤性质、地形属性、气候参数）。因此，SAGA特别擅长快速处理大型（高分辨率）栅格数据集。

"sagang:sagawetnessindex"算法实际上是一个修改过的TWI，其结果是对位于山谷底部的单元格的更现实的土壤湿润潜力。

1
2
3
qgis_show_help("sagang:sagawetnessindex")
#> Saga wetness index (sagang:sagawetnessindex)
#> Description

在这里，我们坚持对所有参数使用默认值。因此，我们只需指定一个参数——输入的DEM。当然，在应用此算法时，您应确保默认值与您的研究目标相符。^[有关“sagang：sagawetnessindex”的其他参数的解释详见https://gis.stackexchange.com/a/323454/20955。]"sagang:sagawetnessindex"不是返回一个而是四个栅格——集水区域、集水坡度、修改后的集水区域和地形湿润指数。

在从QGIS中运行SAGA算法之前，我们将默认的栅格输出格式从.tif更改为SAGA的本地栅格格式.sdat。因此，从现在开始，我们自己不指定的所有输出栅格都将写入.sdat格式。根据您使用的软件版本（SAGA、GDAL），这可能不是必需的，但通常这会在尝试读取使用SAGA创建的输出栅格时减少不必的麻烦。

1
2
3
4
options(qgisprocess.tmp_raster_ext = ".sdat")
dem_wetness = qgis_run_algorithm("sagang:sagawetnessindex",
  DEM = dem
)

我们可以通过在qgis_as_terra()函数中提供输出名称来读取选定的输出。既然我们已经从 QGIS 内部完成了 SAGA 处理，我们将栅格输出格式更改回.tif。

1
2
3
dem_wetness_twi = qgis_as_terra(dem_wetness$TWI)
# plot(dem_wetness_twi)
options(qgisprocess.tmp_raster_ext = ".tif")

你可以在下图的左侧面板上看到TWI地图。地形湿润指数是无单位的：其低值代表不会积聚水的区域，而较高的值显示将以增加的水平积聚水的区域。

数字高程模型的信息也可以分类，例如，到geomorphons——代表地形形态的10个类别的地貌现象，例如斜坡、山脊或山谷。这些表型在许多研究中被使用，包括滑坡敏感性、生态系统服务、人类流动性和数字土壤制图。

geomorphons算法的原始实现是在GRASS GIS中创建的，我们可以在qgisprocess列表中找到它，名为"grass7:r.geomorphon"：

1
2
3
4
grep("geomorphon", qgis_algo$algorithm, value = TRUE)
#> [1] "grass7:r.geomorphon" "sagang:geomorphons" 
qgis_show_help("grass7:r.geomorphon")
# output not shown

计算geomorphons需要一个输入的DEM（数字高程模型，elevation），并且可以通过一组可选参数进行自定义。这些包括，search ——用于计算视线长度的值，以及-m ——一个标志，指定搜索值将以米为单位提供（而不是单元格数量）。关于附加参数的更多信息可以在原始论文和GRASS GIS文档中找到。

1
2
3
4
dem_geomorph = qgis_run_algorithm("grass7:r.geomorphon",
  elevation = dem,
  `-m` = TRUE, search = 120
)

我们的输出，dem_geomorph$forms，包含一个有10个类别的栅格文件——每一个都代表一个地形形态。我们可以使用qgis_as_terra()将其读入R中，然后进行可视化或在后续的计算中使用它。

1
dem_geomorph_terra = qgis_as_terra(dem_geomorph$forms)

有趣的是，一些geomorphons与TWI（地形湿润指数）值之间存在联系，如图所示。最大的TWI值主要出现在山谷和凹地中，而最低的值则如预期那样出现在山脊上。


Topographic wetness index (TWI, left panel) and geomorphons (right panel) derived for the Mongón study area.

SAGA GIS

自动地质科学分析系统（System for Automated Geoscientific Analyses，简称SAGA；表提供了通过命令行界面（在Windows下是saga_cmd.exe，在Linux下是saga_cmd）执行SAGA模块的可能性（参见SAGA模块的wiki页面）。此外，还有一个Python接口（SAGA Python API）。Rsagacmd使用前者从R内部运行SAGA。

在本节中，我们将使用Rsagacmd来划分秘鲁Mongón研究区域在2000年9月22日的归一化植被指数（NDVI）的相似值区域，方法是使用SAGA GIS中的种子区域生长算法。^[阅读章节了解如何从遥感图像计算NDVI的详细信息。]

1
ndvi = rast(system.file("raster/ndvi.tif", package = "spDataLarge"))

要开始使用Rsagacmd，我们需要运行saga_gis()函数。它有两个主要用途：

  • 它动态地^[这意味着可用的库将取决于安装的SAGA GIS版本。]创建一个新对象，该对象包含指向所有有效SAGA-GIS库和工具的链接
  • 它设置了一般的包选项，例如raster_backend（用于处理栅格数据的R包）、vector_backend（用于处理矢量数据的R包）和cores（用于处理的CPU核心的最大数量，默认值：全部）
1
2
library(Rsagacmd)
saga = saga_gis(raster_backend = "terra", vector_backend = "sf")

我们的saga对象包含指向所有可用SAGA工具的连接。它被组织为库列表（工具组），库内部有工具列表。我们可以使用$符号访问任何工具（记得使用TAB键进行自动补全）。

种子区域生长算法的工作分为两个主要步骤。首先，通过在指定大小的局部窗口中找到方差最小的单元格来生成初始单元格（“种子”）。然后，使用区域生长算法将种子的相邻像素合并以创建均匀区域。

1
sg = saga$imagery_segmentation$seed_generation

在上述示例中，我们首先指向了imagery_segmentation库，然后是其seed_generation工具。我们还将其分配给了sg对象，以便在下一步中不重新输入整个工具代码。^[你可以在https://saga-gis.sourceforge.io/saga_tool_doc/8.3.0/imagery_segmentation_2.html上了解更多关于该工具的信息。]如果我们只输入sg，我们将得到该工具的快速概述和一个带有其参数、描述和默认值的数据框。你还可以使用tidy(sg)仅提取参数表。seed_generation工具将栅格数据集作为其第一个参数（features）；可选参数包括指定初始多边形大小的band_width。

1
2
ndvi_seeds = sg(ndvi, band_width = 2)
# plot(ndvi_seeds$seed_grid)

我们的输出是三个对象的列表：variance ——局部方差的栅格图，seed_grid ——生成种子的栅格图，以及seed_points ——生成种子的空间矢量对象。

我们使用的第二个SAGA GIS工具是seeded_region_growing。^[你可以在https://saga-gis.sourceforge.io/saga_tool_doc/8.3.0/imagery_segmentation_3.html上了解更多关于该工具的信息。]seeded_region_growing工具需要两个输入：我们在前一步计算的seed_grid和ndvi栅格对象。此外，我们可以指定几个参数，例如normalize来标准化输入特征，neighbour（4或8-邻域）和method。最后一个参数可以设置为0或1（区域生长基于栅格单元格的值及其位置，或仅基于值）。

在这里，我们只会将method更改为1，这意味着我们的输出区域将仅基于它们的NDVI值的相似性创建。

1
2
3
srg = saga$imagery_segmentation$seeded_region_growing
ndvi_srg = srg(ndvi_seeds$seed_grid, ndvi, method = 1)
plot(ndvi_srg$segments)

工具返回了三个对象的列表：segments、similarity、table。similarity对象是一个栅格，显示种子与其他单元之间的相似性，而table是一个数据框，存储有关输入种子的信息。最后，ndvi_srg$segments是一个栅格，显示我们的结果区域。我们可以使用as.polygons()和st_as_sf()将其转换为多边形。

1
2
3
ndvi_segments = ndvi_srg$segments |>
  as.polygons() |>
  st_as_sf()


Normalized difference vegetation index (NDVI, left panel) and NDVI-based segments derived using the seeded region growing algorithm for the Mongón study area.

生成的多边形（段）代表具有相似值的区域。它们还可以使用各种技术（例如，k-均值聚类、区域化（例如，SKATER）或监督分类方法）进一步聚合成更大的多边形。您可以在练习中尝试这样做。

R还具有其他工具来实现用相似值创建多边形（所谓的段）的目标。它包括SegOptim包，该包允许运行几种图像分割算法，以及supercells，实现了用于处理地理空间数据的超像素算法SLIC。

GRASS GIS

美国陆军工程研究实验室（USA-CERL）于1982年至1995年创建了地理资源分析支持系统（GRASS）的核心。自1997年以来，学术界一直在继续这项工作。与SAGA类似，GRASS最初专注于栅格处理，但自GRASS 6.0以来，才添加了先进的矢量功能。

GRASS将输入数据存储在GRASS GIS数据库中。关于矢量数据，GRASS默认是一个拓扑GIS，即它只存储一次相邻要素的几何形状。SQLite是矢量属性管理的默认数据库驱动程序，属性与几何形状（即GRASS GIS数据库）通过键链接（GRASS GIS矢量管理）。

在使用GRASS之前，必须设置GRASS GIS数据库（也可在R内部设置），用户可能会发现这个过程开始时有点吓人。首先，GRASS数据库需要自己的目录，其中又包含一个位置（有关详细信息，请参见GRASS GIS数据库帮助页面grass.osgeo.org）。位置存储一个项目或一个区域的地理数据。在一个位置内，可以存在几个地图集，通常是指不同的用户或不同的任务。每个位置还有一个PERMANENT地图集——一个自动创建的强制性地图集。为了与项目的所有用户共享地理数据，数据库所有者可以将空间数据添加到PERMANENT地图集。此外，PERMANENT地图集存储栅格数据的投影、空间范围和默认分辨率。总而言之，GRASS GIS数据库可能包含许多位置（同一位置的所有数据具有相同的CRS），每个位置可以存储许多地图集（数据集组）。有关GRASS空间数据库系统的更多信息，请参阅GRASS GIS快速入门。要在R中快速使用GRASS，我们将使用link2GI包，不过，人们也可以逐步设置GRASS GIS数据库。参见GRASS within R了解如何操作。请注意，如果第一次使用GRASS，接下来的段落中的代码说明可能难以理解，但通过逐行运行代码并检查中间结果，背后的推理应该会变得更加清晰。

在此，我们通过地理信息科学中最有趣的问题之一——旅行推销员问题引入rgrass。假设一位旅行推销员想要拜访24个客户。此外，他想在家开始和结束他的旅程，总共有25个地点，同时以尽可能短的距离覆盖。解决这个问题有一个最佳解决方案；然而，对于现代计算机来说，检查所有可能的解决方案（大部分）是不可能的。在我们的案例中，可能的解决方案数量对应于(25 - 1)! / 2，即24的阶乘除以2（因为我们不区分前进或后退方向）。即使一次迭代可以在纳秒内完成，这仍相当于9837145年。幸运的是，有聪明、近乎最优的解决方案，可以在这难以想象的时间的微小部分内运行。GRASS GIS提供了其中的一个解决方案（有关详细信息，请参阅v.net.salesman）。在我们的用例中，我们想找到伦敦街道上前25个自行车站（而不是客户）之间的最短路径（我们简单地假设第一个自行车站对应于我们的旅行推销员的家）。

1
2
data("cycle_hire", package = "spData")
points = cycle_hire[1:25, ]

除了自行车租赁点数据外，我们还需要该区域的街道网络。我们可以使用osmdata包的帮助从OpenStreetMap下载它。为此，我们将街道网络的查询（在OSM语言中称为“高速公路”）限制为points的边界框，并将相应的数据附加为sf对象。osmdata_sf()返回一个包含几个空间对象（点、线、多边形等）的列表，但在这里，我们只保留与其相关ID的线对象。^[为了方便读者，可以使用data("london_streets", package = "spDataLarge")将london_streets附加到全局环境中。]

1
2
3
4
5
6
7
library(osmdata)
b_box = st_bbox(points)
london_streets = opq(b_box) |>
  add_osm_feature(key = "highway") |>
  osmdata_sf()
london_streets = london_streets[["osm_lines"]]
london_streets = select(london_streets, osm_id)

现在我们有了数据，可以开始并启动GRASS会话了。幸运的是，link2GI包中的linkGRASS()让我们只用一行代码就可以设置GRASS环境。您唯一需要提供的是一个确定空间数据库的投影和范围的空间对象。首先，linkGRASS() 会在您的计算机上找到所有的GRASS安装。由于我们将ver_select设置为TRUE，我们可以交互式地选择找到的GRASS安装之一。如果只有一个安装，linkGRASS()会自动选择它。其次，linkGRASS() 建立了与GRASS GIS的连接。

1
2
library(rgrass)
link2GI::linkGRASS(london_streets, ver_select = TRUE)

在我们能使用GRASS地理算法之前，我们需要将数据添加到GRASS的空间数据库中。幸运的是，方便的函数write_VECT()为我们做到了这一点。（对于栅格数据，请使用write_RAST()。）在我们的例子中，我们在只使用第一个属性列的情况下添加街道和自行车租赁点数据，并在GRASS中将它们命名为london_streets和points。

1
2
write_VECT(terra::vect(london_streets), vname = "london_streets")
write_VECT(terra::vect(points[, 1]), vname = "points")

rgrass包期望其输入和输出为terra对象。因此，我们需要使用vect()函数将我们的sf空间矢量转换为terra的SpatVector，以便能够使用write_VECT()。^[您可以通过阅读在R中不同空间类之间的转换博客文章和对象格式之间的强制转换小册子了解更多关于在R中转换空间类的信息]

现在，这两个数据集都存在于GRASS GIS数据库中。要执行我们的网络分析，我们需要一个拓扑清洁的街道网络。GRASS的"v.clean"负责去除重复项、小角度和悬挂等。在这里，我们在每个交叉点处断开线条，以确保后续的路由算法实际上可以在交叉口右转或左转，并将输出保存在名为streets_clean的GRASS对象中。

1
2
3
4
execGRASS(
  cmd = "v.clean", input = "london_streets", output = "streets_clean",
  tool = "break", flags = "overwrite"
)

📌要了解GRASS GIS模块的可能参数和标志，您可以使用help标志。
例如，尝试execGRASS("g.region", flags = "help")。

我们的一些自行车站点可能不会完全位于街道段上。然而，为了找到它们之间的最短路线，我们需要将它们连接到最近的街道段上。"v.net"的连接操作符就是这样做的。我们将其输出保存在streets_points_con中。

1
2
3
4
5
execGRASS(
  cmd = "v.net", input = "streets_clean", output = "streets_points_con",
  points = "points", operation = "connect", threshold = 0.001,
  flags = c("overwrite", "c")
)

得到的清洁数据集作为"v.net.salesman"算法的输入，该算法最终找到所有自行车租赁站之间的最短路线。其中一个参数是center_cats，它需要一个数字范围作为输入。这个范围代表应该计算最短路径的点。由于我们想计算所有自行车站的路线，所以我们将其设置为1-25。要访问旅行推销员算法的GRASS帮助页面，请运行execGRASS("g.manual", entry = "v.net.salesman")。

1
2
3
4
5
execGRASS(
  cmd = "v.net.salesman", input = "streets_points_con",
  output = "shortest_route", center_cats = paste0("1-", nrow(points)),
  flags = "overwrite"
)

要查看我们的结果，我们将结果读入R，将其转换为仅保留几何图形的sf对象，并借助mapview包进行可视化。

1
2
3
4
route = read_VECT("shortest_route") |>
  st_as_sf() |>
  st_geometry()
mapview::mapview(route) + points


Shortest route (blue line) between 24 cycle hire stations (blue dots) on the OSM street network of London.

在此过程中有一些重要的注意事项：

  • 我们本可以使用GRASS的空间数据库，从而允许更快的处理。这意味着我们仅在开始时导出了地理数据。然后我们创建了新对象，但只将最终结果导入回R。要找出当前可用的数据集，请运行execGRASS("g.list", type = "vector,raster", flags = "p")。
  • 我们也可以从R内部访问已经存在的GRASS空间数据库。在将数据导入R之前，您可能想要进行一些（空间）子集选择。对于矢量数据，请使用"v.select"和"v.extract"。"db.select"让您选择矢量图层的属性表的子集，而不返回相应的几何图形。
  • 您也可以从正在运行的GRASS会话中启动R。
  • 请参考出色的GRASS在线帮助或execGRASS("g.manual", flags = "i")了解有关每个可用GRASS地理算法的更多信息。

何时使用什么?

推荐单一的R-GIS界面是困难的，因为使用取决于个人偏好、手头的任务和您对不同GIS软件包的熟悉程度，这反过来可能取决于您的研究领域。正如之前提到的，SAGA特别擅长快速处理大型（高分辨率）栅格数据集，常被水文学家、气候学家和土壤科学家使用。另一方面，GRASS GIS是此处介绍的唯一支持拓扑空间数据库的GIS，这对于网络分析甚至模拟研究特别有用。与GRASS和SAGA-GIS相比，QGIS更加用户友好，尤其对于初次使用GIS的用户，可能是最受欢迎的开源GIS。因此，qgisprocess是大多数用例的合适选择。
其主要优势包括：

  • 对多个GIS的统一访问，因此提供了>1000个地理算法包括重复功能，例如，您可以使用QGIS-，SAGA-或GRASS-地理算法执行覆盖操作
  • 自动数据格式转换（SAGA使用.sdat网格文件，GRASS使用其自己的数据库格式，但QGIS将处理相应的转换）
  • 它自动将地理R对象传递给QGIS地理算法并返回到R
  • 支持命名参数和自动默认值检索的便利功能（受rgrass的启发）

当然，有些情况下您肯定应该使用其他R-GIS桥梁。虽然QGIS是唯一提供对几个GIS软件包统一界面的GIS，但它只提供对相应第三方地理算法的子集的访问。因此，要使用SAGA和GRASS的完整功能集，请坚持使用Rsagacmd和rgrass。此外，如果您想借助地理数据库运行模拟，请直接使用rgrass，因为qgisprocess每次调用都会启动新的GRASS会话。最后，如果您需要拓扑正确的数据和/或空间数据库管理功能，例如多用户访问，我们推荐使用GRASS。

请注意，还有一些GIS软件包具有脚本接口，但没有专门的R包可以访问它们：gvSig、OpenJump和Orfeo Toolbox。^[请注意，link2GI提供了与Orfeo Toolbox的部分集成，并且您还可以通过qgisprocess访问Orfeo Toolbox地理算法。还要注意，可以通过R包traudem从R中访问TauDEM。]

其他桥梁

本章到目前为止的重点是R与桌面GIS软件的接口。我们强调这些桥梁，因为专用的GIS软件是众所周知的，是理解地理数据的常见“入门方式”。它们还提供了许多地理算法的访问途径。

其他的“桥梁”包括指向空间库的接口、空间数据库以及网络映射服务。本节只提供了部分可能性的剖析。由于R的灵活性，通过其调用系统中的其他程序和与其他语言的集成（尤其是通过Rcpp和reticulate），许多其他的桥梁是可能的。目的不是全面覆盖，而是展示在本章开头引用的那句话中访问“灵活性和力量”的其他方式。

GDAL 桥梁

GDAL 是一个支持许多地理数据格式的底层库。GDAL 如此有效，以至于大多数GIS程序在后台使用GDAL 来导入和导出地理数据，而不是重新发明轮子并使用定制的读写代码。但是，GDAL 提供的不仅仅是数据I/O。它拥有用于矢量和栅格数据的地理处理工具，用于在线服务栅格数据的瓦片创建功能，以及对矢量数据的快速栅格化功能。由于GDAL是一个命令行工具，因此可以通过R中的system()命令访问其所有命令。

下面的代码段演示了这一功能：linkGDAL()在计算机上搜索工作的GDAL安装，并将可执行文件的位置添加到PATH变量中，从而允许调用GDAL（通常仅在Windows下需要）。

1
link2GI::linkGDAL()

现在我们可以使用system()函数调用GDAL的任何工具。例如，ogrinfo提供矢量数据集的元数据。在这里，我们将使用两个附加标志调用此工具：-al列出所有层的所有特征，-so仅获取摘要（而不是完整的几何列表）。

1
2
3
4
5
6
7
8
9
10
11
12
our_filepath = system.file("shapes/world.gpkg", package = "spData")
cmd = paste("ogrinfo -al -so", our_filepath)
system(cmd)
#> INFO: Open of `.../spData/shapes/world.gpkg'
#>       using driver `GPKG' successful.
#>
#> Layer name: world
#> Geometry: Multi Polygon
#> Feature Count: 177
#> Extent: (-180.000000, -89.900000) - (179.999990, 83.645130)
#> Layer SRS WKT:
#> ...

其他常用的GDAL工具包括：

  • gdalinfo：提供栅格数据集的元数据
  • gdal_translate：在不同的栅格文件格式之间转换
  • ogr2ogr：在不同的矢量文件格式之间转换
  • gdalwarp：重新投影、转换和剪裁栅格数据集
  • gdaltransform：转换坐标

访问https://gdal.org/programs/ 查看GDAL工具的完整列表，并阅读它们的帮助文件。

通过link2GI提供的对GDAL的“链接”，可以用作从R或系统CLI进行更高级GDAL工作的基础。TauDEM (http://hydrology.usu.edu/taudem) 和Orfeo Toolbox (https://www.orfeo-toolbox.org/) 是提供命令行界面的其他空间数据处理库/程序——上述示例展示了如何通过R通过系统命令行访问这些库。反过来，这可能是以新的R包形式为这些库创建适当接口的起点。

然而，在深入研究创建新桥梁的项目之前，重要的是要意识到现有R包的功能，以及system()调用可能不是平台无关的（可能在某些计算机上失败）。此外，sf通过Rcpp提供的R/C++接口，将GDAL、GEOS和PROJ提供的大部分功能带到R中，从而避免了system()调用。

连接到空间数据库的桥梁

空间数据库管理系统（空间DBMS）以结构化的方式存储空间和非空间数据。它们可以通过唯一标识符（主键和外键）以及空间来组织大量的数据集合到相关表（实体）中（例如，考虑空间连接）。这很有用，因为地理数据集很容易变得庞大和混乱。数据库能够基于空间和非空间字段有效地存储和查询大型数据集，并提供多用户访问和拓扑支持。

最重要的开源空间数据库是PostGIS。^[SQLite/SpatiaLite肯定也很重要，但实际上我们已经介绍了这种方法，因为GRASS在后台使用SQLite。]
R到空间DBMS（如PostGIS）的桥梁很重要，允许访问庞大的数据存储，而无需将几个千兆字节的地理数据加载到RAM中，从而可能导致R会话崩溃。本节的其余部分将展示如何基于PostGIS in Action, Second Edition 中的“Hello real world”从R调用PostGIS。^[感谢Manning Publications、Regina Obe和Leo Hsu允许使用这个例子。]

以下代码需要有效的互联网连接，因为我们正在访问一个位于QGIS Cloud (https://qgiscloud.com/)中的PostgreSQL/PostGIS数据库。^[QGIS Cloud让您将地理数据和地图存储在云中。在后台，它使用QGIS Server和PostgreSQL/PostGIS。这样，读者可以在无需在本地机器上安装PostgreSQL/PostGIS的情况下跟随PostGIS示例。感谢QGIS Cloud团队托管这个例子。]
我们的第一步是通过提供数据库的名称、主机名和用户信息来创建与数据库的连接。

1
2
3
4
5
6
library(RPostgreSQL)
conn = dbConnect(
  drv = PostgreSQL(),
  dbname = "rtafdf_zljbqm", host = "db.qgiscloud.com",
  port = "5432", user = "rtafdf_zljbqm", password = "d3290ead"
)

我们的新对象conn只是我们的R会话与数据库之间的已建立的链接。它不存储任何数据。

经常的第一个问题是：“数据库中可以找到哪些表？”可以使用dbListTables()如下回答：

1
2
3
dbListTables(conn)
#> [1] "spatial_ref_sys" "topology"        "layer"           "restaurants"
#> [5] "highways"

答案是五个表。在这里，我们只对restaurants和highways表感兴趣。前者代表美国的快餐餐厅的位置，后者是美国的主要公路。要了解表中可用的属性，我们可以运行dbListFields：

1
2
3
dbListFields(conn, "highways")
#> [1] "qc_id"        "wkb_geometry" "gid"          "feature"
#> [5] "name"         "state"

现在，由于我们知道了可用的数据集，我们可以执行一些查询——向数据库提出一些问题。查询需要用数据库能理解的语言提供——通常是SQL。第一个查询将从highways表中选择马里兰州（MD）的US Route 1。注意，如果read_sf()提供了与数据库的开放连接和查询，它就可以从数据库中读取地理数据。此外，read_sf()需要知道哪一列代表几何体（这里是：wkb_geometry）。

1
2
3
4
5
6
query = paste(
  "SELECT *",
  "FROM highways",
  "WHERE name = 'US Route 1' AND state = 'MD';"
)
us_route = read_sf(conn, query = query, geom = "wkb_geometry")

这将产生一个名为us_route的MULTILINESTRING类型的sf对象。

如前所述，我们还可以不仅提出非空间问题，而且可以基于数据集的空间属性查询。为了展示这一点，下一个示例在选定的公路周围增加了35公里（35,000米）的缓冲区。

1
2
3
4
5
6
query = paste(
  "SELECT ST_Union(ST_Buffer(wkb_geometry, 35000))::geometry",
  "FROM highways",
  "WHERE name = 'US Route 1' AND state = 'MD';"
)
buf = read_sf(conn, query = query)

注意，这是一个使用您应该已经熟悉的函数（ST_Union()，ST_Buffer()）的空间查询。您还可以在sf包中找到它们，尽管在这里它们是用小写字母写的（st_union()，st_buffer()）。实际上，sf包的函数名称在很大程度上遵循PostGIS的命名惯例。^[前缀st代表空间/时间。]

最后的查询将找到所有Hardee’s餐厅（标记为HDE），它们位于所选公路周围的35公里缓冲区内。

1
2
3
4
5
6
7
8
9
10
11
12
13
query = paste(
  "SELECT *",
  "FROM restaurants r",
  "WHERE EXISTS (",
  "SELECT gid",
  "FROM highways",
  "WHERE",
  "ST_DWithin(r.wkb_geometry, wkb_geometry, 35000) AND",
  "name = 'US Route 1' AND",
  "state = 'MD' AND",
  "r.franchise = 'HDE');"
)
hardees = read_sf(conn, query = query)

最后，按照以下方式关闭数据库连接是一种良好的做法：^[在这里关闭连接非常重要，因为QGIS云（免费版本）仅允许同时进行十个连接。]

1
RPostgreSQL::postgresqlCloseConnection(conn)


Visualization of the output of previous PostGIS commands showing the highway (black line), a buffer (light yellow) and four restaurants (red points) within the buffer.

与PostGIS不同，sf只支持空间矢量数据。

要查询和操作存储在PostGIS数据库中的栅格数据，请使用rpostgis包，或者使用诸如rastertopgsql这样的命令行工具，该工具是PostGIS安装的一部分。

这个小节仅是PostgreSQL/PostGIS的简要介绍。尽管如此，我们还是想鼓励在空间数据库管理系统（DBMS）中存储地理和非地理数据的做法，而只将那些需要进行进一步的地理统计分析的数据子集附加到R的全局环境中。有关所提及的SQL查询的更详细描述以及关于PostgreSQL/PostGIS的更全面介绍。PostgreSQL/PostGIS是一种出色的开源空间数据库选择。SQLite/SpatiaLite轻量级数据库引擎以及在背景中使用SQLite的GRASS同样也是如此。

如果您的数据集对于PostgreSQL/PostGIS来说过大，并且您需要大规模的空间数据管理和查询性能，那么探索分布式计算系统上的大规模地理查询可能是值得的。这样的系统超出了本书的范围，但值得一提的是，提供此功能的开源软件确实存在。在此领域的重要项目包括GeoMesa和Apache Sedona。apache.sedona包提供了对后者的接口。

云技术和服务的桥梁

近年来，云技术在互联网上越来越突出。这也包括使用云技术来存储和处理空间数据。主要的云计算提供商（如亚马逊网络服务、微软Azure / Planetary Computer、谷歌云平台等）在其平台上提供了大量的开放地球观测数据，例如完整的Sentinel-2档案。我们可以使用R直接连接并处理这些档案中的数据，理想情况下是从同一云和地区的机器中进行。

三个有前途的发展使得在云平台上使用此类图像档案变得更加简单和高效，它们分别是时空资产目录（STAC）、云优化的GeoTIFF（COG）图像文件格式，以及数据立方体的概念。

除了托管大型数据档案，近几年还推出了许多基于云的服务来处理地球观测数据。其中包括OpenEO倡议——一种在编程语言（包括R）和各种基于云的服务之间的统一接口。

STAC，COGs，和云端的数据立方体

时空资产目录（STAC）是一个用于描述时空数据的通用描述格式，用于描述云平台上的各种数据集，包括图像、合成孔径雷达（SAR）数据和点云。除了简单的静态目录描述，STAC-API还提供了一项网络服务，以通过空间、时间和其他属性查询目录中的项目（例如图像）。在R中，rstac包允许连接到STAC-API端点并搜索项目。在下面的示例中，我们请求与预定义区域和时间相交的亚马逊网络服务上的Sentinel-2云优化GeoTIFF（COG）数据集中的所有图像。结果包含所有找到的图像及其元数据（例如云覆盖）和指向AWS上实际文件的URL。

1
2
3
4
5
6
7
8
9
10
11
12
library(rstac)
# Connect to the STAC-API endpoint for Sentinel-2 data
# and search for images intersecting our AOI
s = stac("https://earth-search.aws.element84.com/v0")
items = s |>
  stac_search(
    collections = "sentinel-s2-l2a-cogs",
    bbox = c(7.1, 51.8, 7.2, 52.8),
    datetime = "2020-01-01/2020-12-31"
  ) |>
  post_request() |>
  items_fetch()

云存储与本地硬盘不同，传统图像文件格式在基于云的地理处理中表现不佳。广义而言，云优化GeoTIFF格式是GeoTIFF的一种特定类型，使得可以高效地从云存储中仅读取图像的某些部分。
因此，以较低分辨率读取图像或读取图像的矩形子集变得更加高效。作为R用户，您无需安装任何东西即可使用COG，因为GDAL（以及使用它的任何包）已经可以与COG一起工作。然而，请注意，在浏览数据提供商的目录时，COG的可用性是一个很大的优势。

对于较大的感兴趣区域，请求的图像仍然相对难以处理：它们可能使用不同的地图投影，可能在空间上重叠，空间分辨率通常取决于光谱带。gdalcubes包可用于从单个图像中抽象出来，并创建和处理图像集合为四维数据立方体。

下面的代码展示了一个最小示例，用于从前一个STAC-API搜索返回的Sentinel-2图像创建较低分辨率（250m）的最大NDVI合成。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
library(gdalcubes)
# Filter images by cloud cover and create an image collection object
collection = stac_image_collection(items$features,
  property_filter = function(x) {
    x[["eo:cloud_cover"]] < 10
  }
)
# Define extent, resolution (250m, daily) and CRS of the target data cube
v = cube_view(
  srs = "EPSG:3857", extent = collection, dx = 250, dy = 250,
  dt = "P1D"
) # "P1D" is an ISO 8601 duration string
# Create and process the data cube
cube = raster_cube(collection, v) |>
  select_bands(c("B04", "B08")) |>
  apply_pixel("(B08-B04)/(B08+B04)", "NDVI") |>
  reduce_time("max(NDVI)")
# gdalcubes_options(parallel = 8)
# plot(cube, zlim = c(0, 1))

为了通过云覆盖来过滤图像，我们提供一个属性过滤函数，在创建图像集合时应用于每个STAC结果项。该函数接收图像的可用元数据作为输入列表，并返回一个逻辑值，只有函数返回TRUE的图像才会被考虑。在这种情况下，我们忽略了10%或更多云覆盖的图像。有关更多详细信息，请参考2021年OpenGeoHub夏季学校上展示的教程。

STAC、COGs和数据立方体的组合形成了一种云——本地工作流程，用于在云中分析（大型）卫星图像集合。这些工具已经构成了某些工具的支柱，例如sits R包，该包允许对大型地球观测数据进行土地利用和土地覆盖分类。该包从云服务中可用的图像集合中构建EO数据立方体，并使用各种机器学习算法对数据立方体进行土地分类。有关sits的更多信息，请访问https://e-sensing.github.io/sitsbook/ 或阅读相关文章。

openEO

是一个倡议，通过定义用于处理数据的通用语言来支持云服务之间的互操作性。初始想法已在r-spatial.org博客文章中描述，旨在让用户能够轻松地在云服务之间切换，并尽可能减少代码更改。标准化流程使用多维数据立方体模型作为数据的接口。
用户可以通过R、Python、JavaScript、QGIS或网络编辑器连接到8个不同的后端（请参见https://hub.openeo.org），并定义（和链接）集合上的流程。由于后端之间的功能和数据可用性存在差异，因此openeo R包会动态加载来自连接后端的可用流程和集合。然后，用户可以加载图像集合，应用和链接流程，提交作业，以及探索和绘制结果。

以下代码将连接到openEO平台后端，请求可用数据集、流程和输出格式，定义一个流程图来从Sentinel-2数据计算最大NDVI图像，并在登录到后端后最终执行该图。openEO平台后端包括一个免费层次，可以通过现有的机构或社交平台帐户进行注册。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
library(openeo)
con = connect(host = "https://openeo.cloud")
p = processes() # load available processes
collections = list_collections() # load available collections
formats = list_file_formats() # load available output formats
# Load Sentinel-2 collection
s2 = p$load_collection(
  id = "SENTINEL2_L2A",
  spatial_extent = list(
    west = 7.5, east = 8.5,
    north = 51.1, south = 50.1
  ),
  temporal_extent = list("2021-01-01", "2021-01-31"),
  bands = list("B04", "B08")
)
# Compute NDVI vegetation index
compute_ndvi = p$reduce_dimension(
  data = s2, dimension = "bands",
  reducer = function(data, context) {
    (data[2] - data[1]) / (data[2] + data[1])
  }
)
# Compute maximum over time
reduce_max = p$reduce_dimension(
  data = compute_ndvi, dimension = "t",
  reducer = function(x, y) {
    max(x)
  }
)
# Export as GeoTIFF
result = p$save_result(reduce_max, formats$output$GTiff)
# Login, see https://docs.openeo.cloud/getting-started/r/#authentication
login(
  login_type = "oidc", provider = "egi",
  config = list(client_id = "...", secret = "...")
)
# Execute processes
compute_result(graph = result, output_file = tempfile(fileext = ".tif"))

练习

E1. 使用 qgisprocess 通过GRASS GIS的 r.sun 计算 system.file("raster/dem.tif", package = "spDataLarge") 区域在3月21日上午11点的全球太阳辐射。

E2. 使用 Rsagacmd 计算 system.file("raster/dem.tif", package = "spDataLarge") 的集水区域和集水区坡度。

E3. 继续在SAGA GIS部分创建的 ndvi_segments 对象上工作。从 ndvi 栅格中提取平均NDVI值，并使用 kmeans() 将它们分组到六个集群中。可视化结果。

E4. 附加 data(random_points, package = "spDataLarge") 并读取 system.file("raster/dem.tif", package = "spDataLarge") 到R中。从 random_points 随机选择一个点，并找到从该点可以看到的所有 dem 像素（提示：可以使用GRASS GIS计算视域）。可视化结果。例如，绘制一个山体阴影，你的视域输出的数字高程模型和点。另外，尝试使用 mapview。

E5. 通过系统调用为磁盘上存储的光栅文件使用 gdalinfo。你能在那里找到什么信息？

E6. 使用 gdalwarp 降低你的光栅文件的分辨率（例如，如果分辨率是0.5，将其改为1）。注意：在此练习中将使用 -tr 和 -r 标志。

E7. 从本章中介绍的QGIS Cloud中的PostgreSQL/PostGIS数据库中查询所有加利福尼亚州的高速公路。

E8. ndvi.tif 光栅（system.file("raster/ndvi.tif", package = "spDataLarge")）包含了基于2000年9月22日的Landsat数据为Mongón研究区域计算的NDVI。使用 rstac，gdalcubes 和 terra 从2020-08-01到2020-10-31为同一区域下载Sentinel-2图像，计算其NDVI，然后将其与 ndvi.tif 的结果进行比较。

-------------已经到底啦-------------