source: cows/branches/cows_qesdi/cows/service/imps/csmlbackend/wms_csmllayer.py @ 5493

Subversion URL: http://proj.badc.rl.ac.uk/svn/ndg/cows/branches/cows_qesdi/cows/service/imps/csmlbackend/wms_csmllayer.py@5493
Revision 5493, 19.9 KB checked in by pnorton, 12 years ago (diff)

Implemented the changes made while working on the qesdi wms. Also changed the csmlbackend wms so that it works with the new WMSContoller.

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