AHN WCS opvragen

Datasets AHN
#1

Ik probeer via python in QGIS de AHN op te vragen via WCS.

Via een script wordt de bounding box vast gesteld, die is dus elke keer van een andere grootte. De url die ik vervolgens gebruik om de wcs te raadplegen is:

https://geodata.nationaalgeoregister.nl/ahn3/wcs?SERVICE=WCS&VERSION=1.0.0&REQUEST=GetCoverage&FORMAT=GEOTIFF_FLOAT32&COVERAGE=ahn3_05m_dtm&crs=EPSG:28992&response_crs=EPSG:28992&BBOX=76481.0581597727577901,438908.8791885763639584,81054.0107515793206403,440199.2156986797926947&WIDTH=9146&HEIGHT=2581

Op dit moment haal ik de width en height uit de coördinaten van de bounding box (*2 want het is een 0,5m raster), maar dit werkt helaas niet. Wat zie ik over het hoofd?

Is er een manier om van te voren de width en height te bepalen aan de hand van de bounding box? Ik wil de afbeelding in de hoogst mogelijke resolutie (dus 0,5m) ophalen.

#2

Hoi @MennodR,

Het probleem is dat de afmeting van het opgevraagde raster groter is dan de door de service toegestane afmeting. Dit is te zien in de response van je request:

$ curl "https://geodata.nationaalgeoregister.nl/ahn3/wcs?SERVICE=WCS&VERSION=1.0.0&REQUEST=GetCoverage&FORMAT=GEOTIFF_FLOAT32&COVERAGE=ahn3_05m_dtm&crs=EPSG:28992&response_crs=EPSG:28992&BBOX=76481.0581597727577901,438908.8791885763639584,81054.0107515793206403,440199.2156986797926947&WIDTH=9146&HEIGHT=2581"

<?xml version='1.0' encoding="UTF-8" ?>
<ServiceExceptionReport version="1.2.0"
xmlns="http://www.opengis.net/ogc" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:schemaLocation="http://www.opengis.net/ogc https://geodata.nationaalgeoregister.nl/wcs/1.0.0/OGC-exception.xsd">
  <ServiceException code="InvalidParameterValue" locator="width/height">msWCSGetCoverage(): WCS server error. Raster size out of range, width and height of resulting coverage must be no more than MAXSIZE=4000.
  </ServiceException>
</ServiceExceptionReport>

Raster size out of range, width and height of resulting coverage must be no more than MAXSIZE=4000.

Verder zie ik geen problemen met deze aanpak, het zou gewoon moeten werken bij kleinere afmetingen.
Als je bbox groter is dan 4000 pixels in ten minste 1 dimensie dien je de bbox verder op te splitsen in kleinere requests waarbij de width en height niet groter zijn dan 4000 pixels.

Verder zou je in je python script de Content-type response header kunnen checken om te verifieren of je een geotiff terugkrijgt of een service-exception. Want in het geval van een service exceptie krijg je gewoon een HTTP 200 OK statuscode terug.

  • service exceptie respone: Content-Type: application/vnd.ogc.se_xml; charset=UTF-8
  • geotif response: Content-Type: image/tiff

Daarnaast raad ik je aan gebruik te maken van de nieuwe AHN3 WCS url: https://service.pdok.nl/rws/ahn3/wcs/v1_0?service=wcs&request=getcapabilities&version=1.0.0 . De oude wordt over iets minder dan een half jaar uitgezet, zie ook Nieuwe URL’s AHN3 en Luchtfoto labels na cloudmigratie - Datasets - Geoforum

3 likes
#3

Bedankt voor de snelle reactie Anton!

Het werkt inderdaad wel bij afmetingen onder de 4000 pixels. Ik loop alleen tegen een nieuw probleem aan, wat meer QGIS en wat minder PDOK gerelateerd: ik kan de tif afbeelding wel gewoon ophalen als ik de url via mijn browser bezoek, maar het is mij nog niet gelukt om deze in QGIS met python toe te voegen.

Mijn code is nu:
uri = “https://service.pdok.nl/rws/ahn3/wcs/v1_0?service=wcs&VERSION=1.0.0&REQUEST=GetCoverage&FORMAT=GEOTIFF_FLOAT32&COVERAGE=ahn3_05m_dtm&crs=EPSG:28992&response_crs=EPSG:28992&BBOX=88013.7238158045802265,441478.6266455391887575,90700.2187192221754231,442123.9031226849183440&WIDTH=2685&HEIGHT=644”
rlayer = QgsRasterLayer(uri, ‘ahn’, ‘wcs’)
QgsProject.instance().addMapLayer(rlayer)

Dit resulteert echter in een laag met ‘onbekende’ databron. Heeft er iemand een tip over hoe ik dit probleem op los?

#4

Ik houd me ver van python, maar

Het lijkt er op dat je een raster laag wilt toevoegen. Kun je niet beter een WCS laag gebruiken? in de UI is dat een aparte optie, vandaar dat ik er van uit ga dat dat in python ook zo is. Maar zoals gezegd, ik doe niks met python, dus kan het heel goed mis hebben.

#5

Volgens mij wil @MennodR juist rasters downloaden met van tevoren gedefinieerde afmetingen, om deze te combineren in raster bestanden die groter zijn dan 4000x4000 pixels. Maar corrigeer me als ik het bij het verkeerde eind heb. In dat geval is een WCS laag juist niet handig, en zal je de rasters eerst moeten downloaden, en dat de gedownloade bestanden als raster laag aan je project toevoegen. Je kan namelijk niet direct een WCS raster URL als rasterlaag-bron gebruiken.

Ik heb hier de de code van de pdokserviceplugin gebruikt. Het voordeel daarvan is dat de requests door de netwerk stack van QGIS gaan, wat handig is als je bijv. een proxy hebt ingesteld in QGIS. Ook kan je de requests inspecteren in de QGIS Network Logger plugin, wat erg handig is voor debugging:

from qgis.core import QgsBlockingNetworkRequest
from qgis.PyQt.QtNetwork import QNetworkRequest
from qgis.PyQt.QtCore import QUrl
from tempfile import NamedTemporaryFile
"""
Download file to tempfile with extension. Based on https://github.com/rduivenvoorde/pdokservicesplugin/blob/master/pdokservicesplugin/lib/http_client.py
"""
def dowload_file(url, extension):
    qgs_request = QgsBlockingNetworkRequest()
    request = QNetworkRequest(QUrl(url))
    _ = qgs_request.get(request, True) # TODO: implement error handling
    reply = qgs_request.reply()
    _ = reply.error() # TODO: implement error handling
    f = NamedTemporaryFile(suffix=extension, delete=False) # https://docs.python.org/3/library/tempfile.html#tempfile.NamedTemporaryFile -> Whether the name can be used to open the file a second time, while the named temporary file is still open, varies across platforms (it can be so used on Unix; it cannot on Windows) [kan voor problemen zorgen, niet getest op windows]
    f.write(reply.content())
    return f.name

url = "https://service.pdok.nl/rws/ahn3/wcs/v1_0?service=wcs&VERSION=1.0.0&REQUEST=GetCoverage&FORMAT=GEOTIFF_FLOAT32&COVERAGE=ahn3_05m_dtm&crs=EPSG:28992&response_crs=EPSG:28992&BBOX=88013.7238158045802265,441478.6266455391887575,90700.2187192221754231,442123.9031226849183440&WIDTH=2685&HEIGHT=644"
f_name = dowload_file(url, ".tif")
iface.addRasterLayer(f_name, "ahn3")
2 likes
#6

Bedankt Anton! Dat was inderdaad precies wat ik nodig had.

2 likes