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

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

csml wms layer returns empty image if there is no data for a region (avoids pink tiles))

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