source: cows/branches/wcsmerge/cows/service/imps/csmlbackend/wcs_csmllayer.py @ 4633

Subversion URL: http://proj.badc.rl.ac.uk/svn/ndg/cows/branches/wcsmerge/cows/service/imps/csmlbackend/wcs_csmllayer.py@4633
Revision 4633, 17.3 KB checked in by domlowe, 12 years ago (diff)

More progress on merging.

Line 
1"""
2implementation of ILayerMapper, IwmsLayer, IwmsDimension, ILayerSlab interfaces, as defined in wms_iface.py & wxs_iface.py
3
4"""
5import os, string
6import csml
7try:
8    import cdms2 as cdms
9except:
10    import cdms
11from pylons import config, request, session
12import Image
13from copy import copy
14from matplotlib import cm
15import cdtime
16import logging
17log = logging.getLogger(__name__)
18
19from cows.model.contents import DatasetSummary
20from cows.service.wms_iface import IwmsLayer, IwmsDimension, IwmsLayerSlab
21from cows.service.imps.csmlbackend.csmlcommon import CSMLLayerMapper, CSMLConnector, extractToNetCDF
22from cows.service.imps.csmlbackend.wms_csmllayer import CSMLwmsLayerMapper, CSMLwmsLayer
23
24
25class CSMLwcsCoverageMapper(CSMLLayerMapper):
26    """
27    Map keyword arguments to a collection of layers (coverages)
28    Supports the retrieval of coverages according to arbitrary
29    keyword/value pairs. 
30    """
31    def __init__(self):
32        super(CSMLwcsCoverageMapper, self).__init__()
33       
34   
35    def getInfo(self, feature):
36        ''' given a csml feature, return info about the layer/feature
37        @return:   title, abstract, dimensions, units, crss '''
38
39        title, abstract = super(CSMLwcsCoverageMapper, self).getInfo(feature)
40        units=feature.getDomainUnits()
41        dimensions={}
42        tmpunits=copy(units)
43        tmpunits.reverse()
44        domain = feature.getDomain()
45        for dim in feature.getAxisLabels():
46            nextdim=CSMLwmsDimension(domain, dim, tmpunits.pop())
47            if dim not in ['latitude', 'longitude']:
48                dimensions[dim]=nextdim
49        crs=feature.getNativeCRS()
50        crss=[self._crscat.getCRS(crs).twoD]
51        if 'EPSG:4326' in crss:
52            crss.append('CRS:84')
53            crss.append('WGS84')
54        return title, abstract, dimensions, units, crss
55       
56       
57     
58    def map(self, **kwargs):
59        """
60        Given csml.parser.Dataset object list the names of
61        all layers available.
62       
63        @return: A mapping of layer names to ILayer implementations.
64        @raise ValueError: If no layers are available for these keywords.
65        """
66        fileoruri=kwargs['fileoruri']
67        if fileoruri in self.layermapcache.keys():
68            #we've accessed this layer map before, get it from the cache dictionary
69            return self.layermapcache[fileoruri]
70         
71        ds = self.connector.getCsmlDoc(fileoruri)
72        layermap={}
73        self._crscat=csml.csmllibs.csmlcrs.CRSCatalogue()
74        for feature in csml.csmllibs.csmlextra.listify(ds.featureCollection.featureMembers):
75            title, abstract, dimensions, units, crss=self.getInfo(feature)
76            layermap[feature.id]=CSMLCoverage(title,abstract, dimensions, units, crss, feature)
77        if len(layermap) > 0:
78            self.layermapcache[fileoruri]=layermap
79            return layermap
80        else:
81            raise ValueError
82
83class CoverageDescription(DatasetSummary):
84    """
85    Extends DatasetSummary from cows.model.contents
86    """
87    def __init__(self, identifier, timeDomain, crs=None, **kw):
88        super(CoverageDescription, self).__init__(**kw)
89        self.identifier=identifier
90        self.timeLimits=timeDomain
91        self.crs=crs
92
93
94class CSMLCoverage(): #TODO implements ICoverage
95    """ represents a WCS Coverage. Implements ICoverage """
96   
97    def __init__(self, title, abstract, dimensions, units, crss, feature):
98        self.title=title
99        self.abstract=abstract
100        self.dimensions=dimensions
101        self.description='TO DO - coverage description'
102        self.timeLimits='t1, t2'
103        self.units=units
104        self.crss=crss
105        self._feature=feature
106        self.id=feature.id
107        bb= self._feature.getCSMLBoundingBox().getBox()
108        #convert 0 - 360 to -180, 180 as per common WXS convention
109        if abs(bb[2]-bb[0]) >= 359 and abs(bb[2]-bb[0]) < 361:
110            bb[0], bb[2]=-180, 180
111        self.wgs84BBox = bb
112        self.featureInfoFormats = ['text/html']
113       
114    def getDescription(self):
115        pass
116   
117    def getBBox(self, crs):
118        """
119        @return: A 4-typle of the bounding box in the given coordinate
120            reference system.
121        """
122        #TODO: make this generic
123        return self.wgs84BBox
124   
125    def getCvg(self, bbox, time=None, crs=None, response_crs=None):
126        """
127        Creates a subset of the layer in a particular CRS and set of
128        dimensions.
129        #TODO: this is all out of synch
130        @param crs: The coordinate reference system.
131        @param dimValues: A mapping of dimension names to dimension values
132            as specified in the IDimension.extent
133        @param renderOpts: A generic mapping object for passing rendering
134            options
135        @return: An object implementing ILayerSlab
136        #create netcdf for whole lat/lon for given dimValues, use to init slab
137        """
138
139        log.debug('WCS: getSubset(%s, %s, %s)' % (bbox, time, crs))
140        #unpack the Boundingbox.
141       
142       
143        ############################# from old WCS stack #################
144        boundingbox=bbox
145        lon=self._feature.getLongitudeAxis()
146        lat=self._feature.getLatitudeAxis()
147        t=self._feature.getTimeAxis()
148        if None in [lon, lat, t]:
149            #TODO need to return a suitable wcs error.
150            print 'warning, could not get correct axis info'
151            #best guess!
152            if t is None:
153                t='time'
154            if lon is None:
155                lon = 'longitude'
156            if lat is None:
157                lat = 'latitude'
158       
159        #create selection dictionary:
160        sel={}
161        sel[lat]=(boundingbox[1], boundingbox[3])
162        sel[lon]=(boundingbox[0], boundingbox[2])
163        if time is not None:
164            if  type(time) is unicode:
165                sel[t]=str(time)
166            else:
167                sel[t]=time
168        #z is the 4th axis (eg height or pressure).
169        #NOTE, need to decide whether bounding box is of form: x1,y1,z1,x2,y2,z2 or x1,y1,x2,y2,z1,z2
170        #currently the latter is implemented.
171           
172        if len(boundingbox)  == 6:
173            for ax in self._feature.getAxisLabels():
174                if ax not in [lat, lon, t]:
175                    #must be Z 
176                    z=str(ax)
177                    sel[z]=(boundingbox[4], boundingbox[5])
178           
179           
180        axisNames=self._feature.getAxisLabels()
181        ##################################################################
182       
183        filename = extractToNetCDF(self._feature, sel)
184        return filename
185       
186       
187       
188#       
189#        #unicode conversion
190#        for dimval in dimValues:
191#            if dimval != 'time':
192#                dimValues[dimval]=float(dimValues[dimval])
193#            else:
194#                #remove any trailing Zs from time string
195#                if dimValues[dimval] [-1:] in ['Z', 'z']:
196#                    dimValues[dimval]=dimValues[dimval][:-1]
197#        if type(self._feature) == csml.parser.GridSeriesFeature:
198#            randomname= csml.csmllibs.csmlextra.getRandomID() + '.nc'
199#            result= self._feature.subsetToGridSeries(config['tmpfilebuffer'], ncname=randomname, **dimValues)
200#            return result[1]
201#        else:
202#            raise NotImplementedError
203   
204   
205    def getCacheKey(self, crs, dimValues=None):
206        """
207        Create a unique key for use in caching a slab.
208
209        The intention here is that most of the work should be done when
210        instantiating an ILayerSlab object.  These can be cached by the
211        server for future use.  The server will first call getCacheKey()
212        for the slab creation arguments and if the key is in it's cache
213        it will use a pre-generated ILayerSlab object.
214
215        """
216        dimList = list(dimValues.items())
217        dimList.sort()
218        return '%s:%s' % (crs, dimList)
219   
220#class CSMLwcsLayer(CSMLwmsLayer):
221#    """
222#     representing a WMS layer.    Implements IwmsLayer
223#
224#    @ivar title: The layer title.  As seen in the Capabilities document.
225#    @ivar abstract:  Abstract as seen in the Capabilities document.
226#    @ivar dimensions: A dictionary of IDimension objects.
227#    @ivar units: A string describing the units.
228#    @ivar crss: A sequence of SRS/CRSs supported by this layer.
229#
230#    @todo: Do we need minValue/maxValue?
231#
232#    """
233#   
234#    def __init__(self, title, abstract, dimensions, units, crss, feature):
235#        self.featureInfoFormats=None #NotImplemented
236#        self.title=title
237#        self.abstract=abstract
238#        self.dimensions=dimensions
239#        self.units=units
240#        self.crss=crss
241#        self._feature=feature
242#        self.legendSize=(30,100)
243#        bb= self._feature.getCSMLBoundingBox().getBox()
244#        #convert 0 - 360 to -180, 180 as per common WMS convention
245#        if abs(bb[2]-bb[0]) >= 359 and abs(bb[2]-bb[0]) < 361:
246#            bb[0], bb[2]=-180, 180
247#        self.wgs84BBox = bb
248#        self.featureInfoFormats = ['text/html']
249#        try:
250#            self.wgs84BBox = self.getBBox('EPSG:4326')
251#        except:
252#            raise ValueError("Layer must provide a bounding box in EPSG:4326 "
253#                             "coordinates for compatibility with WMS-1.3.0")
254#        self.featureinfofilecache={} #used for caching netcdf file in getFeatureInfo
255#       
256#    def getBBox(self, crs):
257#        """
258#        @return: A 4-typle of the bounding box in the given coordinate
259#            reference system.
260#        """
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#        return self.wgs84BBox
267#        #raise NotImplementedError
268#   
269#    def getSubset(self, crs, dimValues=None):
270#        """
271#        Creates a subset of the layer in a particular CRS and set of
272#        dimensions.
273#
274#        @param crs: The coordinate reference system.
275#        @param dimValues: A mapping of dimension names to dimension values
276#            as specified in the IDimension.extent
277#        @param renderOpts: A generic mapping object for passing rendering
278#            options
279#        @return: An object implementing ILayerSlab
280#        #create netcdf for whole lat/lon for given dimValues, use to init slab
281#        """
282#
283#        log.debug('getSlab(%s, %s)' % (crs, dimValues))
284#       
285#        #unicode conversion
286#        for dimval in dimValues:
287#            if dimval != 'time':
288#                dimValues[dimval]=float(dimValues[dimval])
289#            else:
290#                #remove any trailing Zs from time string
291#                if dimValues[dimval] [-1:] in ['Z', 'z']:
292#                    dimValues[dimval]=dimValues[dimval][:-1]
293#        if type(self._feature) == csml.parser.GridSeriesFeature:
294#            randomname= csml.csmllibs.csmlextra.getRandomID() + '.nc'
295#            result= self._feature.subsetToGridSeries(config['tmpfilebuffer'], ncname=randomname, **dimValues)
296#            return result
297#        else:
298#            raise NotImplementedError
299#       
300##    def getSlab(self, crs, dimValues=None, renderOpts={}):
301##        """
302##        Creates a slab of the layer in a particular CRS and set of
303##        dimensions.
304##
305##        @param crs: The coordinate reference system.
306##        @param dimValues: A mapping of dimension names to dimension values
307##            as specified in the IDimension.extent
308##        @param renderOpts: A generic mapping object for passing rendering
309##            options
310##        @return: An object implementing ILayerSlab
311##        #create netcdf for whole lat/lon for given dimValues, use to init slab
312##        """
313##
314##        log.debug('getSlab(%s, %s)' % (crs, dimValues))
315##       
316##        #unicode conversion
317##        for dimval in dimValues:
318##            if dimval != 'time':
319##                dimValues[dimval]=float(dimValues[dimval])
320##            else:
321##                #remove any trailing Zs from time string
322##                if dimValues[dimval] [-1:] in ['Z', 'z']:
323##                    dimValues[dimval]=dimValues[dimval][:-1]
324##        if type(self._feature) == csml.parser.GridSeriesFeature:
325##            randomname= csml.csmllibs.csmlextra.getRandomID() + '.nc'
326##            result= self._feature.subsetToGridSeries(config['tmpfilebuffer'], ncname=randomname, **dimValues)
327##            #for now have to read netcdf back from
328##            #disk (limitiation of CSML api)
329##            netcdf=cdms.open(result[1])
330##            #and then delete the temporary file
331##            os.system('rm %s'%result[1])
332##            bbox=self.getBBox(crs)
333##            return CSMLwmsLayerSlab(netcdf, self, crs, dimValues, renderOpts, bbox)
334##        else:
335##            raise NotImplementedError
336#       
337#    def getCacheKey(self, crs, dimValues=None, renderOpts={}):
338#        """
339#        Create a unique key for use in caching a slab.
340#
341#        The intention here is that most of the work should be done when
342#        instantiating an ILayerSlab object.  These can be cached by the
343#        server for future use.  The server will first call getCacheKey()
344#        for the slab creation arguments and if the key is in it's cache
345#        it will use a pre-generated ILayerSlab object.
346#
347#        """
348#        dimList = list(dimValues.items())
349#        dimList.sort()
350#        return '%s:%s' % (crs, dimList)
351#
352#    def getFeatureInfo(self, format, crs, point, dimValues):
353#        """
354#        Return a response string descibing the feature at a given
355#        point in a given CRS.
356#
357#        Currently only "html" is supported as output format
358#
359#        @param format: One of self.featureInfoFormats.  Defines which
360#            format the response will be in.
361#        @param crs: One of self.crss
362#        @param point: a tuple (x, y) in the supplied crs of the point
363#            being selected.
364#        @param dimValues: A mapping of dimension names to dimension values.
365#        @return: A string containing the response.
366#
367#        """
368#       
369#        #cached netcdf is indexed by a tuple of feature id and dimvalues - i.e. typically a particular time and Z value for that feature.
370#        #look in dictionary for cached copy, and if so use that as the file object.
371#        dictindex=str((self._feature.id, dimValues))
372#        if dictindex in self.featureinfofilecache:
373#            log.debug('calling cache')
374#            f=self.featureinfofilecache[dictindex]
375#        else: #else, use the csml api to subset the feature afresh
376#            log.debug('not calling cache')
377#            randomname= csml.csmllibs.csmlextra.getRandomID() + '.nc'
378#            result= self._feature.subsetToGridSeries(config['tmpfilebuffer'], ncname=randomname, **dimValues)
379#            #for now have to read netcdf back from disk (limitation of CSML api)
380#            f=cdms.open(result[1])
381#            #append to cache:
382#            self.featureinfofilecache[dictindex]=f
383#            #and then delete the temporary file
384#            os.system('rm %s'%result[1])
385#       
386#        netcdf = f(self.title)  #netcdf here is a cdms transient variable
387#       
388#       
389#        #Now grab the netCDF object for the point specified.
390#        #The reason for the 'cob' option is so that if the grid the data
391#        #is defined on does not have a grid point at the point specified,
392#        #we should  still get the nearest location
393#       
394#        t_point = netcdf(latitude=(point[1], point[1], 'cob'), longitude=(point[0], point[0], 'cob'))
395#        #now get the value recorded at this location
396#        value = t_point.getValue().tolist()
397#        log.debug(value)
398#        log.debug(t_point.fill_value())
399#        #and the fill_value too
400#        fill_value = t_point.fill_value()
401#        #value is actually embedded in a multi dimensional list,
402#        #so we need to extract the actual value from the list
403#        while type(value) is list:
404#                value = value[0]
405#
406#        #now check if the value is actually the fill_value rather than
407#        #a value recorded at the point specified
408#        log.debug('%s %s' % (value, fill_value))
409#        if (2*fill_value) == value:
410#                value = "No value found at position: "+str(point[1])+", "+str(point[0])
411#        else:
412#                value = "Value found at position: "+str(point[1])+", "+str(point[0])+" is: "+str(value)
413#        # finally return the value
414#        return value
415
416
417class CSMLwmsDimension(IwmsDimension):
418    """
419    implements IDimension
420    @ivar units: The units string.
421    @ivar extent: Sequence of extent values.
422
423    """
424   
425    def __init__(self, domain, dimname, unit):
426        self.units = unit
427        self.extent = []
428        #patch to handle current limitations of multiple time dimension scanning in csml.
429        if string.lower(self.units)[:10] in ['days_since', 'seconds_si', 'minutes_si', 'hours_sinc','months _sin', 'years_sinc']:
430            if type(domain[dimname][0]) is not str   :
431                tunits=self.units.replace('_', ' ')
432                for val in domain[dimname]:
433                    csmltime= csml.csmllibs.csmltime.UDtimeToCSMLtime(cdtime.reltime(float(val), tunits).tocomp())
434                    self.extent.append(csmltime)
435            else:
436                for val in domain[dimname]:
437                    self.extent.append(str(val))
438        else:
439            for val in domain[dimname]:
440                self.extent.append(str(val))
441        #for time axis replace units with iso string
442        if dimname == 'time':
443            self.units='ISO8601'
444       
Note: See TracBrowser for help on using the repository browser.