Wat zijn de ahn3_bladindex in de Hoogtekaart WFS ?

[Diede]

@Jeroen_PDOK,

Dankjewel voor je reactie! Ik zag dat ik niet goed had gezocht verder, er staan al best wat vragen over de hoogtekaart op het forum als je zoekt op “hoogte”.

Ik kom nu iets verder met de WCS, alleen krijg ik een Bad Gateway op de dsm kaart en met de dtm kaart krijg ik de melding dat een slice subset niet gemaakt mag worden. (hier mijn request).

https://geodata.nationaalgeoregister.nl/ahn3/wcs?
request=GetCoverage
&service=wcs
&version=2.0.1
&coverageid=ahn3_05m_dtm
&subset=x,http://www.opengis.net/def/crs/EPSG/0/28992(186691.7070)
&subset=y,http://www.opengis.net/def/crs/EPSG/0/28992(426869.7570)

Is er een manier om in deze WCS service op basis van een hoogte (x,y of lat, lon) een response te krijgen? Ik ben niet thuis in deze protocollen dus kan met voorstellen dat ik iets heel basis over het hoofd zie.

Verder, gaat dit een getal als response geven ? Ik ben namelijk bewust niet op zoek naar het inladen van kaarten. Ik wil een scriptje schrijven voor bestaande applicaties van derden. De applicaties moeten als input een punt doorgeven en een hoogte getal terug krijgen. Zij moeten zelf bepalen hoe ze dit cijfer zichtbaar maken in hun UI.

[Jeroen_PDOK]

@Diede graag gedaan! Een collega heeft ooit een stukje over geschreven op het forum over het starten met een WCS: Starten met AHN WCS service - #2 door daneng - AHN - Geoforum wellicht dat het nog wat aanvullende informatie biedt. Qua request zal mijn collega nog even reageren.

[antonbakker]

Hoi @Diede,

Er zit een bug in MapServer (wat PDOK gebruikt als server backend voor de WCS) die ervoor zorgt dat het niet mogelijk is een GetCoverage request te doen met de dekking van 1 pixel. Zie het draadje op de Mapserver mailing list ([mapserver-users] WCS GetCoverage request one by one pixel).

Een workaround kan zijn om een coverage van 2x2 pixels op te vragen en vervolgens de waarde uit te lezen van de pixel die overlapt met je POI.

Dit request werkt namelijk wel:

https://geodata.nationaalgeoregister.nl/ahn3/wcs?service=WCS&Request=GetCoverage&version=2.0.1&CoverageId=ahn3_05m_dsm&format=image/tiff&subset=x(196396.75,196397.75)&subset=y(450031.75,450032.75)

Merk overigens op dat de WCS een multipart response terug geeft. Ik heb hier een bash script voor het splitten van multipart responses (Bash script for splitting multipart response files in seperate files #wcs #bash #multipart · GitHub).

Zie het onderstaande bash script voor het uitlezen van de hoogte van een punt locatie (in RD coordinaten).

#!/usr/bin/env bash
x=196396.84 
y=450032.55
x_lower_bound=$(python3 -c "print(($x//0.25)*0.25)")
# retrieving coverage 2x2, pixel size is 0.5
x_upper_bound=$(python3 -c "print($x_lower_bound + 1.0)")
y_lower_bound=$(python3 -c "print(($y//0.25)*0.25)")
y_upper_bound=$(python3 -c "print($y_lower_bound + 1.0)")
curl "https://geodata.nationaalgeoregister.nl/ahn3/wcs?service=WCS&Request=GetCoverage&version=2.0.1&CoverageId=ahn3_05m_dsm&format=image/tiff&subset=x($x_lower_bound,$x_upper_bound)&subset=y($y_lower_bound,$y_upper_bound)" -o wcs_response
split-mulitpart wcs_response wcs_response_split
gdallocationinfo  wcs_response_split/coverage/out.tif -geoloc $x $y

Zie overigens ook deze vraag hier: point - How to query elevation data from "Digitales Geländemodell (DGM) NRW", maybe using WCS? - Geographic Information Systems Stack Exchange
Vermoedelijk gebruikt die WCS service ook MapServer als service backend.

[antonbakker]

@Diede zie ook deze QGIS Processing Tool (geschreven in Python) voor een voorbeeld hoe de WCS te bevragen voor een specifieke locatie.

[Diede]

@antonbakker Dankjewel voor beide antwoorden ik ga het zeker eens aandachtiger doornemen. Ik weet nog niet precies wanneer maar zal mijn resultaten weer delen.

[Diede]

@Jeroen_PDOK weer bedankt! Ik zal hem ook als start gaan gebruiken zodat ik wat meer van de context snap.

[Diede]

@antonbakker & @Jeroen_PDOK Sorry dat ik een tijd niks heb laten horen, ik ben bezig met dit experiment tussen het werken door dus het blijft nog al eens liggen. Het ideaal beeld wat ik had was dat ik een WFS of andere geo service een coördinaat kon opsturen met het onderstaande request en ik dan een hoogte cijfer kon terug krijgen (Zoals je hieronder in een versimpelt voorbeeld ziet).

Dit kan nu niet en ik vraag me af of dat niet bestaat omdat het gewoon een dom idee van mij is. Ik ben niet bekend in de geo web wereld of de GIS servers dus ik kan me voorstellen dat dit een oplossing is waar niemand op zit te wachten. Alternatief kan ik me goed voorstellen dat WFS nooit word gebruikt voor data responses en ik dus aan de verkeerde techniek zit te denken.

Een request naar een lat,lon coordinaat in Nederland (illustratief)

 https://geodata.nationaalgeoregister.nl/ahn3/
wfs?
service=wfs
&version=2.0.0
&request=GetFeature
&typeName=Hoogte
&Point=5.896700,51.827033
&outputFormat=json

Response zou zijn:

{
"Point": "5.896700,51.827033",
"hoogte": "70m"
}

[antonbakker]

@Diede nu heb ik dit topic op mijn beurt weer laten liggen, excuses daarvoor. Hoogte data wordt over het algemeen als raster data (coverage) ontsloten, aangezien voor elke locatie een hoogte waarde te bepalen is. Omdat het onmogelijk is voor elke locatie in Nederland een hoogtewaarde vast te leggen, wordt dit in een raster vastgelegd; voor elke rastercel een hoogtewaarde. Rasterbestandsformaten zijn geoptimaliseerd om rasterdata op deze wijze vast te leggen. Deze informatie als vector data op slaan zou erg inefficiënt zijn.

Hierdoor is er geen WFS (Web Feature Service) beschikbaar voor de AHN3 maar wel een WCS (Web Coverage Service). In plaats van een GetFeature request kan je een GetCoverage request doen op de WCS, zoals eerder beschreven in dit topic.

Voor zover ik weet zit er in de WCS standaard geen voorziening om op basis van een puntlocatie een waarde op te vragen, en zul je dus eerst op basis van de raster cell size en het raster origin punt de bounding box van de raster cel te berekenen van je point of interest (POI). En dan met die bbox een getcoverage request te doen. De informatie over het raster (cell size en origin), is op te vragen door middel van het describeCoverage request. Toegegeven, echt gebruiksvriendelijk is het WCS koppelvlak niet.

[Diede]

Hoi Anton, geen enkel punt dat je nu reageert, heel erg fijn om je reactie te zien.

Ik had zelf nooit stil gestaan bij het feit dat hoogte in raster cellen word vast gelegd maar met jou uitleg (en de WCS wikipedia lezen) snap ik waarom dit word gedaan. Ik snap nu dat WFS simpelweg niet bedoeld is voor hoogte data en je dus met de WCS service wil werken omdat je werkt met coverage.

Nu kan ik mijn vraag weer wat meer scherper stellen . Waarom word de informatie in format=image/tiff opgevraagd ? Ik zie in je functie get_ahn_val ook terug dat je een image/tiff moet ophalen die moet openen en heel wat stappen met nemen tot je een float value kan terug geven. Is het standaard zo dat elke WCS alleen images terug geeft die je zelf moet uitpakken om de onderliggende data te bereiken ? Althans zo leest de python code nu voor mij, maar ik ben dus een python en geo-leek.

Natuurlijk is het makkelijk praten vanaf de zijlijn voor mij, maar dat is vooral waar ik gebruiksonvriendelijkheid vind, ik snap niet waarom je als je de hoogte waarde van een raster cell wil weten je een image moet downloaden en uit moet kleden en de WCS je niet direct de float waarde laat opvragen ? Ik zie vast een hoop over het hoofd.

Dank nog voor de uitleg dat je niet op punten kan zoeken, klinkt als een design keuze van het protocol. Het is minder toegankelijk omdat je als leek niet snel verschil ziet in discrete data en coverage data, maar de methode sluit weer beter aan bij het model er achter waardoor je daar meer mee kan.

[antonbakker]

Hoi Diede, om data op te halen uit een WCS, dien je een GetCoverage request te doen. Zoals de naam al suggereert krijg je dan een coverage (raster data) in de response terug. Om in te gaan op je eerste vraag:

@Diede

Waarom word de informatie in format=image/tiff opgevraagd ?

Hier wordt het format image/tiff gebruikt omdat dit het enige formaat is wat de service ondersteund om de hoogtewaardes (de data) op te halen:

Nu zul je denken, huh, wat is het verschil tussen image/tiff en die andere output formaten? Het verschil is dat een WCS op een image/tiff request altijd met een GeoTIFF zal antwoorden. Omdat een GeoTIFF altijd ook een geldige TIFF bestand is hebben ze mimetype image/tiff gehandhaafd, wat nogal verwarrend is (zie ). Zie ook deze discussie in het repository van de GeoTIFF spec (waarschuwing: gaat de diepte in).

Dus bij het opvragen van een GeoTIFF (in een multipart response) krijg je een raster bestand terug met de daadwerkelijke hoogtewaarde’s per raster cel. Bij het opvragen van de overige bestandsformaten, bijv. een PNG, krijg je een raster bestand terug met kleurwaardes voor elke raster cel (ook wel bekend als een plaatje).

Dat gezegd hebbende, de lijst van outputformats op de WCS is voor verbetering vatbaar. Voor het ontsluiten van plaatjes bieden we in (bijna) alle gevallen ook een WMS service aan, dus al die image outputformats op de WCS zijn overbodig. Het aanbod beperken tot image/tiff zou veel duidelijkheid scheppen.

@Diede

ik snap niet waarom je als je de hoogte waarde van een raster cell wil weten je een image moet downloaden en uit moet kleden en de WCS je niet direct de float waarde laat opvragen ?

Om in te gaan op bovenstaande vraag; je krijgt altijd een coverage (raster data) terug uit een Web Coverage Service op een GetCoverage request. Dat is een consequentie van de aard van de data en het service protocol (WCS) wat daar omheen is gebouwd. Voor zover ik weet ondersteunt WCS geen andere manieren om data op te vragen. Kortom jouw specifieke usecase wordt niet ondersteund door het WCS protocol.

[sbjager]

Als het om 1 punt gaat, gebruik ik eenvoudigweg de WMS GetFeatureInfo. Bijvoorbeeld:

https://geodata.nationaalgeoregister.nl/ahn3/wms?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetFeatureInfo&FORMAT=image%2Fpng&TRANSPARENT=true&QUERY_LAYERS=ahn3_5m_dtm&LAYERS=ahn3_5m_dtm&INFO_FORMAT=application%2Fjson&I=50&J=50&CRS=EPSG%3A28992&STYLES=&WIDTH=101&HEIGHT=101&BBOX=166397.12000000547%2C442571.2000000054%2C209835.1999999946%2C486009.27999999456

levert de volgende json terug:

{"type": "FeatureCollection", "totalFeatures": "unknown", "features": [{"type": "Feature", "id": "", "geometry": null, "properties": {"GRAY_INDEX": 81.7406234741211}}], "crs": null}

Als je daar de GRAY_INDEX uit peutert, heb je de hoogte te pakken…

Ik gebruik dit op mijn web-kaart, een muisklik doet (als de laag een WMS is) een GetFeatureInfo Request. Dat word door mijn javascript automatisch geregeld voor alle zichtbare lagen, en als de AHN aan staat is dit wat ik terug krijg. Misschien is dat wat je zoekt?

[Diede]

Hoi Anton,

weer hartelijk dank voor je reactie. Ik zie nu pas dat jij PDOK ontwikkelaar bent dus jij bent misschien degene om dit aan te vragen.

De oplossing zoals Stefan die aanreikt is inderdaad wat ik zocht, en dit is ook wat er via bijvoorbeeld de GeoTIFF dan op te halen zou zijn.

https://geodata.nationaalgeoregister.nl/ahn3/wms?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetFeatureInfo&FORMAT=image%2Fpng&TRANSPARENT=true&QUERY_LAYERS=ahn3_5m_dtm&LAYERS=ahn3_5m_dtm&INFO_FORMAT=application%2Fjson&I=50&J=50&CRS=EPSG%3A28992&STYLES=&WIDTH=101&HEIGHT=101&BBOX=166397.12000000547%2C442571.2000000054%2C209835.1999999946%2C486009.27999999456 

Ik vroeg me hier alleen nog bij af, is het bij de WMS service mogelijk om nog verder filters toe te voegen ? zodat ik bijvoorbeeld alleen de GRAY_INDEX ophaal ?

[Diede]

@sbjager Dankjewel, dit is dus wat ik ik zocht blijkt nu :).

[sbjager]

Dat vermoedde ik al.

mbt je vraag over het antwoord: als je als info_format text/plain gebruikt, krijg je een tekst bestandje terug:

Results for FeatureType 'http://ahn3.geonovum.nl:ahn3_5m_dtm':
--------------------------------------------
GRAY_INDEX = 81.7406234741211
--------------------------------------------

Is misschien voor de software die je gebruikt wat makkelijker te parsen. Volgens de GetCapabilities van de wms heb je de keuze tussen text/plain, text/html, text/xml en application/json. Persoonlijk vind ik json het makkelijkst, maar dat ligt helemaal aan wat je wil bereiken en hoe je dat wil bereiken. Als je het antwoord direct wilt weergeven, is html misschien ook nog een optie. Dat zou je dan direct in een tekst-component kunnen gooien bijvoorbeeld.

[Anton]

@sbjager

Als het om 1 punt gaat

Deze mag wel vetgedrukt
Ik denk niet dat het handig is om dit in een loop te gebruiken of geautomatiseerd duizenden punten te verwerken in korte tijd. Dan kun je weer beter een tiff downloaden of met de WCS gaan werken.

[sbjager]

Aangepast.

En serieus: Volledig mee eens, mijn techniek is ook alleen maar bedoeld ter informatie: dus een muisklik op de kaart geeft je de hoogte ter plaatse van die muisklik. Als je een heel gebied wil hebben (iets wat ik voor WIeringen en het Waterloopbos heb gedaan) moet je gewoon de betreffende kaartbladen downloaden, da’s veel makkelijker en sneller. Maar soms is het interessant om te weten hoe hoog het een bepaalde plek ligt (als ik daar wandel of zo, of m’n broer’s achtertuin, of zoiets), en dan werkt dit prima.

[Diede]

Hoi @sbjager en @Anton,

Ik wil het graag voor een bekende van me mogelijk maken om in de applicatie Kikker op elke positie op de kaart te kunnen klikken en in de applicatie de hoogte terug te krijgen. Dit zou helpen bij zijn werk voor het plaatsen van putten voor riolering.

Ik wil voor hem ook graag uitzoeken hoe de algemene voorziening werkt en kijken of we een how-to konden maken voor verschillende GIS applicaties. In principe zullen punten hiervoor dus afdoende werken en zal de hoeveelheid requests beperkt blijven.

[antonbakker]

Wat betreft je vraag of het mogelijk is om te filteren op velden in het resultaat; dit is niet mogelijk. Een minimaler outputformaat opvragen is hier je enige optie, zoals @sbjager aangaf.

Het voordeel van GeoJSON opvragen is overigens dat de output format formeel gespecificeerd is, en dus “machine readable”. Het outputformaat info_format text/plain heeft (voor zover ik weet) geen formele spec (alhoewel er wel een PDOK conventie is).

Overigens interessant om op te merken dat ahn3 WMS invalide GeoJSON produceert: https://geodata.nationaalgeoregister.nl/ahn3/wms?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetFeatureInfo&FORMAT=image%2Fpng&TRANSPARENT=true&QUERY_LAYERS=ahn3_5m_dtm&LAYERS=ahn3_5m_dtm&INFO_FORMAT=application%2Fjson&I=50&J=50&CRS=EPSG%3A28992&STYLES=&WIDTH=101&HEIGHT=101&BBOX=166397.12000000547%2C442571.2000000054%2C209835.1999999946%2C486009.27999999456

De geometrie ontbreekt namelijk. Ik vermoed dat dit wordt veroorzaakt door een bug in de WMS programmatuur en het feit dat er rasterdata als databron wordt gebruikt in WMS.

[Diede]

Hoi @antonbakker, nog dank voor de extra aanvulling hier kan ik goed mee aan de slag