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

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

Splitting CSML backend into common and service specific elements.

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