source: cows/trunk/cows/service/imps/csmlbackend/wfs_csmllayer.py @ 6305

Subversion URL: http://proj.badc.rl.ac.uk/svn/ndg/cows/trunk/cows/service/imps/csmlbackend/wfs_csmllayer.py@6305
Revision 6305, 17.1 KB checked in by domlowe, 10 years ago (diff)

Changes to wfs backend to support point series subsetting in wfs interface

Line 
1# BSD Licence
2# Copyright (c) 2009, Science & Technology Facilities Council (STFC)
3# All rights reserved.
4#
5# See the LICENSE file in the source distribution of this software for
6# the full license text.
7
8"""
9implementation of ILayerMapper, IwfsLayer, IDimension, ILayerSlab interfaces, as defined in wfs_iface.py & wxs_iface.py
10
11"""
12import cows.service.imps.csmlbackend.wfs_csmlstoredqueries as CSMLQueries
13from cows.service.imps.csmlbackend.csmlcommon import CSMLLayerMapper, CSMLConnector
14from cows.service.wfs_iface import *
15from cows.model.storedquery import *
16import csml.parser
17from xml.etree import ElementTree as etree
18from shapely.geometry import Polygon, Point
19from cows.bbox_util import union
20from csml.csmllibs.csmlextra import listify, stringify, cleanString2
21
22import logging
23log = logging.getLogger(__name__)
24
25class CSMLwfsLayerMapper(CSMLLayerMapper):
26    """
27    Map keyword arguments to a collection of layers.
28    Supports the retrieval of sets of layers according to arbitrary
29    keyword/value pairs.
30    Implements  ILayerMapper (does it? TODO: check)
31   
32    WFS differs from WMS/WCS in that the 'layers' are feature types, not instances.
33    So the CSMLwfsLayerMapper map method returns both a map of feature types and a map
34    of feature instances as both are needed for the GetFeature method to work.
35   
36    """
37    def __init__(self):
38        super(CSMLwfsLayerMapper, self).__init__()
39        self.featureinstancecache={}
40        self.queryDescriptions=CSMLStoredQueries()
41   
42
43    def map(self, **kwargs):
44        """
45        Given csml.parser.Dataset object list the names of
46        all layers available.
47       
48        @return: A mapping of layer names to ILayer implementations.
49        @raise ValueError: If no layers are available for these keywords.
50        """
51        fileoruri=kwargs['fileoruri']
52        if fileoruri in self.layermapcache.keys():
53            #we've accessed this layer map before, get it from the cache dictionary
54            return self.layermapcache[fileoruri], self.featureinstancecache[fileoruri]
55         
56        ds = self.connector.getCsmlDoc(fileoruri)
57        try:
58            self.datasetName=ds.name.CONTENT
59        except AttributeError:
60            self.datasetName = 'CSML WFS Service'
61       
62        #The WFS differs from WMS & WCS in that the contents layer is a list of
63        #feature *types* not *instances*. However a record of instances is also
64        #needed to fulfil GetFeature requests:
65        featureset=CSMLFeatureSet() #holds feature instances               
66        layermap={} #feature types       
67        aggregatedBBoxes={}#for aggregations of bounding boxes.
68       
69
70        for feature in csml.csmllibs.csmlextra.listify(ds.featureCollection.featureMembers):
71            title, abstract=self.getInfo(feature)
72            featureset.featureinstances[feature.id]=CSMLFeatureInstance(title, abstract, feature)             
73        for id, instance in featureset.featureinstances.iteritems():
74            ftype=instance.featuretype #namespaced type e.g. 'csml:PointSeriesFeature'
75            if ftype not in layermap.keys():
76                layermap[ftype]=CSMLwfsLayer(ftype, instance.wgs84BBox)
77                #instantiate an aggregator to compare future bounding boxes with.
78                aggregatedBBoxes[ftype]= instance.wgs84BBox               
79            else:
80                #the featuretype has already been stored in the dictionary.
81                #but, the bounding box may need changing to accommodate this new feature instance.
82#                log.debug('Checking bbox for feature id: %s and title: %s'%(instance._feature.id, instance.title))
83                currentbbox=aggregatedBBoxes[ftype]
84                newbbox=union(currentbbox, instance.wgs84BBox)
85                aggregatedBBoxes[ftype]= newbbox
86                layermap[ftype]=CSMLwfsLayer(ftype, newbbox)
87           
88
89        if len(layermap) > 0:
90            #cache results
91            self.layermapcache[fileoruri]=layermap
92            self.featureinstancecache[fileoruri]=featureset
93            return layermap, featureset
94        else:
95            raise ValueError
96
97class CSMLStoredQueries(object):
98    """ Holds definitions of supported WFS stored queries and mappings to functional implementations"""
99    def __init__(self):
100        #self.queries is a dictionary of form: {id:(StoredQuery, func)}
101        #where func is the name of the implementation of the stored query (typically) in wfs_csmlstoredqueries     
102        self.queries={}
103       
104        #simple example query
105        qid='queryOne'       
106        qet=QueryExpressionText('xyzml:SomeFeatureType') 
107        pex1=ParameterExpression('arg1', 'xsd:string')
108        pex2=ParameterExpression('arg2', 'xsd:string')
109        query=StoredQuery(qid, title='test query', abstract = 'my test query', queryExpressionText = qet, parameter=[pex1, pex2])
110        self.queries[qid]=(query, CSMLQueries.query_one_func)
111     
112         
113        #mandatory get feature by id
114        qid='urn-x:wfs:StoredQueryId:ISO:GetFeatureById'
115        qet=QueryExpressionText('csml:AbstractFeature')
116        pex=ParameterExpression('id', 'xsd:anyURI')
117        query=StoredQuery(qid, title='GetFeatureById', abstract = 'Get any feature by id', queryExpressionText = qet, parameter=[pex])
118        self.queries[qid]=(query, CSMLQueries.query_getFeatureById)
119
120        #query to select features by phenomenon
121        qid='urn-x:wfs:StoredQueryId:badc.nerc.ac.uk:phenomenonQuery'
122        qet=QueryExpressionText('csml:AbstractFeature')#TODO, this could definitely differ for different features...
123        pex=ParameterExpression('phenomenon', 'xsd:string')
124        query=StoredQuery(qid, title='SelectFeaturesByPhenomenon', abstract = 'Select features based on their phenomenon type, e.g. "temperature"', queryExpressionText = qet, parameter=[pex])
125        self.queries[qid]=(query, CSMLQueries.query_getFeatureByPhenomenon)
126       
127        #query to extract a PointFeature from a PointSeriesFeature
128        qid='urn-x:wfs:StoredQueryId:badc.nerc.ac.uk:extractPointFromPointSeries'
129        qet=QueryExpressionText('csml:PointFeature')
130        pex1=ParameterExpression('featureid', 'xsd:anyURI')
131        pex2=ParameterExpression('timeinstance', 'gml:TimePositionUnion')
132        query=StoredQuery(qid, title='ExtractPointFromPointSeries', abstract = 'Extract a csml:PointFeature for a single time instance from a csml:PointSeriesFeature', queryExpressionText = qet, parameter=[pex1, pex2])
133        self.queries[qid]=(query, CSMLQueries.query_extractPointFromPointSeries)
134
135        #query to extract a PointSeriesFeature from a PointSeriesFeature
136        qid='urn-x:wfs:StoredQueryId:badc.nerc.ac.uk:extractPointSeriesFromPointSeries'
137        qet=QueryExpressionText('csml:PointSeriesFeature')
138        pex1=ParameterExpression('featureid', 'xsd:anyURI')
139        pex2=ParameterExpression('mintime', 'gml:TimePositionUnion')
140        pex3=ParameterExpression('maxtime', 'gml:TimePositionUnion')
141        query=StoredQuery(qid, title='ExtractPointSeriesFromPointSeries', abstract = 'Extract a csml:PointSeriesFeature for a range of times from a csml:PointSeriesFeature', queryExpressionText = qet, parameter=[pex1, pex2, pex3])
142        self.queries[qid]=(query, CSMLQueries.query_extractPointSeriesFromPointSeries)
143
144        #query to extract a GridSeriesFeature from a GridSeriesFeature
145        qid='urn-x:wfs:StoredQueryId:badc.nerc.ac.uk:extractGridSeriesFromGridSeries'
146        qet=QueryExpressionText('csml:GridSeriesFeature')
147        pex1=ParameterExpression('featureid', 'xsd:anyURI')
148        pex2=ParameterExpression('mintime', 'gml:TimePositionUnion')
149        pex3=ParameterExpression('maxtime', 'gml:TimePositionUnion')
150        pex4=ParameterExpression('bbox', 'xsd:string')
151        paramlist=[pex1,pex2,pex3, pex4]
152        query=StoredQuery(qid, title='ExtractGridSeriesFromGridSeries', abstract = 'Extract a csml:GridSeries from a csml:GridSeriesFeature', queryExpressionText = qet, parameter=paramlist)
153        self.queries[qid]=(query, CSMLQueries.query_extractGridSeriesFromGridSeries)
154
155        #query to extract a PointSeriesFeature from a GridSeriesFeature
156        qid='urn-x:wfs:StoredQueryId:badc.nerc.ac.uk:extractPointSeriesFromGridSeries'
157        qet=QueryExpressionText('csml:PointSeriesFeature')
158        pex1=ParameterExpression('featureid', 'xsd:anyURI')
159        pex2=ParameterExpression('mintime', 'gml:TimePositionUnion')
160        pex3=ParameterExpression('maxtime', 'gml:TimePositionUnion')
161        pex4=ParameterExpression('latitude', 'xsd:string')
162        pex5=ParameterExpression('longitude', 'xsd:string')
163        paramlist=[pex1,pex2,pex3, pex4, pex5]
164        query=StoredQuery(qid, title='ExtractPointSeriesFromGridSeries', abstract = 'Extract a csml:PointSeries from a csml:GridSeriesFeature', queryExpressionText = qet, parameter=paramlist)
165        self.queries[qid]=(query, CSMLQueries.query_extractPointSeriesFromGridSeries)
166       
167
168class CSMLFeatureSet(IFeatureSet):
169    """ A set of features available via a WFS. Supports querying methods as used by OGG filters """
170    def __init__(self):
171        self.featureinstances={}
172
173    def getFeatureList(self):
174        results=[]
175        for featureid,feature in self.featureinstances.iteritems():
176            results.append(feature)
177        return results
178
179    def getFeatureByGMLid(self, gmlid):
180        return self.featureinstances[gmlid]
181   
182    def _checkforOverlap(self, filterbbox, featurebbox):
183        """ Uses Shapely Polygons to calculate bounding box intersections """
184        log.debug('comparing against %s'%str(filterbbox))       
185        filterpolygon=Polygon(((filterbbox[0],filterbbox[1]),(filterbbox[0],filterbbox[3]),(filterbbox[2],filterbbox[3]),(filterbbox[2],filterbbox[1])))
186        featurepolygon=Polygon(((featurebbox[0],featurebbox[1]),(featurebbox[0],featurebbox[3]),(featurebbox[2],featurebbox[3]),(featurebbox[2],featurebbox[1])))       
187        log.debug(dir(filterpolygon))
188        log.debug('intersect result%s'%featurepolygon.intersects(filterpolygon))
189        return filterpolygon.intersects(featurepolygon)
190   
191    def _checkforPointinBox(self, filterbbox, location):
192        """ Uses Shapely Polygons to calculate bounding box intersections """
193        log.debug('comparing against %s'%str(filterbbox))       
194        filterpolygon=Polygon(((filterbbox[0],filterbbox[1]),(filterbbox[0],filterbbox[3]),(filterbbox[2],filterbbox[3]),(filterbbox[2],filterbbox[1])))
195        featurepoint=Point(float(location[0]), float(location[1]))
196        log.debug(featurepoint.within(filterpolygon))
197        log.debug('intersect result%s'%featurepoint.intersects(filterpolygon))
198        return featurepoint.intersects(filterpolygon)
199   
200    def getFeaturesByBBox(self,bboxtuple, srsname):         
201        log.debug('GET FEATURES BY BBOX')
202        result=[]
203        log.debug('Request BBOX: %s %s'%(bboxtuple,srsname))
204        for featureid,feature in self.featureinstances.iteritems():
205            if hasattr(feature._feature, 'boundedBy'):
206                log.debug('Checking bbox for feature %s: '%feature._feature.id)
207                lowercorner=feature._feature.boundedBy.lowerCorner.CONTENT.split()
208                uppercorner=feature._feature.boundedBy.upperCorner.CONTENT.split()
209                featurebbox=(float(lowercorner[0]), float(lowercorner[1]), float(uppercorner[0]), float(uppercorner[1]))               
210                if self._checkforOverlap(bboxtuple, featurebbox):
211                    result.append(feature) 
212            elif hasattr(feature._feature, 'location'):
213                log.debug('Checking location for feature %s: '%feature._feature.id)
214                log.debug(feature._feature.location)
215                location=feature._feature.location.CONTENT.split()
216                if self._checkforPointinBox(bboxtuple, location):
217                    result.append(feature)                       
218        return result
219
220    def getFeaturesByPropertyBetween(self, propertyname, minvalue, maxvalue):
221        log.debug('GET FEATURES BY PropertyBetween')
222        log.debug('propertyname, min, max: %s, %s, %s'%(propertyname, minvalue, maxvalue))
223        result = []
224        #Identifies times overlapping between filter and feature times
225        if propertyname=='csml:value/csml:PointSeriesCoverage/csml:pointSeriesDomain/csml:TimeSeries/csml:timePositionList':
226            for featureid, feature in self.featureinstances.iteritems():
227                featuretimes=feature._feature.value.pointSeriesDomain.timePositionList.CONTENT.split()
228                featureMinTime=featuretimes[0]
229                featureMaxTime=featuretimes[len(featuretimes)-1]
230                if csml.csmllibs.csmltime.compareTimes(featureMinTime, minvalue ,featureMaxTime) == 1:
231                    result.append(feature)
232                elif csml.csmllibs.csmltime.compareTimes(featureMinTime, maxvalue ,featureMaxTime) == 1:
233                    result.append(feature)   
234        return result
235       
236    def getFeaturesByPropertyEqualTo(self, propertyname, value):
237        log.debug('GET FEATURES BY PropertyEqualTo, value=%s'%value)
238        result=[]
239        value=value #value may be unicode
240        if propertyname == 'csml:parameter':
241            log.debug('filtering on csml:parameter')         
242            for featureid,feature in self.featureinstances.iteritems():
243                if feature._feature.parameter.getStandardName() == value:
244                    result.append(feature)     
245                elif feature._feature.parameter.getNonStandardName() == value:                   
246                    result.append(feature) 
247            return result
248   
249class CSMLFeatureInstance(IFeatureInstance):
250    def __init__(self, title, abstract, feature):
251        """ representing a CSML Feature Instance
252        @ivar title: The title of the feature instance
253        @ivar abstract: The abstract of the feature instance
254        @ivar feature: the csml feature instance object
255         """
256        self.title=title
257        self.abstract=abstract
258        self._feature=feature
259        self.featuretype='csml:'+self._feature.__class__.__name__
260        try:
261            bb= self._feature.getCSMLBoundingBox().getBox()
262            #convert 0 - 360 to -180, 180 as per common WMS convention
263            if abs(bb[2]-bb[0]) >= 359 and abs(bb[2]-bb[0]) < 361:
264                bb[0], bb[2]=-180, 180
265            self.wgs84BBox = bb
266        except AttributeError:
267            log.debug('there is a problem getting the bounding box for feature id %s'%self._feature.id)
268            self.wgs84BBox=()
269   
270    def _fixXlinks(self):
271        ''' replaces xlinks in feature domain with inline content'''
272        for att in ['gridSeriesDomain', 'pointDomain', 'profileSeriesDomain','sectionDomain','trajectoryDomain']:
273            if hasattr(self._feature.value, att):
274                domain=getattr(self._feature.value,att)
275                if hasattr(domain, 'coordTransformTable'):
276                    for ordinate in listify(domain.coordTransformTable.gridOrdinates):
277                        #insertedExtract is the 'hidden' resolved xlink data, expose inline
278                        if hasattr(ordinate.coordAxisValues, 'insertedExtract'):
279                            ordinate.coordAxisValues.CONTENT=cleanString2(str(ordinate.coordAxisValues.insertedExtract.getData()[0].tolist()))                           
280   
281    def toGML(self):
282        """ Create a GML (CSML) representation of the feature """
283        nonamespaceFType=self.featuretype.split(':')[1] #remove the csml: prefix
284        qualifiedFeatureType='{http://ndg.nerc.ac.uk/csml}' + nonamespaceFType
285        emptyelem=etree.Element(qualifiedFeatureType)
286        #modify self._feature so that any references to xlinks are replaced with inline content
287        self._fixXlinks()       
288        csmlelem=self._feature.toXML(emptyelem)
289        return etree.tostring(csmlelem)
290       
291       
292     
293class CSMLwfsLayer(IwfsLayer):
294    """ representing a WFS FeatureType (termed layer here). Implements IwfsLayer
295    @ivar featuretype: The namespaced name of the feature type, e.g. csml:PointSeriesFeature
296   
297    """
298    def __init__(self, featuretype, wgs84bb):
299        self.type=featuretype
300        #Have to hard code some definitions for CSML feature types to
301        #use in the capabilities document as they don't exist anywhere else.
302        #Hardcoding is okay as this is a CSML specific interface, so will only ever deal
303        #with CSML feature types.
304        #TODO: However, might be better to move to some sort of schema.cfg file?
305        self.wgs84BBox=wgs84bb     
306        self.outputformats=['text/xml: subtype=csml/2.0']   
307        if self.type=='csml:GridSeriesFeature':
308            self.title='GridSeriesFeature as defined in Climate Science Modelling Language.'
309            self.abstract='The CSML GridSeriesFeature is used to represent 4D gridded data such as atmospheric model output.'
310            self.keywords=['Grid', 'Climate', 'CSML']
311        elif self.type=='csml:PointSeriesFeature':
312            self.title='PointSeriesFeature as defined in Climate Science Modelling Language.'
313            self.abstract='The CSML PointSeriesFeature represents a time series of measurements at a single point in space.'
314            self.keywords=['Point', 'Timeseries', 'Climate', 'CSML']
315        #TODO: definitions for all feature types.
316       
317       
318
319
Note: See TracBrowser for help on using the repository browser.