source: cows/trunk/cows/service/imps/wms_csmllayer.py @ 4228

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

Changed print statements to loggging calls where appropriate.

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