AHN WCS bevragen in EPSG:4326

Datasets AHN
,
#1

Beste @antonbakker, @Jeroen_D,

Ik heb momenteel nog steeds een ‘error’ op het gebied van data uitlezen bij AHN4. Ik heb meedere malen in de materie en XML bestand gedoken, maar begrijp niet waar het verschil zit. Het WCS gedeelte van mijn code is:

url ='https://service.pdok.nl/rws/ahn/wcs/v1_0?request=getcapabilities&service=wcs'
    wcs = WebCoverageService(url, version='2.0.1')
  
    
    # Input for function
    x = float(longitude)
    y = float(latitude)
    radius = 50  # in meters

    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=['dsm_05m'],
                                format='image/tiff',          # Format of the output (image or TIFF file)
                                crs='EPSG:28992',                   # Coordinate reference system         
                                response_crs='EPSG:4326',          
                                bbox=subset_bbox,
                                WIDTH=radius*2,
                                HEIGHT=radius*2)
    
    with open('test4.tif', 'wb') as f:
        f.write(output.read())

    data = Image.open('test4.tif')```

In precies zou de werking erachter moeten kloppen volgens mij, maar ik krijg daarbij steeds de volgende error:

ServiceException: <?xml version="1.0" encoding="UTF-8"?>
<ows:ExceptionReport xmlns:ows=“http://www.opengis.net/ows/2.0” xmlns:xsi=“http://www.w3.org/2001/XMLSchema-instance” version=“2.0.1” xml:lang=“en-US” xsi:schemaLocation=“http://www.opengis.net/ows/2.0 http://schemas.opengis.net/ows/2.0/owsExceptionReport.xsd”>
<ows:Exception exceptionCode=“InvalidParameterValue” locator=“size”>
ows:ExceptionTextmsWCSGetCoverage20(): WCS server error. Raster size out of range, width and height of resulting coverage must be no more than MAXSIZE=4000.</ows:ExceptionText>
</ows:Exception>
</ows:ExceptionReport>```

Iemand die mij hierbij kan helpen? Bij de ‘oude’ url voor AHN3 lukt het wel gek genoeg (met aanpassingen van de identifier etc. natuurlijk).

Ik hoor het heel graag! Of eventueel ook andere tips.

Alvast bedankt.

Groet,
Sanne

1 like
Nieuwe versie AHN beschikbaar via PDOK
#2

Hoi @sbroek, ik zie zo direct twee problemen met je code:

  1. Afgaande op bovenstaande code bevat subset_bbox een boundingbox in degrees (EPSG:4326), maar in je aanroep getCoverage specificeer je EPSG:28992.

  2. response_crs keywordt arg bestaat niet in getCoverage zie de OWSLIB broncode. Je kan wel extra keyword args meegeven die aan het request worden meegegeven als query parameters. Dan zou je OUTPUTCRS moeten gebruiken.

Ik heb een voorbeeld script geschreven met debug output van OWSLIB zodat eventuele problemen eenvoudiger te diagnosticeren zijn. Verder was het nodig om de spatial subset mee te geven met subsets = [('X', x_min, x_max), ('Y', y_min, y_max)], in plaats van het bbox kw arguments. Is mij niet helemaal duidelijk wat het bbox argument dan doet in getCoverage, als je er geen spatial subset mee kan opvragen…

#!/usr/bin/env python3
from owslib.wcs import WebCoverageService
import logging

WCS_URL ='https://service.pdok.nl/rws/ahn/wcs/v1_0?request=getcapabilities&service=wcs'

# enable logging of OWSLib -- see https://github.com/geopython/OWSLib#logging
owslib_log = logging.getLogger('owslib')
ch = logging.StreamHandler()
ch.setLevel(logging.DEBUG)
ch.setFormatter(logging.Formatter('%(asctime)s - %(name)s - %(levelname)s - %(message)s'))
owslib_log.addHandler(ch)
owslib_log.setLevel(logging.DEBUG)

def get_coverage(x, y, output_filepath, bbox_size=100):
    """WCS GetCoverage function by x, y coordinates and bbox_size in RD (EPSG:28992)
    Arguments:
        x {float} -- x coordinate in EPSG:28992
        y {float} -- x coordinate in EPSG:28992
        output_filepath {string} -- output file path output coverage
    Keyword Arguments:
        bbox_size {integer} -- bounding box size in meters, constructed around (x,y)
        output_coverage_size {integer} -- output coverage size in pixels 
    """

    x_min=x-(bbox_size/2)
    x_max=x+(bbox_size/2)
    y_min=y-(bbox_size/2)
    y_max=y+(bbox_size/2)

    wcs = WebCoverageService(WCS_URL, version='2.0.1')    
    output = wcs.getCoverage(identifier=['dsm_05m'], # see https://github.com/geopython/OWSLib/blob/9530e923697badd5e5962e8f0243f3861150edd0/owslib/coverage/wcs201.py#L117C1-L117C1 for all function arguments (OWSLIB 0.28)
                             format='image/tiff',
                             crs='EPSG:28992',
                             subsets = [('X', x_min, x_max), ('Y', y_min, y_max)],
                             OUTPUTCRS='EPSG:4326'
                            )
    
    with open(output_filepath, 'wb') as f:
        f.write(output.read())
        print(f"output coverage saved in {output_filepath}")

get_coverage(190455.3,445219.9, "test_4326.tif", bbox_size=200)
1 like
#3

Beste @antonbakker,

Enorm bedankt voor de hulp. Het is mij nu duidelijker waar de inputs uit bestaan en hoe het samen komt. Ik ben momenteel bezig met het verwerken van mijn data, echter zoals je al aangeeft is de EPSG:28992 gesprecifieerd in de getCoverage. Ik heb geprobeerd deze om te zetten (in de CRS) naar EPSG:4326 omdat ik hierbij coordinaten in DD kan invoegen), echter zonder succes. Ik krijg dat steeds dezelfde foutmelding…
Is dit eenvoudig te ‘converten’ binnen de WCS GetCoverage request functie of is dit wel de standaard waarin coordinaten gegeven moeten worden?

Ik hoor het graag. Bedankt,

Groet,
Sanne

#4

Voor het bevragen van de WCS in EPSG:4326 is het nodig om de query param SUBSETTINGCRS mee te geven; zie het volgende code voorbeeld:

#!/usr/bin/env python3
from owslib.wcs import WebCoverageService
import logging

WCS_URL ='https://service.pdok.nl/rws/ahn/wcs/v1_0?request=getcapabilities&service=wcs'

# enable logging of OWSLib -- see https://github.com/geopython/OWSLib#logging
owslib_log = logging.getLogger('owslib')
ch = logging.StreamHandler()
ch.setLevel(logging.DEBUG)
ch.setFormatter(logging.Formatter('%(asctime)s - %(name)s - %(levelname)s - %(message)s'))
owslib_log.addHandler(ch)
owslib_log.setLevel(logging.DEBUG)

def get_coverage_lonlat(lon,lat, output_filepath, bbox_size=(0.0015,0.0009)): # bbox_size corresponds to approximately 100x100m within Netherlands (European) land boundaries
    """WCS GetCoverage function by x, y coordinates and bbox_size in degrees (EPSG:4326)
    Arguments:
        lon {float} -- lon coordinate in EPSG:4326, centre point boundingbox
        lat {float} -- lat coordinate in EPSG:4326, centre point boundingbox        
        output_filepath {string} -- output file path output coverage
    Keyword Arguments:
        bbox_size {(float,float)} -- bounding box size in degrees, bbox_size[0]; lon, bbox_size[1];lat. Bounding box constructed around (lon,lat) - 1m ~ 0.000146 degrees
    """
 
    lon_min=lon-(bbox_size[0]/2)
    lon_max=lon+(bbox_size[0]/2)
    lat_min=lat-(bbox_size[1]/2)
    lat_max=lat+(bbox_size[1]/2) 

    print(f"{lon_min},{lat_min},{lon_max},{lat_max}")
    wcs = WebCoverageService(WCS_URL, version='2.0.1')
    crs='EPSG:4326'

    # see https://github.com/geopython/OWSLib/blob/9530e923697badd5e5962e8f0243f3861150edd0/owslib/coverage/wcs201.py#L117C1-L117C1 for all function arguments (OWSLIB 0.28)
    # last two keyword arguments are passed through to wcs getcoverage request as query params
    # see mapserver docs for more info on wcs request parameters: https://mapserver.org/ogc/wcs_server.html#wcs-2-0-kvp-request-parameters
    output = wcs.getCoverage(identifier=['dsm_05m'], 
                             format='image/tiff',
                             crs=crs,
                             subsets = [('X', lon_min, lon_max), ('Y', lat_min, lat_max)],
                             OUTPUTCRS=crs,
                             SUBSETTINGCRS=crs
                            )
    with open(output_filepath, 'wb') as f:
        f.write(output.read())
        print(f"output coverage saved in {output_filepath}")

get_coverage_lonlat(5.926126,52.009562, "test_4326.tif")