source: cows/trunk/cows/service/imps/csmlbackend/wms_csmllayer.py @ 5365

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

Changes to cows csml backend to render wms legend bar.

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