如何在 GRASS GIS 中下载 AHN 数据

AHN

“Actueel Hoogtebestand Nederland”（AHN，版本 4）是荷兰最新的数字高程模型 (DEM)。 它包括数字地形模型（DTM）和数字表面模型（DSM）。 前者代表没有任何特征（例如建筑物、植被或其他结构）的裸露地面， 而后者则包括这些特征。

两个图层都能够以 6.5 x 5 公里的图块以 0.5 或 5 米的分辨率下载。 下载一个图块（大约 500 MB）是可行的。 但是，如果需要更多图块， 可以说更简单的方法是使用 WCS （网络覆盖服务）协议下载数据， 该协议允许根据空间限制选择下载可用数据的哪些部分。

导入 DTM

此示例显示了为 Land van Cuijk 下载 0.5 米分辨率 DTM 的步骤， 研究小组 Climate-robust Landscapes 正在该市开展多项研究。 第一步是获取 r.in.wcs 插件。 通过此插件， 可使用 WCS 协议下载数据。

步骤1

安装 r.in.wcs 插件

第一步是导入 Python 库。 请注意， 接下来的脚本不会重复此操作。

import grass.script as gs

现在，可 r.in.wcs 使用 g.extension 函数安装插件。

gs.run_command("g.extension", extension="r.in.wcs")

r.in.wcs 插件使用当前区域的分辨率和范围导入栅格图层。 因此， 关键的一步是设置该区域的范围和分辨率， 使其与 AHN 层的范围和分辨率完美对齐。

图形用户界面

需要 g.extension 安装插件的功能。 在主菜单中， 转至 Settings > Addons extension > Install extension from addon。g.extension或者， 在命令行中输入， 这将打开如下所示的窗口。

g.extension 函数， 从命令行启动。

步骤2

设置范围和分辨率以匹配 AHN 的范围和分辨率。 接下来， 下载荷兰城市的行政边界， 并提取“Land van Cuijk”的边界。

gs.run_command("g.region",
    n=618750, 
    s=306250, 
    w=10000, 
    e=280000, 
    res=0.5)

图形用户界面

g.region 在命令行或控制台中输入。 您还可以在下找到该功能 main menu > Settings > Computational region > set region 。 这将打开以下屏幕：

设定界限

设置分辨率

通过此设置， 可以将该区域与 AHN 数据集对齐。 接下来， 下载荷兰城市的行政边界， 并提取“Land van Cuijk”的边界。

步骤3

下载带有邻里行政边界的图层。 现在，可以设置区域以匹配矢量图层的范围， 同时注意该区域与 DTM 完美对齐。

gs.run_command(
    "v.in.wfs",
    url="https://service.pdok.nl/cbs/wijkenbuurten/2022/wfs/v1_0?",
    output="municipalities",
    name="gemeenten",
)

接下来，提取“Land van Cuijk”市的边界

gs.run_command(
    "v.extract",
    input="municipalities",
    where="naam = 'Land van Cuijk'",
    output="LandvanCuijk",
)

现在，可以设置区域以匹配矢量图层的范围， 同时注意该区域与 DTM 完美对齐。

图形用户界面

v.in.wfs 在命令行或控制台中输入。 还可以在下找到该功能 main menu > File > Import vector data。 这将打开以下屏幕（需要在两个选项卡中填写参数）：

下载带有城市边界的矢量图层。 定义基本 URL 和输出层的名称。

下载带有城市边界的矢量图层， 填写要下载的 WFS 图层的名称。

提取 Land van Cuijk 市的边界。 选择带有城市的矢量图层的名称， 并给出输出层的名称。

步骤4

将区域设置为正确的范围和分辨率。

首先，获取矢量图层边界框和区域范围的北、南、东、西坐标。

region_tiles = gs.parse_command("g.region", 
    flags="gu")
region_munic = gs.parse_command("v.info", 
    flags="g", 
    map="LandvanCuijk")

现在，调整矢量图层边界框的北边界，使其与 AHN 切片对齐。

from math import floor

n = float(region_tiles["n"]) - floor(
    (float(region_tiles["n"]) - float(region_munic["north"]))
)
s = float(region_tiles["s"]) + floor(
    (float(region_munic["south"]) - float(region_tiles["s"]))
)
w = float(region_tiles["w"]) + floor(
    (float(region_munic["west"]) - float(region_tiles["w"]))
)
e = float(region_tiles["e"]) - floor(
    (float(region_tiles["e"]) - float(region_munic["east"]))
)
gs.run_command("g.region", n=n, s=s, e=e, w=w)

最后，导入 DTM 层，这可能需要一些时间。

步骤5

导入定义区域的 DTM。

gs.run_command(
    "r.in.wcs",
    url="https://service.pdok.nl/rws/ahn/wcs/v1_0?",
    coverage="dtm_05m",
    urlparams="GetMetadata",
    output="dtm_05",
)

最大浮点值 ( 3.4028235e+38 ) 不被识别为 NULL 值。 因此，使用 r.nullsetnull 函数中的参数来指定该值必须设置为 NULL。

gs.run_command("r.null", 
    map="dtm_05", 
    setnull=3.4028235e+38)

如果这不起作用， 可以使用 r.mapcalc 函数将此值设置为 NULL。

gs.run_command("r.mapcalc", 
    expression="dtm_05m = if(dtm_05m > 1000,null(),dtm_05m)",
    overwrite=True)

为了使它更容易一些， 将步骤 2、4 和 5 合并到一个新的插件 r.in.ahn 中。 它的边缘有点粗糙（它只适用于具有 CRS RD New 的位置）， 但它达到了它的目标。 已经将其提交到 GRASS GIS 插件存储库， 所以希望它能在一周左右的时间内可用。

脚注

应该熟悉 GRASS GIS 以及GRASS GIS 中使用的区域概念。 如果是 GRASS GIS 的新手， 强烈建议首先查看GRASS GIS 快速入门以及有关 GRASS GIS 数据库的说明。

下载整个城市的 DTM 需要一段时间。 如果想加快速度， 可以使用自己的矢量数据来处理较小的区域。

本文内容来源于网站：https://ecodiv.earth/post/download_AHN/， 由小编整理编译。