mucking around with wcs requests

wcs.py

import requests
from lxml import etree

def hit(url, params):
    "Hit a URL, parse the xml"
    response = requests.get(url, params=params)
    print(response.url)
    if response.ok:
        print(response.headers)
        return etree.fromstring(response.content)
    else:
        print(f'Error: {response.reason}')

class WebCoverageService(object):
    
    namespaces = dict(
        wcs="http://www.opengis.net/wcs/1.1",
        ows="http://www.opengis.net/ows",
        owcs="http://www.opengis.net/wcs/1.1/ows",
        xsi="http://www.w3.org/2001/XMLSchema-instance",
        xlink="http://www.w3.org/1999/xlink",
    )
    
    def __init__(self, url):
        self.urls = {
            'capabilities': url,
            'describe': None,
            'getcoverage': None
        }
    
    @property
    def capabilities(self):
        try:
            assert self._capabilities is not None
            return self._capabilities
        except:
            # Initialize capabilities from the endpoint
            params = {
                'request': 'GetCapabilities', 
                'service': 'WCS'
            }
            self._capabilities = hit(self.urls['capabilities'], params)
            return self._capabilities
        
    @property
    def version(self):
        "The version to use (latest poss from service)"
        try:
            assert self._version is not None
            return self._version
        except:
            caps = self.capabilities
            self._versions = sorted(caps.xpath(
                'owcs:ServiceIdentification/owcs:ServiceTypeVersion/text()',
                namespaces=self.namespaces))
            self._version = self._versions[-1]
            return self._version
    
    @property
    def summary(self):
        "Return WCS summary"
        try:
            assert self._summary is not None
            return self._summary
        except:
            caps = self.capabilities
            self._summary = caps.xpath('./wcs:Contents/wcs:CoverageSummary',
                                       namespaces=self.namespaces)[0]
            return self._summary
        
    @property
    def bounds(self):
        "Return bounds"
        summary = self.summary
        lower = tuple(
            float(p) for p in 
            summary.xpath('ows:WGS84BoundingBox/ows:LowerCorner/text()', 
                          namespaces=self.namespaces)[0].split()
        )
        upper = tuple(
            float(p) for p in 
            summary.xpath('ows:WGS84BoundingBox/ows:UpperCorner/text()', 
                          namespaces=self.namespaces)[0].split()
        )
        return lower, upper
    
    @property
    def operations(self):
        try:
            assert self._operations is not None
            return self._operations
        except:
            self._operations = {}
            elements = self.capabilities.xpath(
                '.'
                '/owcs:OperationsMetadata'
                '/owcs:Operation',
                namespaces=self.namespaces
            )
            for operation in elements:
                name, info = self.parse_operation(operation)
                self._operations[name] = info
            return self._operations
    
    def parse_operation(self, operation):
        "Parse the XML info for a WCS operation"
        # Get URL to hit for operation
        operation_name = operation.attrib['name']
        info = {
            'url': operation.xpath(
                'owcs:DCP/owcs:HTTP/owcs:Get/@xlink:href', 
                namespaces=self.namespaces)[0]
        }
        self.urls[operation_name] = info['url']
                             
        # Pull out parameter info
        info['parameters'] =  {}
        parameters = operation.xpath('owcs:Parameter', namespaces=self.namespaces)
        for parameter in parameters:
            name = parameter.attrib['name'].lower()
            allowed_values = parameter.xpath('owcs:AllowedValues/owcs:Value/text()',
                                             namespaces=self.namespaces)
            info['parameters'][name] = allowed_values
        return operation_name, info
    
    @property
    def coverages(self):
        "Get a list of coverages"
        operations = self.operations
        return operations['DescribeCoverage']['parameters']['identifier']
    
    def describe(self, coverage=None, parameters=None):
        "Describe the coverage (may fail)"
        # Check we have a valid coverage
        if coverage is None:
            coverage = self.coverages[0]
        if coverage not in self.coverages:
            raise ValueError(f'Unknown coverage {coverage}, '
                              'allowed values are {self.coverages}')
        
        # Get description
        info = self.operations['DescribeCoverage']
        _params = {k: v[0] for k, v in info['parameters'].items()}
        _params['version'] = self.version
        if parameters is not None:
            _params.update(parameters)
        description = hit(self.urls['DescribeCoverage'], _params)
        return description
    
    def get_coverage(self, coverage=None, lower=None, upper=None, fmt=None, epsg=4326, **parameters):
        "Get data inside the given bounding box"
        # Check we have a valid coverage
        if coverage is None:
            coverage = self.coverages[0]
        if coverage not in self.coverages:
            raise ValueError(f'Unknown coverage {coverage}, '
                              'allowed values are {self.coverages}')
        
        # Create subset - default bounds are the bounds of the layer
        lower = lower or self.bounds[0]
        upper = upper or self.bounds[1]
        epsg = f'urn:ogc:def:crs:EPSG::{epsg}'
        bbox = f'{lower[0]},{lower[1]},{upper[0]},{upper[1]},{epsg}'
        
        # Default format is GeoTIFF, or first format if not available
        formats = self.operations['GetCoverage']['parameters']['format']
        if fmt is None:
            if 'image/GeoTIFF' in formats:
                fmt = 'image/GeoTIFF'
            else:
                fmt = formats[0]
        fmt = fmt or 'image/GeoTIFF'
        
        # All other parameters default to initial values
        info = self.operations['GetCoverage']
        _params = {k: v[0] for k, v in info['parameters'].items()}
        _params['version'] = self.version
        _params['store'] = False
        if parameters is not None:
            _params.update(parameters)
        _params['identifier'] = coverage
        _params['format'] = fmt
        _params['boundingbox'] = bbox
        
        # Hit endpoint 
        return hit(info['url'], _params)