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

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

more 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
22from cows.service.imps.csmlbackend.wms_csmllayer import CSMLwmsLayerMapper, CSMLwmsLayer
23
24class CSMLwcsLayerMapper(CSMLLayerMapper):
25    """
26    Map keyword arguments to a collection of layers.
27    Supports the retrieval of sets of layers according to arbitrary
28    keyword/value pairs.
29    Implements  CSMLLayerMapper
30   
31    """
32    def __init__(self):
33        super(CSMLwcsLayerMapper, self).__init__()
34   
35   
36   
37    def getInfo(self, feature):
38        ''' given a csml feature, return info about the layer/feature
39        @return:   title, abstract, dimensions, units, crss '''
40
41        title, abstract = super(CSMLwcsLayerMapper, self).getInfo(feature)
42        units=feature.getDomainUnits()
43        dimensions={}
44        tmpunits=copy(units)
45        tmpunits.reverse()
46        domain = feature.getDomain()
47        for dim in feature.getAxisLabels():
48            nextdim=CSMLwmsDimension(domain, dim, tmpunits.pop())
49            if dim not in ['latitude', 'longitude']:
50                dimensions[dim]=nextdim
51        crs=feature.getNativeCRS()
52        crss=[self._crscat.getCRS(crs).twoD]
53        if 'EPSG:4326' in crss:
54            crss.append('CRS:84')
55            crss.append('WGS84')
56        return title, abstract, dimensions, units, crss
57       
58       
59     
60    def map(self, **kwargs):
61        """
62        Given csml.parser.Dataset object list the names of
63        all layers available.
64       
65        @return: A mapping of layer names to ILayer implementations.
66        @raise ValueError: If no layers are available for these keywords.
67        """
68        fileoruri=kwargs['fileoruri']
69        if fileoruri in self.layermapcache.keys():
70            #we've accessed this layer map before, get it from the cache dictionary
71            return self.layermapcache[fileoruri]
72         
73        ds = self.connector.getCsmlDoc(fileoruri)
74        layermap={}
75        self._crscat=csml.csmllibs.csmlcrs.CRSCatalogue()
76        for feature in csml.csmllibs.csmlextra.listify(ds.featureCollection.featureMembers):
77            title, abstract, dimensions, units, crss=self.getInfo(feature)
78            layermap[feature.id]=CSMLCoverage(title,abstract, dimensions, units, crss, feature)
79        if len(layermap) > 0:
80            self.layermapcache[fileoruri]=layermap
81            return layermap
82        else:
83            raise ValueError
84
85class CoverageDescription(DatasetSummary):
86    """
87    Extends DatasetSummary from cows.model.contents
88    """
89    def __init__(self, identifier, timeDomain, crs=None, **kw):
90        super(CoverageDescription, self).__init__(**kw)
91        self.identifier=identifier
92        self.timeLimits=timeDomain
93        self.crs=crs
94
95
96class CSMLCoverage(): #TODO implements ICoverage
97    """ represents a WCS Coverage. Implements ICoverage """
98   
99    def __init__(self, title, abstract, dimensions, units, crss, feature):
100        self.title=title
101        self.abstract=abstract
102        self.dimensions=dimensions
103        self.description='TO DO - coverage description'
104        self.timeLimits='t1, t2'
105        self.units=units
106        self.crss=crss
107        self._feature=feature
108        self.id=feature.id
109        bb= self._feature.getCSMLBoundingBox().getBox()
110        #convert 0 - 360 to -180, 180 as per common WXS convention
111        if abs(bb[2]-bb[0]) >= 359 and abs(bb[2]-bb[0]) < 361:
112            bb[0], bb[2]=-180, 180
113        self.wgs84BBox = bb
114        self.featureInfoFormats = ['text/html']
115       
116    def getDescription(self):
117        pass
118   
119    def getBBox(self, crs):
120        """
121        @return: A 4-typle of the bounding box in the given coordinate
122            reference system.
123        """
124        #TODO: make this generic
125        return self.wgs84BBox
126   
127    def getSubset(self, crs, dimValues=None):
128        """
129        Creates a subset of the layer in a particular CRS and set of
130        dimensions.
131
132        @param crs: The coordinate reference system.
133        @param dimValues: A mapping of dimension names to dimension values
134            as specified in the IDimension.extent
135        @param renderOpts: A generic mapping object for passing rendering
136            options
137        @return: An object implementing ILayerSlab
138        #create netcdf for whole lat/lon for given dimValues, use to init slab
139        """
140
141        log.debug('getSlab(%s, %s)' % (crs, dimValues))
142       
143        #unicode conversion
144        for dimval in dimValues:
145            if dimval != 'time':
146                dimValues[dimval]=float(dimValues[dimval])
147            else:
148                #remove any trailing Zs from time string
149                if dimValues[dimval] [-1:] in ['Z', 'z']:
150                    dimValues[dimval]=dimValues[dimval][:-1]
151        if type(self._feature) == csml.parser.GridSeriesFeature:
152            randomname= csml.csmllibs.csmlextra.getRandomID() + '.nc'
153            result= self._feature.subsetToGridSeries(config['tmpfilebuffer'], ncname=randomname, **dimValues)
154            return result[1]
155        else:
156            raise NotImplementedError
157   
158   
159    def getCacheKey(self, crs, dimValues=None):
160        """
161        Create a unique key for use in caching a slab.
162
163        The intention here is that most of the work should be done when
164        instantiating an ILayerSlab object.  These can be cached by the
165        server for future use.  The server will first call getCacheKey()
166        for the slab creation arguments and if the key is in it's cache
167        it will use a pre-generated ILayerSlab object.
168
169        """
170        dimList = list(dimValues.items())
171        dimList.sort()
172        return '%s:%s' % (crs, dimList)
173   
174class CSMLwcsLayer(CSMLwmsLayer):
175    """
176     representing a WMS layer.    Implements IwmsLayer
177
178    @ivar title: The layer title.  As seen in the Capabilities document.
179    @ivar abstract:  Abstract as seen in the Capabilities document.
180    @ivar dimensions: A dictionary of IDimension objects.
181    @ivar units: A string describing the units.
182    @ivar crss: A sequence of SRS/CRSs supported by this layer.
183
184    @todo: Do we need minValue/maxValue?
185
186    """
187   
188    def __init__(self, title, abstract, dimensions, units, crss, feature):
189        self.featureInfoFormats=None #NotImplemented
190        self.title=title
191        self.abstract=abstract
192        self.dimensions=dimensions
193        self.units=units
194        self.crss=crss
195        self._feature=feature
196        self.legendSize=(30,100)
197        bb= self._feature.getCSMLBoundingBox().getBox()
198        #convert 0 - 360 to -180, 180 as per common WMS convention
199        if abs(bb[2]-bb[0]) >= 359 and abs(bb[2]-bb[0]) < 361:
200            bb[0], bb[2]=-180, 180
201        self.wgs84BBox = bb
202        self.featureInfoFormats = ['text/html']
203        try:
204            self.wgs84BBox = self.getBBox('EPSG:4326')
205        except:
206            raise ValueError("Layer must provide a bounding box in EPSG:4326 "
207                             "coordinates for compatibility with WMS-1.3.0")
208        self.featureinfofilecache={} #used for caching netcdf file in getFeatureInfo
209       
210    def getBBox(self, crs):
211        """
212        @return: A 4-typle of the bounding box in the given coordinate
213            reference system.
214        """
215        #bb= self._feature.getCSMLBoundingBox().getBox()
216        #convert 0 - 360 to -180, 180 as per common WMS convention
217        #if abs(bb[2]-bb[0]) >= 359 and abs(bb[2]-bb[0]) < 361:
218        #    bb[0], bb[2]=-180, 180
219        #self.wgs84BBox = bb
220        return self.wgs84BBox
221        #raise NotImplementedError
222   
223    def getSubset(self, crs, dimValues=None):
224        """
225        Creates a subset of the layer in a particular CRS and set of
226        dimensions.
227
228        @param crs: The coordinate reference system.
229        @param dimValues: A mapping of dimension names to dimension values
230            as specified in the IDimension.extent
231        @param renderOpts: A generic mapping object for passing rendering
232            options
233        @return: An object implementing ILayerSlab
234        #create netcdf for whole lat/lon for given dimValues, use to init slab
235        """
236
237        log.debug('getSlab(%s, %s)' % (crs, dimValues))
238       
239        #unicode conversion
240        for dimval in dimValues:
241            if dimval != 'time':
242                dimValues[dimval]=float(dimValues[dimval])
243            else:
244                #remove any trailing Zs from time string
245                if dimValues[dimval] [-1:] in ['Z', 'z']:
246                    dimValues[dimval]=dimValues[dimval][:-1]
247        if type(self._feature) == csml.parser.GridSeriesFeature:
248            randomname= csml.csmllibs.csmlextra.getRandomID() + '.nc'
249            result= self._feature.subsetToGridSeries(config['tmpfilebuffer'], ncname=randomname, **dimValues)
250            return result
251        else:
252            raise NotImplementedError
253       
254#    def getSlab(self, crs, dimValues=None, renderOpts={}):
255#        """
256#        Creates a slab of the layer in a particular CRS and set of
257#        dimensions.
258#
259#        @param crs: The coordinate reference system.
260#        @param dimValues: A mapping of dimension names to dimension values
261#            as specified in the IDimension.extent
262#        @param renderOpts: A generic mapping object for passing rendering
263#            options
264#        @return: An object implementing ILayerSlab
265#        #create netcdf for whole lat/lon for given dimValues, use to init slab
266#        """
267#
268#        log.debug('getSlab(%s, %s)' % (crs, dimValues))
269#       
270#        #unicode conversion
271#        for dimval in dimValues:
272#            if dimval != 'time':
273#                dimValues[dimval]=float(dimValues[dimval])
274#            else:
275#                #remove any trailing Zs from time string
276#                if dimValues[dimval] [-1:] in ['Z', 'z']:
277#                    dimValues[dimval]=dimValues[dimval][:-1]
278#        if type(self._feature) == csml.parser.GridSeriesFeature:
279#            randomname= csml.csmllibs.csmlextra.getRandomID() + '.nc'
280#            result= self._feature.subsetToGridSeries(config['tmpfilebuffer'], ncname=randomname, **dimValues)
281#            #for now have to read netcdf back from
282#            #disk (limitiation of CSML api)
283#            netcdf=cdms.open(result[1])
284#            #and then delete the temporary file
285#            os.system('rm %s'%result[1])
286#            bbox=self.getBBox(crs)
287#            return CSMLwmsLayerSlab(netcdf, self, crs, dimValues, renderOpts, bbox)
288#        else:
289#            raise NotImplementedError
290       
291    def getCacheKey(self, crs, dimValues=None, renderOpts={}):
292        """
293        Create a unique key for use in caching a slab.
294
295        The intention here is that most of the work should be done when
296        instantiating an ILayerSlab object.  These can be cached by the
297        server for future use.  The server will first call getCacheKey()
298        for the slab creation arguments and if the key is in it's cache
299        it will use a pre-generated ILayerSlab object.
300
301        """
302        dimList = list(dimValues.items())
303        dimList.sort()
304        return '%s:%s' % (crs, dimList)
305
306    def getFeatureInfo(self, format, crs, point, dimValues):
307        """
308        Return a response string descibing the feature at a given
309        point in a given CRS.
310
311        Currently only "html" is supported as output format
312
313        @param format: One of self.featureInfoFormats.  Defines which
314            format the response will be in.
315        @param crs: One of self.crss
316        @param point: a tuple (x, y) in the supplied crs of the point
317            being selected.
318        @param dimValues: A mapping of dimension names to dimension values.
319        @return: A string containing the response.
320
321        """
322       
323        #cached netcdf is indexed by a tuple of feature id and dimvalues - i.e. typically a particular time and Z value for that feature.
324        #look in dictionary for cached copy, and if so use that as the file object.
325        dictindex=str((self._feature.id, dimValues))
326        if dictindex in self.featureinfofilecache:
327            log.debug('calling cache')
328            f=self.featureinfofilecache[dictindex]
329        else: #else, use the csml api to subset the feature afresh
330            log.debug('not calling cache')
331            randomname= csml.csmllibs.csmlextra.getRandomID() + '.nc'
332            result= self._feature.subsetToGridSeries(config['tmpfilebuffer'], ncname=randomname, **dimValues)
333            #for now have to read netcdf back from disk (limitation of CSML api)
334            f=cdms.open(result[1])
335            #append to cache:
336            self.featureinfofilecache[dictindex]=f
337            #and then delete the temporary file
338            os.system('rm %s'%result[1])
339       
340        netcdf = f(self.title)  #netcdf here is a cdms transient variable
341       
342       
343        #Now grab the netCDF object for the point specified.
344        #The reason for the 'cob' option is so that if the grid the data
345        #is defined on does not have a grid point at the point specified,
346        #we should  still get the nearest location
347       
348        t_point = netcdf(latitude=(point[1], point[1], 'cob'), longitude=(point[0], point[0], 'cob'))
349        #now get the value recorded at this location
350        value = t_point.getValue().tolist()
351        log.debug(value)
352        log.debug(t_point.fill_value())
353        #and the fill_value too
354        fill_value = t_point.fill_value()
355        #value is actually embedded in a multi dimensional list,
356        #so we need to extract the actual value from the list
357        while type(value) is list:
358                value = value[0]
359
360        #now check if the value is actually the fill_value rather than
361        #a value recorded at the point specified
362        log.debug('%s %s' % (value, fill_value))
363        if (2*fill_value) == value:
364                value = "No value found at position: "+str(point[1])+", "+str(point[0])
365        else:
366                value = "Value found at position: "+str(point[1])+", "+str(point[0])+" is: "+str(value)
367        # finally return the value
368        return value
369
370
371class CSMLwmsDimension(IwmsDimension):
372    """
373    implements IDimension
374    @ivar units: The units string.
375    @ivar extent: Sequence of extent values.
376
377    """
378   
379    def __init__(self, domain, dimname, unit):
380        self.units = unit
381        self.extent = []
382        #patch to handle current limitations of multiple time dimension scanning in csml.
383        if string.lower(self.units)[:10] in ['days_since', 'seconds_si', 'minutes_si', 'hours_sinc','months _sin', 'years_sinc']:
384            if type(domain[dimname][0]) is not str   :
385                tunits=self.units.replace('_', ' ')
386                for val in domain[dimname]:
387                    csmltime= csml.csmllibs.csmltime.UDtimeToCSMLtime(cdtime.reltime(float(val), tunits).tocomp())
388                    self.extent.append(csmltime)
389            else:
390                for val in domain[dimname]:
391                    self.extent.append(str(val))
392        else:
393            for val in domain[dimname]:
394                self.extent.append(str(val))
395        #for time axis replace units with iso string
396        if dimname == 'time':
397            self.units='ISO8601'
398       
399class CSMLwmsLayerSlab(IwmsLayerSlab):
400    """
401    Implements LayerSlab
402    Represents a particular horizontal slice of a WMS layer.
403
404    ILayerSlab objects are designed to be convenient to cache.
405    They should be pickleable to enable memcached support in the future.
406
407    @ivar layer: The source ILayer instance.
408    @ivar crs: The coordinate reference system.
409    @ivar dimValues: A mapping of dimension values of this view.
410    @ivar renderOpts: The renderOpts used to create this view.
411    @ivar bbox: The bounding box as a 4-tuple.
412    """
413
414    def __init__(self, netcdf, layer, crs, dimValues, renderOpts, bbox):
415        self._netcdf=netcdf
416        self.layer = layer
417        self.crs = crs
418        self.dimValues = dimValues
419        self.renderOpts=renderOpts
420        self.bbox=bbox
421       
422        #set colour map for ALL images from this slab
423        v=self._netcdf(layer.title)
424        tvar=v(squeeze=1)
425        #get the min and max values to use for the colourmapping.
426        #change the fill values to ensure they aren't picked up as false max/mins.
427        tvar.missing_value=999999999
428        value=tvar.getValue()
429        self.minval=min(min(l) for l in value)
430        tvar.missing_value=-0.000000001
431        value=tvar.getValue()
432        self.maxval=max(max(l) for l in value)
433
434       
435    def getImage(self, bbox, width, height):
436        """
437        Create an image of a sub-bbox of a given size.
438
439        @ivar bbox: A bbox 4-tuple.
440        @ivar width: width in pixels.` 
441        @ivar height: height in pixels.
442        @return: A PIL Image object.
443
444        """
445
446        log.debug('getImage(%s, %s, %s)' % (bbox, width, height))
447       
448
449        cmap=eval(config['colourmap']) # renderOpts is hook for colourmap, for now use config
450        grid=Grid(self.layer, self._netcdf, bbox, width, height)
451        #how to handle varmin,varmax? ..read array?
452        #minval, maxval=genutil.minmax(grid.value)
453        #minval=min(min(l) for l in grid.value)
454        #maxval=max(max(l) for l in grid.value)
455        minval=self.minval
456        maxval=self.maxval
457        renderer=RGBARenderer(minval, maxval)         
458        return renderer.renderGrid(grid, bbox, width, height, cmap)
459   
460class Grid(object):
461    """A class encapsulating a simple regularly spaced, rectilinear
462    grid.  This is the only type of grid pywms is expected to
463    understand and adaptors should be provided to connect to
464    underlying implementations such as cdms or csml.
465
466    @cvar crs: Coordinate reference system
467
468    @ivar x0: coordinate of the lower left corner.
469    @ivar y0: coordinate of the lower left corner.
470    @ivar dx: The x grid spacing.
471    @ivar dy: The y grid spacing.
472    @ivar nx: The number of x grid points.
473    @ivar ny: The number of y grid points.
474    @ivar value: A masked array of the grid values.
475    @ivar ix: The dimension index of longidude in value
476    @ivar iy: The dimension index of latitude in value
477    @ivar long_name: The name of the field.
478    @ivar units: The units of the field.
479    """
480    def __init__(self, layer, netcdf, bbox, width, height):
481        #we know the axes are called latitude and longitude as the CSML code has written it:
482        v=netcdf(layer.title)
483        tvar=v(latitude=(bbox[1], bbox[3]), longitude=(bbox[0],bbox[2]),squeeze=1)
484        order=tvar.getOrder()
485        #array of data
486        self.value=tvar.getValue()
487        if order == 'xy':
488            self.ix=0
489            self.iy=1
490        else:
491            self.ix=1
492            self.iy=0
493        lat = tvar.getLatitude()
494        lon = tvar.getLongitude()
495        self.x0=lon[0]
496        self.y0=lat[0]
497        self.dx=abs(lon[0]-lon[1])
498        self.dy=abs(lat[0]-lat[1])
499       
500           
501
502        self.nx=len(lon)
503        self.ny=len(lat)
504        self.long_name=tvar.id  #TODO, get long name from feature
505        self.units=tvar.units
Note: See TracBrowser for help on using the repository browser.