How to query elevation data from "Digitales Geländemodell (DGM) NRW", maybe using WCS?

Question

I use street network data originating from OpenStreetMap, stored in PostGIS, and want to add corresponding height/elevation data from "Digitales Geländemodell (DGM)", available for North Rhine-Westphalia/Germany from "Geobasis NRW": https://www.bezreg-koeln.nrw.de/brk_internet/geobasis/hoehenmodelle/gelaendemodell/index.html

This page lists several links to get data. As far as I understand only the downloadable xyz files and the first WCS link ("wcs_nw_dgm") contain the actual height data in m/cm.

"xyz" are huge text files, 124 MB for each square of 2 x 2 km, in total about 9000 files for North Rhine-Westphalia. I downloaded some of them and wrote a Python script to query elevation of highway nodes in them. This is working, but I would prefer to use less of download and disk capacity.

I found I can bind the WCS service in QGIS and can query single point elevations with QGIS's Info tool, so I hope this could be used to selectively query a few thousands of network nodes? But it seems WCS itself is not usable for point queries, but returns data for rectangles as tif/png/jpg image.

Answer 1

Short answer: Use gdallocationinfo.

Longer answer with more details and an alternative way to use WCS requests without GDAL, that almost works.

The WCS way

It should be possible to read a 1x1 pixel sized image from WCS and check the height value of the only pixel in the image. With this certain WCS server it does not quite work.

First check the server GetCapabilities

https://www.wcs.nrw.de/geobasis/wcs_nw_dgm?service=WCS&request=GetCapabilities

The result shows that server supports WCS version 2.0.1 and following outputformats.

<wcs:ServiceMetadata> 
<wcs:formatSupported>image/tiff</wcs:formatSupported> 
<wcs:formatSupported>image/png</wcs:formatSupported> 
<wcs:formatSupported>image/jpeg</wcs:formatSupported>
<wcs:formatSupported>image/png; mode=8bit</wcs:formatSupported>

The formats are not optimal for point queries but we must just select one and try to handle the case with it. GeoTIFF feels like best option. Next fetch more info about the coverage.

https://www.wcs.nrw.de/geobasis/wcs_nw_dgm?service=WCS&version=2.0.1&request=DescribeCoverage&coverageID=nw_dgm

<gml:Envelope srsName="http://www.opengis.net/def/crs/EPSG/0/25832" axisLabels="x y" uomLabels="m m" srsDimension="2">
<gml:lowerCorner>278000 5560000</gml:lowerCorner>
<gml:upperCorner>536000 5828000</gml:upperCorner>

With this information it is possible to try GetCoverage request that might work.

https://www.wcs.nrw.de/geobasis/wcs_nw_dgm?service=WCS&version=2.0.1&request=GetCoverage&coverageID=nw_dgm&format?image/tiff&SUBSET=x(400000)&SUBSET=y(5700000)

It does not work because server does not support slicing (in this case, taking 0 dimensional subset out of 2 dimensional raster). That's can actually be considered to be impossible with GeoTIFF or other image formats.

<ows:Exception exceptionCode="InvalidSubsetting" locator="subset">
<ows:ExceptionText>msWCSGetCoverage20(): WCS server error. Subset operation 'slice' is not supported.</ows:ExceptionText>
</ows:Exception>

Now try to get an image with equal upper limits and lower limits for x and y. This could yield a 1x1 pixel sized GeoTIFF. It does not. There may be a bug in the server software.

https://www.wcs.nrw.de/geobasis/wcs_nw_dgm?service=WCS&version=2.0.1&request=GetCoverage&coverageID=nw_dgm&format?image/tiff&SUBSET=x(400000,400000)&SUBSET=y(5700000,5700000)

<ows:Exception exceptionCode="">
<ows:ExceptionText>msImageCreate(): Image handling error. Attempt to allocate raw 
image failed, out of memory.</ows:ExceptionText>
</ows:Exception>

It should in any case be possible to request an area that is 1x1 meters. The request works and returns 1x1 pixel sized image. Unfortunately the location and pixel size are wrong when checked with gdalinfo. This must be a bug in the server software.

https://www.wcs.nrw.de/geobasis/wcs_nw_dgm?service=WCS&version=2.0.1&request=GetCoverage&coverageID=nw_dgm&format?image/tiff&SUBSET=x(400000,400001)&SUBSET=y(5700000,5700001)

Origin = (342500.000000000000000,5761000.000000000000000)
Pixel Size = (-129000.000000000000000,134000.000000000000000)

Next trial with 10x10 meter area.

https://www.wcs.nrw.de/geobasis/wcs_nw_dgm?service=WCS&version=2.0.1&request=GetCoverage&coverageID=nw_dgm&format?image/tiff&SUBSET=x(400000,400010)&SUBSET=y(5700000,5700010)

Now the result is in correct place and with 1x1 meter pixel size that is native for the WCS coverage.

Origin = (400000.000000000000000,5700010.000000000000000)
Pixel Size = (1.000000000000000,-1.000000000000000)

When the result is saved on disk it is possible to use gdallocationinfo https://gdal.org/programs/gdallocationinfo.html and query the pixel value at some georeferenced point.

gdallocationinfo wcs_nw_dgm.tif -geoloc 400001 5700001
Report:
  Location: (1P,9L)
  Band 1:
    Value: 121.169998168945

The GDAL way

Fortunately you can query the pixel value directly from WCS service with gdallocationinfo and let it do what needs to be done on the background.

This is the request. You can read the meaning of the parameters from the GDAL documentation.

gdallocationinfo "WCS:https://www.wcs.nrw.de/geobasis/wcs_nw_dgm?service=WCS&version=2.0.1&coverageID=nw_dgm" -geoloc 400001 5700001 --config GDAL_HTTP_UNSAFESSL YES --debug on

This is the debug info that gets printed and it reveals how GDAL is resolving the task.

GDAL: GDALOpen(WCS:https://www.wcs.nrw.de/geobasis/wcs_nw_dgm?service=WCS&version=2.0.1&coverageID=nw_dgm, this=00ED9BF8) succeeds as WCS.
Report:
  Location: (122001P,127999L)
  Band 1:
GDAL: GDAL_CACHEMAX = 102 MB
WCS: Requesting https://www.wcs.nrw.de/geobasis/wcs_nw_dgm?service=WCS&REQUEST=GetCoverage&VERSION=2.0.1&COVERAGEID=nw_dgm&SUBSET=x%28399856,400880%29&SUBSET=y%285700000,5700512%29&Format=image/tiff
HTTP: Fetch(https://www.wcs.nrw.de/geobasis/wcs_nw_dgm?service=WCS&REQUEST=GetCoverage&VERSION=2.0.1&COVERAGEID=nw_dgm&SUBSET=x%28399856,400880%29&SUBSET=y%285700000,5700512%29&Format=image/tiff)
HTTP: libcurl/7.37.1 OpenSSL/1.0.1h zlib/1.2.3
WCS: GDALOpenResult() on content-type: image/tiff
GDAL: GDALOpen(/vsimem/wcs/00ED9BF8/wcsresult.dat, this=0626D708) succeeds as GTiff.
GDAL: GDALClose(/vsimem/wcs/00ED9BF8/wcsresult.dat, this=0626D708)
    Value: 121.169998168945