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

Datasets AHN
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