摘要: AHN “Actueel Hoogtebestand Nederland”（AHN，版本 4）是荷兰最新的数字高程模型 (DEM)。它包括数字地形模型（DTM）和数字表面模型（DSM）。前者代表没有任何特征（例如建筑物、植被或其他结构）的裸露地面，而后者则包...
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/，由小编整理编译。
