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

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

Merging in qesdi changes to cows trunk - still need to merge new backend.

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