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

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

Fixed a pywms class name error. cows.test.util creates a test CSML
dataset.

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