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

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@3903
Revision 3903, 14.8 KB checked in by domlowe, 11 years ago (diff)

csml layer caching added at the layermapper level and at the getfeatureinfo level - this has significantly speeded up the wms

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