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.

1 like