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

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

various small changes to csml backend

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