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

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

Removed the old modis data reader class. Also modified COWS to pass the style parameter to the GetLegendImage? methods on the layer objects.

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
30
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                       , style=None):
285        """
286        Create an image of the colourbar for this layer.
287        @param orientation: Either 'vertical' or 'horizontal'
288        @return: A PIL image with labels
289
290        """
291        width = self.legendSize[0]
292        height = self.legendSize[1]
293        # if width and height specified in the GetLegend request parameters
294        # then use them instead of the default values
295        if 'width' in renderOpts:
296                width = renderOpts['width']
297        if 'height' in renderOpts:
298                height = renderOpts['height']
299        cmap = cm.get_cmap(config['colourmap'])
300        renderer=RGBARenderer(self._minval, self._maxval)
301       
302        log.debug("dimValues = %s" % (dimValues,))
303
304        if orientation =='vertical':
305            legendImage= renderer.renderColourbar(630, 30, cmap, isVertical=True)
306        else:
307            legendImage= renderer.renderColourbar(630, 30, cmap, isVertical=False)       
308        imageWithLabels=Image.new('RGBA', (630, 80), "white")
309        imageWithLabels.paste(legendImage, (0,0))
310        #add minvalue label
311        minvalueImg=self._simpletxt2image(str(self._minval), (49,25))
312        imageWithLabels.paste(minvalueImg,(0,40))
313        #add midvalue  label
314        midvalue=self._minval+(self._maxval-self._minval)/2
315        #add maxvalue label
316        midvalueImg=self._simpletxt2image(str(midvalue),(49,25))
317        imageWithLabels.paste(midvalueImg,(280,40))
318        #add maxvalue label
319        maxvalueImg=self._simpletxt2image(str(self._maxval), (49,25))
320        imageWithLabels.paste(maxvalueImg,(575,40))
321        #add units:
322        unitsImg=self._simpletxt2image('Units of measure: %s'%str(self.units), (200,25))
323        imageWithLabels.paste(unitsImg,(240,60))     
324#        return imageWithLabels                         
325        return imageWithLabels.resize((width, height)) 
326   
327    def _simpletxt2image(self, text, size):
328        image = Image.new('RGBA',size,"white")
329        fontfullpath = config['legendfont']
330        ifo = ImageFont.truetype(fontfullpath,16)
331        draw = ImageDraw.Draw(image)
332        draw.text((0, 0), text, font=ifo,fill=(100, 123, 165))
333        return image
334 
335   
336   
337   
338class CSMLwmsDimension(IwmsDimension):
339    """
340    implements IDimension
341    @ivar units: The units string.
342    @ivar extent: Sequence of extent values.
343
344    """
345   
346    def __init__(self, domain, dimname, unit):
347        self.units = unit
348        self.extent = []
349        #patch to handle current limitations of multiple time dimension scanning in csml.
350        if string.lower(self.units)[:10] in ['days_since', 'seconds_si', 'minutes_si', 'hours_sinc','months _sin', 'years_sinc']:
351            if type(domain[dimname][0]) is not str   :
352                tunits=self.units.replace('_', ' ')
353                for val in domain[dimname]:
354                    csmltime= csml.csmllibs.csmltime.UDtimeToCSMLtime(cdtime.reltime(float(val), tunits).tocomp())
355                    self.extent.append(csmltime)
356            else:
357                for val in domain[dimname]:
358                    self.extent.append(str(val))
359        else:
360            for val in domain[dimname]:
361                self.extent.append(str(val))
362        #for time axis replace units with iso string
363        if dimname == 'time':
364            self.units='ISO8601'
365       
366class CSMLwmsLayerSlab(IwmsLayerSlab):
367    """
368    Implements LayerSlab
369    Represents a particular horizontal slice of a WMS layer.
370
371    ILayerSlab objects are designed to be convenient to cache.
372    They should be pickleable to enable memcached support in the future.
373
374    @ivar layer: The source ILayer instance.
375    @ivar crs: The coordinate reference system.
376    @ivar dimValues: A mapping of dimension values of this view.
377    @ivar renderOpts: The renderOpts used to create this view.
378    @ivar bbox: The bounding box as a 4-tuple.
379    """
380
381    def __init__(self, netcdf, layer, crs, dimValues, renderOpts, bbox):
382        self._netcdf=netcdf
383        self.layer = layer
384        self.crs = crs
385        self.dimValues = dimValues
386        self.renderOpts=renderOpts
387        self.bbox=bbox
388       
389        #set colour map for ALL images from this slab
390        v=self._netcdf(layer.title)
391        tvar=v(squeeze=1)
392        #get the min and max values to use for the colourmapping.
393        #change the fill values to ensure they aren't picked up as false max/mins.
394        tvar.missing_value=999999999
395        value=tvar.getValue()
396        self.minval=min(min(l) for l in value)
397        tvar.missing_value=-0.000000001
398        value=tvar.getValue()
399        self.maxval=max(max(l) for l in value)
400
401       
402    def getImage(self, bbox, width, height):
403        """
404        Create an image of a sub-bbox of a given size.
405
406        @ivar bbox: A bbox 4-tuple.
407        @ivar width: width in pixels.` 
408        @ivar height: height in pixels.
409        @return: A PIL Image object.
410
411        """
412
413
414        lbbox = self.layer.getBBox(self.crs)
415        ibbox = bbox_util.intersection(bbox, lbbox)
416   
417        log.debug('bbox = %s' % (bbox,))
418        log.debug('lbbox = %s' % (lbbox,))
419        log.debug('ibbox = %s' % (ibbox,))
420   
421        # If bbox is not within layerObj.bbox then we need to calculate the
422        # pixel offset of the inner bbox, request the right width/height
423        # and paste the image into a blank background
424        if bbox == ibbox:
425            img = self._renderImage(bbox, width, height)
426            log.debug('image.size = %s' % (img.size,))
427                   
428        else:
429           
430            ix0, iy0 = bbox_util.geoToPixel(ibbox[0], ibbox[3], bbox, width, height,
431                                            roundUpY=True)
432            ix1, iy1 = bbox_util.geoToPixel(ibbox[2], ibbox[1], bbox, width, height,
433                                            roundUpX=True)
434            iw = ix1-ix0
435            ih = iy1-iy0
436            log.debug('Deduced inner image: %s, (%d x %d)' % ((ix0, iy0, ix1, iy1), iw, ih))
437            img1 = self._renderImage(ibbox, iw, ih)
438
439            img = Image.new('RGBA', (width, height))
440            img.paste(img1, (ix0, iy0))
441                   
442        return img
443
444
445    def _renderImage(self, bbox, width, height):
446        log.debug('_renderImage(%s, %s, %s)' % (bbox, width, height))
447       
448        cmap = cm.get_cmap(config['colourmap'])
449
450        grid=Grid(self.layer, self._netcdf, bbox, width, height)
451        #If there is no data for the requested area, return a blank image:
452        #TODO this should be included in the overhauled rendering code
453        if grid.ok ==False:
454            width=abs(width)
455            height=abs(height)
456            img = numpy.empty((width,height),numpy.uint32)
457            img.shape=height, width
458            pilImage = Image.frombuffer('RGBA',(width,height),img,'raw','RGBA',0,1)
459            return pilImage
460        #how to handle varmin,varmax? ..read array?
461        #minval, maxval=genutil.minmax(grid.value)
462        #minval=min(min(l) for l in grid.value)
463        #maxval=max(max(l) for l in grid.value)
464        minval=self.minval
465        maxval=self.maxval
466        renderer=RGBARenderer(minval, maxval)         
467        return renderer.renderGrid(grid, bbox, width, height, cmap)
468   
469   
470class Grid(object):
471    """A class encapsulating a simple regularly spaced, rectilinear
472    grid.  This is the only type of grid pywms is expected to
473    understand and adaptors should be provided to connect to
474    underlying implementations such as cdms or csml.
475
476    @cvar crs: Coordinate reference system
477
478    @ivar x0: coordinate of the lower left corner.
479    @ivar y0: coordinate of the lower left corner.
480    @ivar dx: The x grid spacing.
481    @ivar dy: The y grid spacing.
482    @ivar nx: The number of x grid points.
483    @ivar ny: The number of y grid points.
484    @ivar value: A masked array of the grid values.
485    @ivar ix: The dimension index of longidude in value
486    @ivar iy: The dimension index of latitude in value
487    @ivar long_name: The name of the field.
488    @ivar units: The units of the field.
489    """
490    def __init__(self, layer, netcdf, bbox, width, height):
491        #we know the axes are called latitude and longitude as the CSML code has written it:
492        v=netcdf(layer.title)       
493       
494        #sometimes these dimensions get squeezed out so need to get hold of the dx, dy spacing before the variable is subsetted.
495        lat=v.getLatitude()
496        lon=v.getLongitude()
497        self.dx=abs(lon[0]-lon[1])
498        self.dy=abs(lat[0]-lat[1])
499        self.long_name=v.id
500        self.units=v.units
501        #now do the subset.
502        try:
503            tvar=v(latitude=(bbox[1], bbox[3]), longitude=(bbox[0],bbox[2]),squeeze=1)     
504            if type(tvar) == numpy.float32:
505                order ='xy'
506                self.value=numpy.ndarray(tvar)
507                self.ix=0
508                self.iy=1
509                self.x0=bbox[0]
510                self.y0=bbox[1]
511                self.nx=self.ny=1
512            else:    #it's a transient variable.
513                order=tvar.getOrder()
514                #array of data
515                self.value=tvar.getValue()
516                if order == 'xy':
517                    self.ix=0
518                    self.iy=1
519                else:
520                    self.ix=1
521                    self.iy=0
522                lat = tvar.getLatitude()
523                lon = tvar.getLongitude()     
524                if lon is not None:
525                    self.x0=lon[0]
526                    self.nx=len(lon)
527                else: #it's a single longitude value
528                    self.x0=bbox[0]
529                    self.nx= 1     
530                if lat is not None:
531                    self.y0=lat[0]
532                    self.ny=len(lat)
533                else: #it's a single latitude value
534                    self.y0=bbox[1]
535                    self.ny= 1
536                self.ok=True #the grid is correctly initialised
537        except cdms.CDMSError:
538            self.ok=False #the grid failed due to a cdms error (most likely data not available for bounding box)
539           
540       
Note: See TracBrowser for help on using the repository browser.