Hoogtewaarden AHN4

Datasets AHN
1

Beste,

Voor de uitbereiding van mijn applicatie ben ik op zoek naar een manier om via een service de hoogtewaarden uit de AHN4 0,5m op te vragen. Via de WCS DTM_05m krijg ik met mijn beperkte kennis vooralsnog alleen een .tif, jpeg, of .png bestand dat mijn programma (Rhino/grasshopper) niet kan interpeteren.

Kan ik op een of andere manier de ‘originele’ 0,5x0,5 hoogtewaarden achterhalen zonder dat dat meteen tot een afbeelding leidt?

Ik gebruik nu voor https://service.pdok.nl/rws/ahn/wcs/v1_0 GET
{
“SERVICE”: “WCS”,
“REQUEST”: “GetCoverage”,
“VERSION”: “2.0.1”,
“COVERAGEID”: “dtm_05m”,
“CRS”: “http://www.opengis.net/def/crs/EPSG/0/28992”,
“FORMAT”: “image/png”,
“SCALESIZE”: “x(1000),y(1000)”
}

Alvast bedankt voor enige hulp, met vriendelijke groet,

Thomas

1 like
2

Hoi ik zou de wms gebruiken. Hier een voorbeeld in Python

`from shapely.geometry import Polygon
from owslib.wms import WebMapService
from PIL import Image
import io
import numpy as np

Stap 1: Definieer de polygon en bereken de bounding box

polygon_points = [
(126441.75, 502316.28),
(126490.47, 502283.52),
(126483.33, 502271.34),
(126434.61, 502302.84),
(126441.75, 502316.28)
]

polygon = Polygon(polygon_points)
minx, miny, maxx, maxy = polygon.bounds

Stap 2: Gebruik van de bounding box om een GetMap-verzoek naar de WMS-service te sturen

wms_url = "https://service.pdok.nl/rws/ahn/wms/v1_0?"
wms = WebMapService(wms_url)

layer = ‘dtm_05m’
srs = ‘EPSG:28992’
bbox = (minx, miny, maxx, maxy)
size = (256, 256) # Pas deze grootte aan indien nodig

img = wms.getmap(
layers=[layer],
srs=srs,
bbox=bbox,
size=size,
format=‘image/png’,
transparent=True
)

Stap 3: Ophalen van de afbeelding en verwerken van de gegevens

image = Image.open(io.BytesIO(img.read()))
image_array = np.array(image)

Verwerken van de hoogtegegevens (dit is een vereenvoudigd voorbeeld)

Hier gaan we ervan uit dat de hoogtegegevens als grijstinten worden weergegeven

en dat elke pixel een hoogte vertegenwoordigt

heights = image_array[:, :, 0] # Neem de R-kanaal van de afbeelding (bijv. grijstinten)

Analyseren van de hoogteverschillen

min_height = heights.min()
max_height = heights.max()
height_difference = max_height - min_height

print(f"Minimale hoogte: {min_height}")
print(f"Maximale hoogte: {max_height}")
print(f"Hoogteverschil: {height_difference}")`

3

Hoi Carlo,

Dank voor je uitgebreide antwoord.

Ik merk dat ik mijn vraag beter moet duiden. In jouw voorbeeld wordt data uit de png gehaald om die vervolgens te verwerken. Echter, ik zou het graag geheel zonder png/jpg willen doen.

Ik vroeg me dus af of de originele (cloudpoint-)data die tot de png heeft geleid ook via een service ontsloten wordt? Of zo nee, of dit niet voor meer mensen wenselijk zou zijn?

Met vriendelijke groet,

Thomas Broos

4

Hoi @ThomasBroos

Je kunt de ahn het beste downloaden als geotiff. Bijvoorbeeld met de volgende request:

https://service.pdok.nl/rws/ahn/wcs/v1_0?SERVICE=WCS&VERSION=1.0.0&REQUEST=GetCoverage&FORMAT=image/tiff&COVERAGE=dsm_05m&crs=EPSG:28992&BBOX=206964.365,474224.129,207264.365,474524.129&WIDTH=601&HEIGHT=601

De hoogte data kun je met uitlezen met de python bibliotheek rasterio.
https://rasterio.readthedocs.io/en/stable/

5

Het zou mogelijk moeten zijn een puntwaarde met een GetFeatureInfo-verzoek te bevragen. Praktisch probleem is wel dat de coördinaten ten opzichte van de WMS-tegel zijn. Misschien is dat op te lossen door een tegel van 1×1 pixel op te vragen en daarvan de waarde van de enige pixel?

6

@ebeere dank voor de suggestie! Helaas ben ik dan weer met een image format aan de gang, hoewel het inderdaad geen png/jpg is. Mijn verdere verwerking heeft daar allemaal problemen mee. Dus ik ben geloof ik op zoek naar een door mensen begrijpbare string/integer/float output.

@geotiles
Dank! Dit komt tot nu toe het dichtst in de buurt bij wat ik zoek. Ik heb nu voor een locatie in Deventer een hoogtewaarde van ± 8.2 meter gekregen. Daar ben ik erg gelukkig mee.
De situatie is wel zo dat ik eigenlijk voor een oppervlakte van, zeg 200x200m, al die waarden zou willen opvragen. Is daar met deze manier ook een methode voor?
Anders ben ik denk ik genoodzaakt om zeg per m2 een request te doen, wat een heleboel requests oplevert… en wellicht een ban…@pdokbeheer ?

Mvg, Thomas

7

Hoi @ThomasBroos uit de geotiff kun je float destileren met rasterio.

image
image1120×533 92.9 KB

8

Daar is de WCS-service voor. Het antwoord is een GeoTIFF met in float32 de hoogte ten opzichte van NAP.

2 likes
9

Besten, graag zou ik het nog één keer willen proberen:

Ik probeer via een service voor een gebied -van zeg 300mx300m- de hoogtepunten te ontsluiten.

  • het WMS 1x1 pixel antwoord van @ebeere kwam het dichtst bij, maar levert me absurd veel aanvragen en reactietijd op voor iets waarvan ik denk dat het simpeler zou moeten kunnen.

  • Zodra het antwoord op image gebaseerd is, kan ik er niet mee verder. De pythons die worden aangehaald om daar vervolgens de floats uit te destilleren, zijn op mijn systeem (grasshopper op rhino.compute van een derde partij) namelijk niet toegestaan.

Is er een methode om een op text gebaseerd antwoord te krijgen, of liefst zelfs gewoon direct de float-waarden?
Met die simpelere waarden kan ik namelijk wel wat. Als het json of xml is, krijg ik het er vast ook wel uit. Als het format maar voor het menselijk oog begrijpbaar is dan, zonder tussenkomst van een parser dus. Het lijkt me dat deze waarden ergens in een tussenstap de onderlegger voor de image formats moeten zijn geweest…

Graag zou ik ook een reactie krijgen van @PDOKbeheer! Misschien dat zij weten of dit nu of in de toekomst mogelijk is. We zouden er BIM-mend Nederland erg blij mee maken.

Met vriendelijke groet en bij voorbaat dank,

Thomas Broos

10

Goede morgen @ThomasBroos

Als je dit via een derde partij (Packhunt?, ShapDiver?) voor
elkaar wilt krijgen dan wordt het lastiger. Heron (via package manager) kan GEOtiffs uitlezen, maar wordt geloof ik ook niet ondersteund door deze partijen.
Je zult dan zelf een server omgeving moeten opzetten om rhino compute op te draaien.

1 like
11

@ThomasBroos op dit moment zijn de ondersteunde outputformats van de AHN WCS:

<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>
<wcs:formatSupported>image/vnd.jpeg-png</wcs:formatSupported>
<wcs:formatSupported>image/vnd.jpeg-png8</wcs:formatSupported>

Dat zijn dus allemaal image formats. De beste oplossing voor dit probleem zou denk ik zijn om een text based outputformat toe te laten voegen aan de WCS. door @PDOKbeheer. Ik denk dat het mogelijk moet zijn XYZ – ASCII Gridded XYZ als outputformaat toe te voegen. Dit maakt het consumeren van de WCS een stuk gebruiksvriendelijker en ook mogelijk in omgevingen waar geen image parsers beschikbaar zijn. Wat weer voor meer gebruik van de AHN data zorgt.

Dus feature request aan @PDOKbeheer?

3 likes
12

Ha @ebeere en @antonbakker , dank voor jullie antwoord. Dan zien we dat geloof ik hetzelfde. Ik heb inmiddels via het contactformulier van Pdok een featurerequest gedaan. Hopelijk is dat de goede weg.

Dank voor de hulp!

Thomas

1 like
13

@ThomasBroos bij dezen alvast even een reactie vanuit PDOK. We gaan even kijken wat een logische invulling is voor een text gebaseerd antwoord te krijgen op rasterdata. Ik ben het eens dat de oplossing in een coverage service zit (WCS) en niet zozeer in de viewservice (WM(T)S). Het formaat wat @antonbakker voorstelt lijkt me inderdaad het overwegen waard en nemen we mee in onze analyse.

1 like
14

Het is jammer dat van AHN alleen de ruwe data beschikbaar is in .laz scanvorm. .laz werkt ‘voor BIM’ zoveel veelzijdiger dan image-bestanden. Meeste BIM / CAD software het niet veel met geotif of ecw, behalve weergeven. Je kan offline bestanden in CloudCompare bewerken en omzetten naar xyz of .laz.

Wat ook wel intressant is is deze methode om de ruwe AHN data direct te downloaden binnen een bepaald gebied. Maar goed, daarvan kun je geen BIM maaiveld model produceren.

GitHub - HideBa/ahn_cli: Python CLI tool which helps you to download AHN point cloud data easily