Process for averaging over a polygon

flyingpigeon/processes/__init__.py

@@ -1,5 +1,6 @@
 
 # from .wps_say_hello import SayHello
+from .wps_average_wfs_polygon import AverageWFSPolygonProcess
 from .wps_subset_wfs_polygon import SubsetWFSPolygonProcess
 from .wps_subset_bbox import SubsetBboxProcess
 from .wps_subset_continents import SubsetcontinentProcess
@@ -17,6 +18,7 @@
 
 processes = [
     # SayHello(),
+    AverageWFSPolygonProcess(),
     SubsetWFSPolygonProcess(),
     SubsetBboxProcess(),
     SubsetcontinentProcess(),

flyingpigeon/processes/subset_base.py

@@ -46,6 +46,42 @@
                          as_reference=True,
                          supported_formats=[FORMATS.META4])
 
+typename = LiteralInput('typename',
+                        'TypeName',
+                        abstract='Name of the layer in GeoServer.',
+                        data_type='string',
+                        min_occurs=0,
+                        max_occurs=1)
+
+featureids = LiteralInput('featureids',
+                          'Feature Ids',
+                          abstract='fid(s) of the feature in the layer.',
+                          data_type='string',
+                          min_occurs=0,
+                          max_occurs=1000)
+
+geoserver = LiteralInput('geoserver',
+                         'Geoserver',
+                         abstract=('Typically of the form '
+                                   'http://host:port/geoserver/wfs'),
+                         data_type='string',
+                         min_occurs=0)
+
+mosaic = LiteralInput('mosaic',
+                      'Union of Feature Ids',
+                      abstract=('If True, selected regions will be '
+                                'merged into a single geometry.'),
+                      data_type='boolean',
+                      min_occurs=0,
+                      default=False)
+
+shape = ComplexInput('shape', 'Vector Shape',
+                     abstract='An ESRI Shapefile, GML, JSON, GeoJSON, or single layer GeoPackage.'
+                              ' The ESRI Shapefile must be zipped and contain the .shp, .shx, and .dbf.',
+                     min_occurs=1,
+                     max_occurs=1,
+                     supported_formats=[FORMATS.GEOJSON, FORMATS.GML, FORMATS.JSON, FORMATS.SHP])
+
 
 def get_feature(url, typename, features):
     """Return geometry from WFS server."""
@@ -99,7 +135,7 @@ def make_geoms(feature, mosaic=False):
     crs = ocgis.CoordinateReferenceSystem(epsg=crs_code)
     geoms = [
         {'geom': shape(f['geometry']), 'crs': crs,
-         'properties': f['properties']}
+         'properties': f['properties'] if not isinstance(f['properties'], list) else "List not supported by ocgis"}
         for f in feature['features']]
 
     if mosaic:
@@ -171,13 +207,19 @@ def parse_feature(self, request):
                     geoserver, typename, featureids, e)
                 raise Exception(msg) from e
 
+            # Remove properties because it crashes ocgis
+            for geom in geoms.values():
+                if isinstance(geom, dict):
+                    geom.pop("properties")
+
+
             if mosaic:
                 return {'mosaic': geoms}
             else:
                 return dict(zip(featureids, geoms))
 
-        else:
-            return {}
+        elif 'shape' in request.inputs:
+            return {'_shp_': request.inputs['shape'][0].file}
 
     def parse_daterange(self, request):
         """Return [start, end] or None."""

flyingpigeon/processes/wps_average_polygon.py

@@ -0,0 +1,24 @@
+from pywps import Process
+from .subset_base import Subsetter, resource, variable, start, end, output, metalink, shape, mosaic
+
+
+class AveragePolygonProcess(Process, Subsetter):
+    """Subset a NetCDF file using WFS geometry."""
+
+    def __init__(self):
+        inputs = [resource, shape, mosaic, start, end, variable]
+        outputs = [output, metalink]
+
+        super().__init__(
+            self._handler,
+            identifier='average-polygon',
+            title='Average over polygon',
+            version='0.1',
+            abstract=('Return the average of the data for which grid cells intersect the '
+                      'selected polygon for each input dataset as well as'
+                      'the time range selected.'),
+            inputs=inputs,
+            outputs=outputs,
+            status_supported=True,
+            store_supported=True,
+        )

flyingpigeon/processes/wps_average_wfs_polygon.py

@@ -0,0 +1,74 @@
+import tempfile
+from pathlib import Path
+
+from pywps import Process, FORMATS
+from pywps.inout.outputs import MetaFile, MetaLink4
+
+from .subset_base import Subsetter, resource, variable, start, end, output, metalink, typename, \
+    featureids, geoserver, mosaic
+
+import ocgis
+import ocgis.exc
+
+
+class AverageWFSPolygonProcess(Process, Subsetter):
+    """Subset a NetCDF file using WFS geometry."""
+
+    def __init__(self):
+        inputs = [resource, typename, featureids, geoserver, mosaic, start, end, variable]
+        outputs = [output, metalink]
+
+        super(AverageWFSPolygonProcess, self).__init__(
+            self._handler,
+            identifier='average-wfs-polygon',
+            title='Average over polygon',
+            version='0.1',
+            abstract=('Return the average of the data for which grid cells intersect the '
+                      'selected polygon for each input dataset as well as'
+                      'the time range selected.'),
+            inputs=inputs,
+            outputs=outputs,
+            status_supported=True,
+            store_supported=True,
+        )
+
+    def _handler(self, request, response):
+
+        # Gather geometries, aggregate if mosaic is True.
+        geoms = self.parse_feature(request)
+        dr = self.parse_daterange(request)
+
+        ml = MetaLink4('subset', workdir=self.workdir)
+
+        if geoms.get('_shp_', None) is not None:
+            geom = geoms['_shp_']
+        else:
+            geom = [g for g in geoms.values()]
+
+        for res in self.parse_resources(request):
+            variables = self.parse_variable(request, res)
+            rd = ocgis.RequestDataset(res, variables)
+            prefix = Path(res).stem
+
+            try:
+                ops = ocgis.OcgOperations(
+                    dataset=rd, geom=geom,
+                    spatial_operation='clip', aggregate=True,
+                    time_range=dr, output_format='nc',
+                    interpolate_spatial_bounds=True,
+                    prefix=prefix, dir_output=tempfile.mkdtemp(dir=self.workdir))
+
+                out = ops.execute()
+
+                mf = MetaFile(prefix, fmt=FORMATS.NETCDF)
+                mf.file = out
+                ml.append(mf)
+
+            except ocgis.exc.ExtentError:
+                continue
+
+        response.outputs['output'].file = ml.files[0].file
+        response.outputs['metalink'].data = ml.xml
+        response.update_status("Completed", 100)
+
+        return response

flyingpigeon/processes/wps_subset_wfs_polygon.py

@@ -4,7 +4,8 @@
 from pywps import Process, LiteralInput, FORMATS
 from pywps.inout.outputs import MetaFile, MetaLink4
 
-from .subset_base import Subsetter, resource, variable, start, end, output, metalink
+from .subset_base import Subsetter, resource, variable, start, end, output, metalink, typename, \
+    featureids, geoserver, mosaic
 
 import ocgis
 import ocgis.exc
@@ -14,42 +15,14 @@ class SubsetWFSPolygonProcess(Process, Subsetter):
     """Subset a NetCDF file using WFS geometry."""
 
     def __init__(self):
-        inputs = [
-            resource,
-            LiteralInput('typename',
-                         'TypeName',
-                         abstract='Name of the layer in GeoServer.',
-                         data_type='string',
-                         min_occurs=0,
-                         max_occurs=1),
-            LiteralInput('featureids',
-                         'Feature Ids',
-                         abstract='fid(s) of the feature in the layer.',
-                         data_type='string',
-                         min_occurs=0,
-                         max_occurs=1000),
-            LiteralInput('geoserver',
-                         'Geoserver',
-                         abstract=('Typically of the form '
-                                   'http://host:port/geoserver/wfs'),
-                         data_type='string',
-                         min_occurs=0),
-            LiteralInput('mosaic',
-                         'Union of Feature Ids',
-                         abstract=('If True, selected regions will be '
-                                   'merged into a single geometry.'),
-                         data_type='boolean',
-                         min_occurs=0,
-                         default=False),
-            start, end, variable]
-
+        inputs = [resource, typename, featureids, geoserver, mosaic, start, end, variable]
         outputs = [output, metalink]
 
         super(SubsetWFSPolygonProcess, self).__init__(
             self._handler,
             identifier='subset-wfs-polygon',
             title='Subset',
-            version='0.2',
+            version='0.3',
             abstract=('Return the data for which grid cells intersect the '
                       'selected polygon for each input dataset as well as'
                       'the time range selected.'),
@@ -79,7 +52,7 @@ def _handler(self, request, response):
                 try:
                     ops = ocgis.OcgOperations(
                         dataset=rd, geom=geom['geom'],
-                        spatial_operation='clip', aggregate=True,
+                        spatial_operation='clip', aggregate=False,
                         time_range=dr, output_format='nc',
                         interpolate_spatial_bounds=True,
                         prefix=prefix, dir_output=tempfile.mkdtemp(dir=self.workdir))

setup.cfg

@@ -22,7 +22,7 @@ replace = Version="{new_version}"
 addopts =
 	--strict
 	--tb=native
-	tests/
+
 python_files = test_*.py
 markers =
 	online: mark test to need internet connection

tests/test_average_wfs_polygon.py

@@ -0,0 +1,39 @@
+from pywps import Service
+from pywps.tests import client_for, assert_response_success
+
+from .common import get_output
+from flyingpigeon.processes.wps_average_wfs_polygon import AverageWFSPolygonProcess
+
+
+def test_wps_subset():
+    client = client_for(Service(processes=[AverageWFSPolygonProcess()]))
+
+    datainputs = "resource=@xlink:href={nc};" \
+                 "typename={typename};" \
+                 "geoserver={geoserver};" \
+                 "featureids={fid1};" \
+                 "featureids={fid2};" \
+                 "mosaic={mosaic}"
+
+    inputs = dict(nc="https://pavics.ouranos.ca/twitcher/ows/proxy/thredds/dodsC/birdhouse/nrcan"
+                     "/nrcan_northamerica_monthly"
+                     "/tasmax/nrcan_northamerica_monthly_2015_tasmax.nc",
+                  typename="public:USGS_HydroBASINS_lake_na_lev12",
+                  geoserver="https://pavics.ouranos.ca/geoserver/wfs",
+                  fid1="USGS_HydroBASINS_lake_na_lev12.67061",
+                  fid2="USGS_HydroBASINS_lake_na_lev12.67088",
+                  mosaic=False)
+
+    resp = client.get(
+        "?service=WPS&request=Execute&version=1.0.0&identifier=average-wfs-polygon&datainputs={}".format(
+            datainputs.format(**inputs)))
+    assert_response_success(resp)
+
+    out = get_output(resp.xml)
+    assert 'metalink' in out
+
+    inputs['mosaic'] = True
+    resp = client.get(
+        "?service=WPS&request=Execute&version=1.0.0&identifier=average-wfs-polygon&datainputs={}".format(
+            datainputs.format(**inputs)))
+    assert_response_success(resp)