source: TI05-delivery/ows_framework/branches/ows_framework-refactor/ows_common/ows_common/service/imps/wms_csmllayer.py @ 3915

Subversion URL: http://proj.badc.rl.ac.uk/svn/ndg/TI05-delivery/ows_framework/branches/ows_framework-refactor/ows_common/ows_common/service/imps/wms_csmllayer.py@3915
Revision 3915, 15.8 KB checked in by domlowe, 11 years ago (diff)

Added support to WMS for retrieving CSML docs from eXist

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