source: cows/branches/cows-vis/cows/service/imps/csmlbackend/wms_csmllayer.py @ 5265

Subversion URL: http://proj.badc.rl.ac.uk/svn/ndg/cows/branches/cows-vis/cows/service/imps/csmlbackend/wms_csmllayer.py@5265
Revision 5265, 16.2 KB checked in by domlowe, 11 years ago (diff)

adding separate cowsclient pylons app

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