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

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

Improved the legend min and max handling in the csml backend, also changed the controller slightly to pass the dimensions to the getLegendImage function.

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