Een bbox om een specifiek (x,y)-punt/coordinaat voor het uitlezen van informatie uit AHN

Datasets AHN
#1

Hallo,

Ik probeer vanuit een WCS AHN3 (of indien ook aanwezig data vanuit AHN4 het liefst), een DSM met resolutie van 0.5m uit te lezen en hier een bbox (boundary box) om te creeeren. Ik heb al het een en ander gelezen, maar ik loop vast nu op de volgende punten:

• Wat is de juiste link voor het verkrijgen van AHN data? Ik krijg bij dit URL in de code, die voorheen wel gelukt is een paar jaar terug (http://geodata.nationaalgeoregister.nl/ahn3/ows), de volgende foutmelding: {“timestamp”:1683626459532,“status”:404,“error”:“Not Found”,“message”:“Not Found”,“Request”:"",“Server”:null,“Dataset”:null,“Service”:null}
• Zou de onderstaande code genoeg input geven om zo een correcte ‘bbox’ te creeeren om een coordinaat heen?
Zie de onderstaande code. Mijn vraag hier is met name of ik hieruit van dus een gehele website (met werkende link hopelijk), een stukje van de toegewijde .TIF bestand terug kan krijgen om zo verder analyse te doen.

# Import required packages
from PIL import Image
import numpy as np
from owslib.wcs import WebCoverageService
import copy
import json

url = ‘http://geodata.nationaalgeoregister.nl/ahn3/ows’
wcs = WebCoverageService(url, version=‘1.0.0’)

# Input for function
x = float(lng)
y = float(lat)
radius = 50  # in meters
min_radius = 2
height_offset = 1.5
slices = 360

conv = 0.000011  # conversion from meters to degrees
subset_bbox = x-(radius*conv), y-(radius*conv), x+(radius*conv), y+(radius*conv)
grid_dist = 0.5

output = wcs.getCoverage(identifier='ahn3_05m_DSM',
                         format='GEOTIFF_FLOAT32',
                         crs='EPSG:4326',
                         response_crs='EPSG:4326',
                         bbox=subset_bbox,
                         WIDTH=radius*2,
                         HEIGHT=radius*2)

with open('test.tiff', 'wb') as f:
    f.write(output.read())

data = Image.open('test.tiff')

Veel dank alvast.

1 like
#2

Hallo sbroek,

Welkom! Ik heb je code zo gauw niet bekeken of geprobeerd, maar de url voor de AHN3 is inmiddels gewijzigd:

https://service.pdok.nl/rws/ahn3/wcs/v1_0

Dus vandaar de 404 die je terug krijgt.

#3

Welkom @sbroek!

Wat betreft je vragen:

  1. Juiste link: zie indd antwoord @sbjager. Als PDOK publiceren we de up-to-date service url’s via het Nationaalgeorgister en onze website.
  2. Zie hieronder een code voorbeeld om een werkend WCS GetCoverage te doen met OWSLib in Python. Het request en de output is in RD (EPSG:28992) om het simpel te houden. Als het goed is moet je ook een geotiff in EPSG:4326 kunnen terugkrijgen van de WCS, maar ik heb daar niet zo snel een voorbeeld van liggen. Code voorbeeld is op basis van de code van de QGIS Pdokservicesplugin. Lees overigens ook dit antwoord over het uitlijnen van je request met het coverage zoals dat aangeboden wordt in WCS.
#!/usr/bin/env python3
from owslib.wcs import WebCoverageService

def main():
    x = 191313.0  
    y = 440283.3 # coordinaten in RD/EPSG:28992 van een locatie in NL
    origin = [10000,618750] # origin van ahn3 coverage in RD/EPSG:28992
    cell_size = 0.5 # cell size van ahn3 coverage
    bbox_size = 50 # in meters - dit is formaat van coverage dat opgehaald wordt
    bbox_size_pixels=bbox_size/cell_size
    wcs_url = "https://service.pdok.nl/rws/ahn3/wcs/v1_0?request=GetCapabilities&service=WCS"
    
    x_lower_bound = origin[0] + (((x - origin[0]) // cell_size) * cell_size) # n.b. // -> floor operator  https://www.askpython.com/python/python-floor-division-double-slash-operator
    x_upper_bound = x_lower_bound + (bbox_size_pixels * cell_size)
    y_lower_bound = origin[1] + (((y - origin[1]) // cell_size) * cell_size)
    y_upper_bound = y_lower_bound + (bbox_size_pixels * cell_size)
   
    wcs = WebCoverageService(wcs_url, version='2.0.1')       
    output = wcs.getCoverage(
                                identifier=['ahn3_05m_dtm'],
                                format='GEOTIFF_FLOAT32',
                                crs='EPSG:28992',
                                subsets = [('X', x_lower_bound, x_upper_bound), ('Y', y_lower_bound, y_upper_bound)],
                                width=bbox_size_pixels,
                                height=bbox_size_pixels,
                            )
    with open('test.tiff', 'wb') as f:
        f.write(output.read())
        print(f"ouput written to test.tiff")

if __name__ == "__main__":
    main()
2 likes
#4

Beste @sbjager,

Ik zag op Changelog - Esri Nederland Content dat de data van AHN4 beschikbaar is. Kan dit kloppen? En zo ja, wat is de juiste URL voor de vernieuwde AHN4 data?

Ik hoor het graag!

Groet,
Sanne

#5

Beste @antonbakker,

Ik ben momenteel bezig met de implementatie van AHN4 data via de wcs url (https://service.pdok.nl/rws/ahn/wcs/v1_0?request=getcapabilities&service=wcs). Echter om deze data uit te lezen, heb ik nieuwe inputs nodig in mijn wcs.getCoverage code gedeelte.
output = wcs.getCoverage(identifier=‘ahn4_05m_DSM’,
format=‘GEOTIFF_FLOAT32’,
crs=‘EPSG:4326’,
response_crs=‘EPSG:4326’,

Wat zijn deze voor de AHN4 data? met name voor de ‘identifier’ en ‘format’ - ik hoor het graag!

Groet, Sanne

#6

Hoi @sbroek, als je in het capabilitities document van de WCS kijkt kan je zien dat de coverages die aangeboden worden de volgende identifiers hebben:

  • dsm_05m
  • dtm_05m

En dat de format parameter voor Geotiff output image/tiff is. Dus dan wordt de wcs.getCoverage call:

output = wcs.getCoverage(
                                identifier=['dtm_05m'],
                                format='image/tiff',
                                crs='EPSG:28992',
                                subsets=[('X', x_lower_bound, x_upper_bound), ('Y', y_lower_bound, y_upper_bound)],
                                width=bbox_size_pixels,
                                height=bbox_size_pixels,
                            )
Nieuwe versie AHN beschikbaar via PDOK
#7

Hi @antonbakker ,

Bedankt voor je reactie! Nu heb ik het volgende geprobeerd, echter krijg ik hierbij de volgende error:
output = wcs.getCoverage(identifier=[‘dsm_05m’],

File C:\Program Files\Python311\Lib\site-packages\owslib\coverage\wcs201.py:156 in getCoverage
if log.isEnabledFor(logging.DEBUG):

AttributeError: module ‘logging’ has no attribute ‘isEnabledFor’

Kan het zijn dat dit iets te doen heeft met de versie of dat het nog niet is geactiveerd voor deze url/AHN4 bestand?

Ik hoor het graag!

Groet,
Sanne

#8

Hoi @sbroek, ik vermoed dat de foutmelding iets te maken heeft met hoe Python en de OWSLib dependency geïnstalleerd is op jouw systeem. Ik denk dat daar iets fout zit, maar kan op basis van deze info niet zeggen wat. De foutmelding wordt veroorzaakt door OWSLib die de isEnabledFor methode aanroept van de logging module, die op jouw machine niet bestaat - die wel door OWSLib verwacht wordt.

Zie bijvoorbeeld deze stackoverflow vraag met een vergelijkbare foutmelding (maar niet identiek).

De foutmelding heeft iig niets te maken met de nieuwe AHN versie, want “works on my machine”.

#9

Hi Anton,

Het is eventjes geleden, maar echter krijg ik de AHN4 URL nog steeds niet aan de praat. Ik gebruik hiervoor Spyder 5.4.3 (Python 3.11.1). Na veel onderzoek, zou het daar kennelijk niet aan moeten liggen. Heb jij toevallig een voorbeeld code met een bepaalde coordinaat (i.e. 52.006884, 4.357335) waarvoor die het bij jou wel doet?

Ik hoor het heel graag! Het zou mijn afstudeerproject enorm helpen.

Bedankt,

Groet,
Sanne

1 like