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

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

First attempt at trying to implement a folder structure for the csml files.

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           
112        layermap={}
113        self._crscat=csml.csmllibs.csmlcrs.CRSCatalogue()
114       
115        for feature in csml.csmllibs.csmlextra.listify(ds.featureCollection.featureMembers):
116            title, abstract, timepositions, timelimits, units, crss, axisDescriptions=self.getInfo(feature)
117            layermap[feature.id]=CSMLCoverage([title],[abstract], timepositions, timelimits, units, crss, axisDescriptions, feature)
118           
119        if len(layermap) > 0:
120            self.layermapcache[fileoruri]=layermap
121            return layermap
122        else:
123            raise ValueError
124
125class AxisDescription(object):
126    """ represents an axisDescription from the rangeSet (see wcs 1.0.0 describe coverage) """
127    def __init__(self, name, label, values):
128        self.name=name
129        self.label=label
130        self.values=values
131
132class CSMLCoverage(object): #TODO: define ICoverage
133    """ represents a WCS Coverage. Implements ICoverage """
134   
135    def __init__(self, title, abstract, timepositions, timelimits, units, crss, axisDescriptions, feature):
136        self.title=title
137        self.abstract=abstract
138        self.description='TO DO - coverage description'
139        self.timePositions=timepositions
140        self.timeLimits=timelimits
141        self.units=units
142        self.crss=crss
143        self._feature=feature
144        self.id=feature.id
145        self.bboxes=[]
146               
147        csmlbb = self._feature.getCSMLBoundingBox()
148        twoDCRS=csmlbb.get2DCRSName()
149        #default to global bounding box:
150        bb=[-180,90,180,90]
151        if twoDCRS=='EPSG:4326':
152            try:
153                bb=csmlbb.getBox()
154                #convert 0 - 360 to -180, 180 as per common WXS convention
155                if abs(bb[2]-bb[0]) >= 359 and abs(bb[2]-bb[0]) < 361:
156                    bb[0], bb[2]=-180, 180       
157            except:
158                #failed to get better bbox info from CSML file:
159                pass
160        self.wgs84BBox = WGS84BoundingBox(bb[:2],
161                                         bb[2:])
162        #now create list of all other bounding boxes:
163        if twoDCRS not in ['EPSG:4326', None]:
164            #then add this bbox to the list:
165            alternativebbox=csmlbb.getBox()
166            self.bboxes.append(BoundingBox([alternativebbox[0],alternativebbox[1]], [alternativebbox[2],alternativebbox[3]], twoDCRS))
167        self.featureInfoFormats = ['text/html']
168        self.axisDescriptions = axisDescriptions
169   
170    def getBBox(self, crs):
171        """
172        @return: A 4-typle of the bounding box in the given coordinate
173            reference system.
174        """
175        #TODO: make this generic
176        return self.wgs84BBox
177   
178    def getCvg(self, bbox, time=None, crs=None, response_crs=None, **kwargs):
179        """
180        Creates a subset of the layer in a particular CRS and set of
181        dimensions.
182        #TODO: this is all out of synch
183        @param crs: The coordinate reference system.
184        @param dimValues: A mapping of dimension names to dimension values
185            as specified in the IDimension.extent
186        @param renderOpts: A generic mapping object for passing rendering
187            options
188        @return: An object implementing ILayerSlab
189        #create netcdf for whole lat/lon for given dimValues, use to init slab
190        """
191
192        log.debug('WCS: getSubset(%s, %s, %s)' % (bbox, time, crs))
193        #unpack the Boundingbox.
194       
195       
196        ############################# from old WCS stack #################
197        boundingbox=bbox
198        lon=self._feature.getLongitudeAxis()
199        lat=self._feature.getLatitudeAxis()
200        t=self._feature.getTimeAxis()
201        if None in [lon, lat, t]:
202            #TODO need to return a suitable wcs error.
203            log.debug('warning, could not get correct axis info')
204            #best guess!
205            if t is None:
206                t='time'
207            if lon is None:
208                lon = 'longitude'
209            if lat is None:
210                lat = 'latitude'
211       
212        #create selection dictionary:
213        sel={}
214        sel[lat]=(boundingbox[1], boundingbox[3])
215        sel[lon]=(boundingbox[0], boundingbox[2])
216        if time is not None:
217            if  type(time) is unicode:
218                sel[t]=str(time)
219            else:
220                sel[t]=time
221               
222        #z is the 4th axis/band (eg height or pressure or wavelength) requested via kwargs as defined in the rangeset.axisDescriptions.
223        for kw in kwargs:
224            log.debug(kw)
225            log.debug(self._feature.getAxisLabels())
226            for ax in self._feature.getAxisLabels():
227                if ax not in [lat, lon, t]:
228                    if ax == kw:
229                        z=str(ax)
230                        sel[z]=(kwargs[kw])
231                        log.debug('Z axis: %s'%z)     
232        ##################################################################
233        log.debug('Final selection being made to the csml api %s'%str(sel))
234        filename = extractToNetCDF(self._feature, sel)
235        return filename
236       
Note: See TracBrowser for help on using the repository browser.