Python script to create WCS GetCoverage request with OWSLib #owslib #python #wcs #ahn #raster

wcs-getcoverage-request.py

#!/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()