Download AHN data in GRASS GIS

dem
dtm
dsm
osgeo
open data
grass gis
data import
wcs
wfs
In the Netherlands, we are lucky to have access to 0.5 meter resolution DEM. You can download the file in tiles of 6.5 by 5 kilometers, or, arguably easier, you can use the WCS protocol to download the data for arbitrary areas. This post shows you how to do this in GRASS GIS.
Author

Paulo van Breugel

Published

March 17, 2024

The AHN

The ‘Actueel Hoogtebestand Nederland’ (AHN, version 4) is the latest digital elevation model (DEM) of the Netherlands. It includes a digital terrain model (DTM) and a digital surface model (DSM). The former represents the bare ground surface without any features such as buildings, vegetation, or other structures, while the latter includes these features.

Both layers can be downloaded at 0.5 or 5 meter resolution in tiles of 6.5 by 5 km. Downloading one tile, which is approximately 500 MB, is doable. However, if you need more tiles, the arguably easier way is to download the data using the WCS (Web Coverage Service) protocol, which allows you to choose what part of the available data to download based on spatial constraints.

Import the DTM

This example shows the steps to download the 0.5-meter resolution DTM for the Land van Cuijk, a municipality in which our research group Climate-robust Landscapes is carrying out several studies. The first step is to get the r.in.wcs addon. With this addon, you can use the WCS protocol to download data 1 2.

Install the r.in.wcs addon.

The first step is to import the Python libraries. Note, this will not be repeated with the next scripts.

import grass.script as gs

Now, you can install the r.in.wcs addon using the g.extension function.

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

You need the g.extension function to install addons. In the main menu, go to Settings > Addons extension > Install extension from addon. Alternatively, type in g.extension on the command line. This will open the window shown below.

The g.extension function, started on the command line.

The g.extension function, started on the command line.

The r.in.wcs addon imports the raster layers using the resolution and extent of the current region. A crucial step, therefore, is to set the extent and resolution of the region so that it aligns perfectly with that of the AHN layer.

Set the extent and resolution to match those of the AHN.

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

Type in g.region on the command line or in the console. You can also find the function under main menu > Settings > Computational region > set region. This opens the following screen:

Set the bounds

Set the bounds

Set the resolution

Set the resolution

With these settings, you align the region with the AHN data set.

Next, download the administrative boundaries of the Dutch municipalities, and extract the boundaries of the “Land van Cuijk”.

Download layer with administrative boundaries of the neighborhood.

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

Next, extract the boundaries of the municipality of “Land van Cuijk”

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

Type in v.in.wfs on the command line or in the console. You can also find the function under main menu > File > Import vector data. This opens the following screen (you need to fill in parameters in two tabs):

Download the vector layer with the municipality boundaries. Define the base URL and the name of the output layer.

Download the vector layer with the municipality boundaries. Define the base URL and the name of the output layer.

Download the vector layer with the municipality boundaries. Fill in the name of the WFS layer to download.

Download the vector layer with the municipality boundaries. Fill in the name of the WFS layer to download.

Extract the boundaries of the municipality of Land van Cuijk. Select the name of the vector layer with municipalities and give the name of the output layer.

Extract the boundaries of the municipality of Land van Cuijk. Select the name of the vector layer with municipalities and give the name of the output layer.

Extract the boundaries of the municipality of Land van Cuijk. Fill in the query, which defines which features you want to select and save. Tip: Use the convenient query builder.

Extract the boundaries of the municipality of Land van Cuijk. Fill in the query, which defines which features you want to select and save. Tip: Use the convenient query builder.

Now, you can set the region to match the extent of the vector layer, while taking care that the region aligns perfectly with the DTM.

Set the region to the correct extent and resolution.

First, get the north, south, east, and west coordinates of the bounding box of the vector layer and of the region’s extent.

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

Now, adjust the northern boundary of the bounding box of the vector layer so that it aligns with the AHN tiles.

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)

And finally, import the DTM layer. Warming, this may take some time.

Import the DTM for the defined region.

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

The maximum floating-point value (3.4028235e+38) is not recognized as a NULL value. So use the setnull parameter in the r.null function to specify that this value has to be set to NULL.

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

If this does not work (it didn’t for me), you can use the r.mapcalc function to set this value to NULL.

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

Fill in the WCS service URL

Fill in the WCS service URL

Get the overage name to request

Get the overage name to request

To make it a bit easier, I combined steps 2, 4, and 5 in a new addon r.in.ahn. It is a bit rough around the edges (it only works in a Location with CRS RD New), but it serves its goal. See the next post for a brief description how to use this addon.

Footnotes

  1. You are expected to be familiar with GRASS GIS and the concept of region used in GRASS GIS. If you are new to GRASS GIS, you are warmly recommended to first check out the GRASS GIS Quickstart and the explanation about the GRASS GIS database.↩︎

  2. Downloading the DTM for the whole municipality will take a while. If you want to speed things up, you can work with a smaller area by using your own vector data.↩︎