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

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

Splitting wms_iface into wxs_iface and wms_iface interfaces.
Modifying implemenations of wms_iface to inherit from new interfaces.

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