source: cows/trunk/cows/service/imps/csmlbackend/wcs_csmllayer.py @ 6398

Subversion URL: http://proj.badc.rl.ac.uk/svn/ndg/cows/trunk/cows/service/imps/csmlbackend/wcs_csmllayer.py@6398
Revision 6398, 10.3 KB checked in by pnorton, 11 years ago (diff)

Modified the WCS code to work with the csml files in folders in the same way as the WMS code.

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, IwmsLayer, IwmsDimension, ILayerSlab interfaces, as defined in wms_iface.py & wxs_iface.py
10
11"""
12import os, string
13import csml
14try:
15    import cdms2 as cdms
16except:
17    import cdms
18from copy import copy
19
20import logging
21log = logging.getLogger(__name__)
22
23from cows.model.wcs import WcsDatasetSummary
24from cows.model.common import BoundingBox, WGS84BoundingBox
25from cows.service.imps.csmlbackend.csmlcommon import CSMLLayerMapper, CSMLConnector, extractToNetCDF
26
27
28class CSMLwcsCoverageMapper(CSMLLayerMapper): 
29    """
30    Map keyword arguments to a collection of layers (coverages)
31    Supports the retrieval of coverages according to arbitrary
32    keyword/value pairs. 
33    """
34    def __init__(self):
35        super(CSMLwcsCoverageMapper, self).__init__()
36       
37   
38    def getInfo(self, feature):
39        ''' given a csml feature, return info about the layer/feature
40        @return:   title, abstract, units, crss '''
41
42        title, abstract = super(CSMLwcsCoverageMapper, self).getInfo(feature)
43        units=feature.getDomainUnits()
44        tmpunits=copy(units)
45        tmpunits.reverse()
46        domain = feature.getDomain()
47        tax= feature.getTimeAxis()
48       
49        try:
50            timepositions=domain[tax]
51            timelimits=[domain[tax][0]+'Z',domain[tax][len(domain[tax])-1]+'Z']
52        except: #no time available, return empty strings
53            timepositions=['']
54            timelimits=['',''] 
55           
56        crs=feature.getNativeCRS()
57        log.debug('crs %s'%crs)
58        crss=[self._crscat.getCRS(crs).twoD] 
59        if 'EPSG:4326' in crss:
60            crss.append('CRS:84')
61            crss.append('WGS84')
62        log.debug('crss %s'%crss)
63        #build definitions of any Z axes such as air_pressure:
64        axisDescriptions=[]       
65        lon=feature.getLongitudeAxis()
66        lat=feature.getLatitudeAxis()
67        t=feature.getTimeAxis()
68       
69        if None in [lon, lat, t]:
70            #TODO need to return a suitable wcs error.
71            log.debug('warning, could not get correct axis info')
72            #best guess!
73            if t is None:
74                t='time'
75            if lon is None:
76                lon = 'longitude'
77            if lat is None:
78                lat = 'latitude'
79               
80       
81        #get the valid values for the Z dimension e.g. the available pressure levels
82        for ax in feature.getAxisLabels():
83            if ax not in [lat, lon, t]:
84                name=label=ax
85                domain=feature.getDomain()
86                values=domain[name].tolist()
87                axis=AxisDescription(name, label, values)
88                axisDescriptions.append(axis)   
89
90        return title, abstract, timepositions, timelimits, units, crss, axisDescriptions
91           
92    def getCoverageDescription(self):
93        pass
94       
95     
96    def map(self, **kwargs):
97        """
98        Given csml.parser.Dataset object list the names of
99        all layers available.
100       
101        @return: A mapping of layer names to ILayer implementations.
102        @raise ValueError: If no layers are available for these keywords.
103        """
104        fileoruri=kwargs['fileoruri']
105        if fileoruri in self.layermapcache.keys():
106            #we've accessed this layer map before, get it from the cache dictionary
107            return self.layermapcache[fileoruri]
108       
109        if self.connector.isGroup(fileoruri):
110            self.datasetName = fileoruri
111        else:
112           
113            ds = self.connector.getCsmlDoc(fileoruri)
114            try:
115                self.datasetName=ds.name.CONTENT
116            except AttributeError:
117                self.datasetName = 'CSML WCS Service'
118           
119           
120        self._crscat=csml.csmllibs.csmlcrs.CRSCatalogue()
121       
122        if self.connector.isGroup(fileoruri):
123            coverages = self._getCoveragesFromFolder(folderPath=fileoruri)
124        else:
125            coverages = self._getCoveragesFromCSMLFile(fileoruri)
126       
127        layermap = {}
128        for c in coverages:
129            layermap[c.name] = c
130       
131#        for feature in csml.csmllibs.csmlextra.listify(ds.featureCollection.featureMembers):
132#            title, abstract, timepositions, timelimits, units, crss, axisDescriptions=self.getInfo(feature)
133#           
134#            name = feature.id
135#            layermap[feature.id]=CSMLCoverage(name, [title],[abstract],
136#                                               timepositions, timelimits, units,
137#                                               crss, axisDescriptions, feature)
138#           
139        if len(layermap) > 0:
140            self.layermapcache[fileoruri]=layermap
141            return layermap
142        else:
143            raise ValueError
144   
145    def _getCoveragesFromFolder(self, folderPath=None):
146        l = []
147       
148        for fc in self.connector.list(folder=folderPath):
149           
150            if self.connector.isGroup(fc):
151                l += self._getCoveragesFromFolder(folderPath=fc)
152            else: 
153                l += self._getCoveragesFromCSMLFile(fc)
154               
155        return l
156   
157    def _getCoveragesFromCSMLFile(self, fileoruri):
158        l = []
159        ds = self.connector.getCsmlDoc(fileoruri)
160       
161        for feature in csml.csmllibs.csmlextra.listify(ds.featureCollection.featureMembers):
162            title, abstract, timepositions, timelimits, units, crss, axisDescriptions=self.getInfo(feature)
163           
164            name = fileoruri.replace('/','_') + '_' +feature.id
165            c = CSMLCoverage(name, [title],[abstract], timepositions, timelimits, 
166                             units,crss, axisDescriptions, feature)
167           
168            l.append(c)
169   
170        return l
171                         
172       
173class AxisDescription(object):
174    """ represents an axisDescription from the rangeSet (see wcs 1.0.0 describe coverage) """
175    def __init__(self, name, label, values):
176        self.name=name
177        self.label=label
178        self.values=values
179
180class CSMLCoverage(object): #TODO: define ICoverage
181    """ represents a WCS Coverage. Implements ICoverage """
182   
183    def __init__(self, name, title, abstract, timepositions, timelimits, units, crss, axisDescriptions, feature):
184        self.name = name
185        self.title=title
186        self.abstract=abstract
187        self.description='TO DO - coverage description'
188        self.timePositions=timepositions
189        self.timeLimits=timelimits
190        self.units=units
191        self.crss=crss
192        self._feature=feature
193        self.id=feature.id
194        self.bboxes=[]
195               
196        csmlbb = self._feature.getCSMLBoundingBox()
197        twoDCRS=csmlbb.get2DCRSName()
198        #default to global bounding box:
199        bb=[-180,90,180,90]
200        if twoDCRS=='EPSG:4326':
201            try:
202                bb=csmlbb.getBox()
203                #convert 0 - 360 to -180, 180 as per common WXS convention
204                if abs(bb[2]-bb[0]) >= 359 and abs(bb[2]-bb[0]) < 361:
205                    bb[0], bb[2]=-180, 180       
206            except:
207                #failed to get better bbox info from CSML file:
208                pass
209        self.wgs84BBox = WGS84BoundingBox(bb[:2],
210                                         bb[2:])
211        #now create list of all other bounding boxes:
212        if twoDCRS not in ['EPSG:4326', None]:
213            #then add this bbox to the list:
214            alternativebbox=csmlbb.getBox()
215            self.bboxes.append(BoundingBox([alternativebbox[0],alternativebbox[1]], [alternativebbox[2],alternativebbox[3]], twoDCRS))
216        self.featureInfoFormats = ['text/html']
217        self.axisDescriptions = axisDescriptions
218   
219    def getBBox(self, crs):
220        """
221        @return: A 4-typle of the bounding box in the given coordinate
222            reference system.
223        """
224        #TODO: make this generic
225        return self.wgs84BBox
226   
227    def getCvg(self, bbox, time=None, crs=None, response_crs=None, **kwargs):
228        """
229        Creates a subset of the layer in a particular CRS and set of
230        dimensions.
231        #TODO: this is all out of synch
232        @param crs: The coordinate reference system.
233        @param dimValues: A mapping of dimension names to dimension values
234            as specified in the IDimension.extent
235        @param renderOpts: A generic mapping object for passing rendering
236            options
237        @return: An object implementing ILayerSlab
238        #create netcdf for whole lat/lon for given dimValues, use to init slab
239        """
240
241        log.debug('WCS: getSubset(%s, %s, %s)' % (bbox, time, crs))
242        #unpack the Boundingbox.
243       
244       
245        ############################# from old WCS stack #################
246        boundingbox=bbox
247        lon=self._feature.getLongitudeAxis()
248        lat=self._feature.getLatitudeAxis()
249        t=self._feature.getTimeAxis()
250        if None in [lon, lat, t]:
251            #TODO need to return a suitable wcs error.
252            log.debug('warning, could not get correct axis info')
253            #best guess!
254            if t is None:
255                t='time'
256            if lon is None:
257                lon = 'longitude'
258            if lat is None:
259                lat = 'latitude'
260       
261        #create selection dictionary:
262        sel={}
263        sel[lat]=(boundingbox[1], boundingbox[3])
264        sel[lon]=(boundingbox[0], boundingbox[2])
265        if time is not None:
266            if  type(time) is unicode:
267                sel[t]=str(time)
268            else:
269                sel[t]=time
270               
271        #z is the 4th axis/band (eg height or pressure or wavelength) requested via kwargs as defined in the rangeset.axisDescriptions.
272        for kw in kwargs:
273            log.debug(kw)
274            log.debug(self._feature.getAxisLabels())
275            for ax in self._feature.getAxisLabels():
276                if ax not in [lat, lon, t]:
277                    if ax == kw:
278                        z=str(ax)
279                        sel[z]=(kwargs[kw])
280                        log.debug('Z axis: %s'%z)     
281        ##################################################################
282        log.debug('Final selection being made to the csml api %s'%str(sel))
283        filename = extractToNetCDF(self._feature, sel)
284        return filename
285       
Note: See TracBrowser for help on using the repository browser.