Skip to content

使用 GDAL 工具进行加权多标准叠加分析

参考来源: (https://spatialthoughts.com/2022/01/31/multi-criteria-overlay-gdal/)

本文仅供学习交流，版权归原作者所有。

多标准叠加分析（Multi-criteria Overlay Analysis）是根据选定区域应具备的各种属性，分配土地以适应特定目标的过程。虽然这是一种常见的 GIS 操作，但在栅格空间中执行效果最佳。

本文概述了典型的操作工作流：获取源矢量数据，将其转换为适当的栅格，进行重分类，并执行数学运算以进行加权适宜性分析。本文使用的是开源 GDAL 库提供的命令行实用程序。如果你想在 QGIS 中进行此类分析，请查看我的教程：多标准叠加分析 (QGIS3)。

我们将使用伦敦市的犯罪和基础设施数据，寻找适合建设新停车设施的区域，以帮助减少自行车盗窃。我们的分析将应用以下 3 个标准。拟建的停车场必须：

  1. 位于自行车盗窃高发区
  2. 靠近自行车路线
  3. 远离现有的停车设施

问题陈述：确定建设新自行车停车设施的合适位置

想跟随本教程一起操作吗？

  • 安装 GDAL。我倾向于 使用 conda 安装。
  • 源数据图层包含在 multicriteria.gpkg 中。下载该文件并保存到你的计算机。打开终端并切换到包含该 GeoPackage 的文件夹。

第一步：矢量图层栅格化

叠加分析的第一步是将每个矢量数据图层转换为栅格。这是使用 gdal_rasterize 命令完成的。让我们先检查一下数据图层的内容。

ogrinfo multicriteria.gpkg

一个重要的考虑因素是所有栅格必须具有相同的范围（extent）。我们可以使用边界图层（boundary layer）并找出该图层 CRS 中的范围。我们将使用 gdal_rasterize 命令的 -te 选项指定此范围。

ogrinfo -so multicriteria.gpkg boundary

让我们从栅格化 cycling_routes（自行车路线）图层开始。我们希望创建一个输出栅格，其中有路线的像素值为 1，没有路线的像素值为 0。该值使用 -burn 选项指定。我们可以使用 -tr 选项将输出栅格的分辨率设置为 10 米。命令完成后，每个像素的值将根据路线是否与其相交而被设置为 1 或 0。

注意：\ 字符用于将长命令断成多行。如果你在复制/粘贴时遇到问题，请将其删除并在单行中运行命令。

gdal_rasterize -ot Int16 -burn 1 -tr 10 10 \
  -te 523843.7 177847.3 531169.1 183893.8 multicriteria.gpkg \
  -l cycling_routes routes.tif

我们可以使用相同的选项栅格化 cycle_parking（自行车停车点）图层。

gdal_rasterize -ot Int16 -burn 1 -tr 10 10 \
  -te 523843.7 177847.3 531169.1 183893.8 multicriteria.gpkg \
  -l cycle_parking existing_parking.tif

现在只需要栅格化 bicycle_thefts（自行车盗窃）图层。这是一个密集的点图层，包含个别盗窃案的位置——点与点之间存在重叠。我们可以通过将每个像素内的点总数相加并计算盗窃密度来创建栅格图层。这种映射点密度的方法称为分箱（Binning）或聚合（Aggregation）。我们可以使用 -add 选项为每个输入要素增加像素的值。

gdal_rasterize -burn 1 -add -tr 100 100 \
  -te 523843.7 177847.3 531169.1 183893.8 multicriteria.gpkg \
  -l bicycle_thefts thefts.tif -ot Int16

第二步：生成邻近度（欧几里得距离）栅格

现在我们已经转换了所有输入图层，我们可以生成邻近度栅格。这也称为欧几里得距离栅格——输出栅格中的每个像素代表到输入栅格中最近像素的距离。可以使用 -distunits 选项以 像素（pixels） 或 地理单位（geographic units） 指定距离。生成的栅格随后可用于确定距离输入一定范围内的合适区域。

gdal_proximity.py routes.tif routes_proximity.tif \
 -ot Int16 -distunits GEO


gdal_proximity.py existing_parking.tif \ 
 existing_parking_proximity.tif -ot Int16 -distunits GEO

第三步：重分类栅格值

要在叠加分析中使用邻近度栅格，我们必须首先对其进行重分类以创建离散值。这里我们使用 0-100 的范围，并根据输入范围分配适当的值。

对于代表 盗窃（Thefts） 的点密度栅格，我们为更多的盗窃分配更高的值。对于 路线（Routes） 栅格，我们为靠近路线的像素分配更高的值。最后，对于 现有停车点（Existing Parking） 栅格，我们为距离较远的像素分配更高的值。最终分配如下所示。

盗窃 (Thefts) 路线 (Routes) 现有停车点 (Existing Parking) 重分类值 (Re-Class Value)
> 20 0-50m > 100m 100
10-20 50-100m 50-100m 50
< 10 > 100m < 50m 10

以下命令使用 gdal_calc.py 获取输入栅格的连续值，并创建具有离散值的重分类栅格。

gdal_calc.py -A thefts.tif \
  --outfile thefts_reclass.tif \
  --calc="100*(A>20) + 50*(A>10)*(A<=20) + 10*(A<10)"

gdal_calc.py -A routes_proximity.tif \
  --outfile routes_reclass.tif \
  --calc="100*(A<=50) + 50*(A>50)*(A<=100) + 10*(A>100)"

gdal_calc.py -A existing_parking_proximity.tif \
  --outfile existing_parking_reclass.tif \
  --calc="100*(A>100) + 50*(A>50)*(A<=100) + 10*(A<50)"

第四步：叠加分析

现在我们准备进行最终的叠加分析。在继续之前，确保所有栅格具有相同的像素分辨率非常重要。我们对 existing_parking_reclass.tif 和 routes_reclass.tif 图层进行下采样，使其分辨率与 thefts_reclass.tif 图层匹配。

gdalwarp -r average -tr 100 100 \ 
 existing_parking_reclass.tif parking_reclass_resampled.tif
gdalwarp -r average -tr 100 100 \ 
 routes_reclass.tif routes_reclass_resampled.tif

我们现在准备好运行叠加分析并生成适宜性地图。

简单叠加 (Simple OVERLAY)

如果我们选择给予每个标准相同的权重（重要性）——我们可以简单地将所有 3 个栅格相加并将像素值归一化，使它们处于 0-100 的范围内。我们使用 --NoDataValue 选项设置输出中的无数据值。

gdal_calc.py \
  -A thefts_reclass.tif \
  -B routes_reclass_resampled.tif \
  -C parking_reclass_resampled.tif --outfile suitability.tif \
  --calc="(A + B + C)/3" --NoDataValue=0

加权叠加 (Weighted Overlay)

在加权叠加分析中，根据相对重要性为不同的标准分配权重。让我们假设靠近自行车路线并不像最终决策中的其他两个因素那么重要。我们可以给予其他标准更高的重要性，如下所示。

  • 盗窃密度：权重=2
  • 路线：权重=1
  • 现有停车点：权重=2

我们可以使用类似的表达式，但必须注意通过将总和除以所有权重的总和来归一化输出。

gdal_calc.py \
  -A thefts_reclass.tif \
  -B routes_reclass_resampled.tif \
  -C parking_reclass_resampled.tif --outfile suitability.tif \
  --calc="(2*A + B + 2*C)/5" --NoDataValue=0

生成的栅格代表了根据所选标准及其相对权重计算出的新自行车停车设施的适宜性。

你还可以通过将命令放入 Shell 脚本并运行它来作为批处理过程运行此分析。下面是运行加权多标准叠加分析的完整脚本。

gdal_rasterize -burn 1 -add -tr 100 100 -te 523843.7 177847.3 531169.1 183893.8 multicriteria.gpkg -l bicycle_thefts thefts.tif -ot Int16
gdal_rasterize -ot Int16 -burn 1 -tr 10 10 -te 523843.7 177847.3 531169.1 183893.8 multicriteria.gpkg -l cycling_routes routes.tif
gdal_rasterize -ot Int16 -burn 1 -tr 10 10 -te 523843.7 177847.3 531169.1 183893.8 multicriteria.gpkg -l cycle_parking existing_parking.tif
gdal_proximity.py routes.tif routes_proximity.tif -ot Int16 -distunits GEO
gdal_proximity.py existing_parking.tif existing_parking_proximity.tif -ot Int16 -distunits GEO
gdal_calc.py -A thefts.tif --outfile thefts_reclass.tif --calc="100*(A>20) + 50*(A>10)*(A<=20) + 10*(A<10)"
gdal_calc.py -A routes_proximity.tif --outfile routes_reclass.tif --calc="100*(A<=50) + 50*(A>50)*(A<=100) + 10*(A>100)"
gdal_calc.py -A existing_parking_proximity.tif --outfile existing_parking_reclass.tif --calc="100*(A>100) + 50*(A>50)*(A<=100) + 10*(A<50)"
gdalwarp -r average -tr 100 100 existing_parking_reclass.tif parking_reclass_resampled.tif
gdalwarp -r average -tr 100 100 routes_reclass.tif routes_reclass_resampled.tif
gdal_calc.py -A thefts_reclass.tif -B routes_reclass_resampled.tif -C parking_reclass_resampled.tif --outfile suitability.tif --calc="(2*A + B + 2*C)/5" --NoDataValue=0

有兴趣学习 GDAL 工具吗？查看我的免费课程 Mastering GDAL Tools。

数据致谢

本文中使用的数据集来自以下开放数据门户。

  • 自行车停车点 (Bicycle Parking): Cycling Infrastructure data: Public Tfl data. 包含 OS 数据 © Crown copyright and database rights 2016 以及 Geomni UK Map data © and database rights [2019]
  • 自行车盗窃 (Cycling Thefts): 伦敦个人犯罪和反社会行为 (ASB) 事件，包括街道级位置信息以及与犯罪相关的后续警察和法院结果。由 Single Online Home National Digital Team 根据 Open Government Licence v3.0 发布。下载自 https://data.police.uk/
  • 威斯敏斯特边界 (Westminster Boundary): Statistical GIS Boundary Files for London. 包含 National Statistics data © Crown copyright and database right [2015], 包含 Ordnance Survey data © Crown copyright and database right [2015]