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

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

Modifying wcs and wfs layermappers to provide dataset name

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        try:
49            timepositions=domain[tax]
50            timelimits=[domain[tax][0]+'Z',domain[tax][len(domain[tax])-1]+'Z']
51        except: #no time available, return empty strings
52            timepositions=['']
53            timelimits=['',''] 
54        crs=feature.getNativeCRS()
55        log.debug('crs %s'%crs)
56        crss=[self._crscat.getCRS(crs).twoD] 
57        if 'EPSG:4326' in crss:
58            crss.append('CRS:84')
59            crss.append('WGS84')
60        log.debug('crss %s'%crss)
61        #build definitions of any Z axes such as air_pressure:
62        axisDescriptions=[]       
63        lon=feature.getLongitudeAxis()
64        lat=feature.getLatitudeAxis()
65        t=feature.getTimeAxis()
66        if None in [lon, lat, t]:
67            #TODO need to return a suitable wcs error.
68            log.debug('warning, could not get correct axis info')
69            #best guess!
70            if t is None:
71                t='time'
72            if lon is None:
73                lon = 'longitude'
74            if lat is None:
75                lat = 'latitude'
76               
77       
78        #get the valid values for the Z dimension e.g. the available pressure levels
79        for ax in feature.getAxisLabels():
80            if ax not in [lat, lon, t]:
81                name=label=ax
82                domain=feature.getDomain()
83                values=domain[name].tolist()
84                axis=AxisDescription(name, label, values)
85                axisDescriptions.append(axis)   
86
87        return title, abstract, timepositions, timelimits, units, crss, axisDescriptions
88           
89    def getCoverageDescription(self):
90        pass
91       
92     
93    def map(self, **kwargs):
94        """
95        Given csml.parser.Dataset object list the names of
96        all layers available.
97       
98        @return: A mapping of layer names to ILayer implementations.
99        @raise ValueError: If no layers are available for these keywords.
100        """
101        fileoruri=kwargs['fileoruri']
102        if fileoruri in self.layermapcache.keys():
103            #we've accessed this layer map before, get it from the cache dictionary
104            return self.layermapcache[fileoruri]
105         
106        ds = self.connector.getCsmlDoc(fileoruri)
107        try:
108            self.datasetName=ds.name.CONTENT
109        except AttributeError:
110            self.datasetName = 'CSML WCS Service'
111        layermap={}
112        self._crscat=csml.csmllibs.csmlcrs.CRSCatalogue()
113        for feature in csml.csmllibs.csmlextra.listify(ds.featureCollection.featureMembers):
114            title, abstract, timepositions, timelimits, units, crss, axisDescriptions=self.getInfo(feature)
115            layermap[feature.id]=CSMLCoverage([title],[abstract], timepositions, timelimits, units, crss, axisDescriptions, feature)
116        if len(layermap) > 0:
117            self.layermapcache[fileoruri]=layermap
118            return layermap
119        else:
120            raise ValueError
121
122class AxisDescription(object):
123    """ represents an axisDescription from the rangeSet (see wcs 1.0.0 describe coverage) """
124    def __init__(self, name, label, values):
125        self.name=name
126        self.label=label
127        self.values=values
128
129class CSMLCoverage(object): #TODO: define ICoverage
130    """ represents a WCS Coverage. Implements ICoverage """
131   
132    def __init__(self, title, abstract, timepositions, timelimits, units, crss, axisDescriptions, feature):
133        self.title=title
134        self.abstract=abstract
135        self.description='TO DO - coverage description'
136        self.timePositions=timepositions
137        self.timeLimits=timelimits
138        self.units=units
139        self.crss=crss
140        self._feature=feature
141        self.id=feature.id
142        self.bboxes=[]
143               
144        csmlbb = self._feature.getCSMLBoundingBox()
145        twoDCRS=csmlbb.get2DCRSName()
146        #default to global bounding box:
147        bb=[-180,90,180,90]
148        if twoDCRS=='EPSG:4326':
149            try:
150                bb=csmlbb.getBox()
151                #convert 0 - 360 to -180, 180 as per common WXS convention
152                if abs(bb[2]-bb[0]) >= 359 and abs(bb[2]-bb[0]) < 361:
153                    bb[0], bb[2]=-180, 180       
154            except:
155                #failed to get better bbox info from CSML file:
156                pass
157        self.wgs84BBox = WGS84BoundingBox(bb[:2],
158                                         bb[2:])
159        #now create list of all other bounding boxes:
160        if twoDCRS not in ['EPSG:4326', None]:
161            #then add this bbox to the list:
162            alternativebbox=csmlbb.getBox()
163            self.bboxes.append(BoundingBox([alternativebbox[0],alternativebbox[1]], [alternativebbox[2],alternativebbox[3]], twoDCRS))
164        self.featureInfoFormats = ['text/html']
165        self.axisDescriptions = axisDescriptions
166   
167    def getBBox(self, crs):
168        """
169        @return: A 4-typle of the bounding box in the given coordinate
170            reference system.
171        """
172        #TODO: make this generic
173        return self.wgs84BBox
174   
175    def getCvg(self, bbox, time=None, crs=None, response_crs=None, **kwargs):
176        """
177        Creates a subset of the layer in a particular CRS and set of
178        dimensions.
179        #TODO: this is all out of synch
180        @param crs: The coordinate reference system.
181        @param dimValues: A mapping of dimension names to dimension values
182            as specified in the IDimension.extent
183        @param renderOpts: A generic mapping object for passing rendering
184            options
185        @return: An object implementing ILayerSlab
186        #create netcdf for whole lat/lon for given dimValues, use to init slab
187        """
188
189        log.debug('WCS: getSubset(%s, %s, %s)' % (bbox, time, crs))
190        #unpack the Boundingbox.
191       
192       
193        ############################# from old WCS stack #################
194        boundingbox=bbox
195        lon=self._feature.getLongitudeAxis()
196        lat=self._feature.getLatitudeAxis()
197        t=self._feature.getTimeAxis()
198        if None in [lon, lat, t]:
199            #TODO need to return a suitable wcs error.
200            log.debug('warning, could not get correct axis info')
201            #best guess!
202            if t is None:
203                t='time'
204            if lon is None:
205                lon = 'longitude'
206            if lat is None:
207                lat = 'latitude'
208       
209        #create selection dictionary:
210        sel={}
211        sel[lat]=(boundingbox[1], boundingbox[3])
212        sel[lon]=(boundingbox[0], boundingbox[2])
213        if time is not None:
214            if  type(time) is unicode:
215                sel[t]=str(time)
216            else:
217                sel[t]=time
218               
219        #z is the 4th axis/band (eg height or pressure or wavelength) requested via kwargs as defined in the rangeset.axisDescriptions.
220        for kw in kwargs:
221            log.debug(kw)
222            log.debug(self._feature.getAxisLabels())
223            for ax in self._feature.getAxisLabels():
224                if ax not in [lat, lon, t]:
225                    if ax == kw:
226                        z=str(ax)
227                        sel[z]=(kwargs[kw])
228                        log.debug('Z axis: %s'%z)     
229        ##################################################################
230        log.debug('Final selection being made to the csml api %s'%str(sel))
231        filename = extractToNetCDF(self._feature, sel)
232        return filename
233       
Note: See TracBrowser for help on using the repository browser.