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

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@3904
Revision 3904, 14.7 KB checked in by domlowe, 12 years ago (diff)

removing prints

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       
216        #cached netcdf is indexed by a tuple of feature id and dimvalues - i.e. typically a particular time and Z value for that feature.
217        #look in dictionary for cached copy, and if so use that as the file object.
218        dictindex=str((self._feature.id, dimValues))
219        if dictindex in self.featureinfofilecache:
220            print 'calling cache'
221            f=self.featureinfofilecache[dictindex]
222        else: #else, use the csml api to subset the feature afresh
223            print 'not calling cache'
224            randomname= csml.csmllibs.csmlextra.getRandomID() + '.nc'
225            result= self._feature.subsetToGridSeries(config['tmpfilebuffer'], ncname=randomname, **dimValues)
226            #for now have to read netcdf back from disk (limitation of CSML api)
227            f=cdms.open(result[1])
228            #append to cache:
229            self.featureinfofilecache[dictindex]=f
230            #and then delete the temporary file
231            os.system('rm %s'%result[1])
232       
233        netcdf = f(self.title)  #netcdf here is a cdms transient variable
234       
235       
236        #Now grab the netCDF object for the point specified.
237        #The reason for the 'cob' option is so that if the grid the data
238        #is defined on does not have a grid point at the point specified,
239        #we should  still get the nearest location
240       
241        t_point = netcdf(latitude=(point[1], point[1], 'cob'), longitude=(point[0], point[0], 'cob'))
242        #now get the value recorded at this location
243        value = t_point.getValue().tolist()
244        print value
245        print t_point.fill_value()
246        #and the fill_value too
247        fill_value = t_point.fill_value()
248        #value is actually embedded in a multi dimensional list,
249        #so we need to extract the actual value from the list
250        while type(value) is list:
251                value = value[0]
252
253        #now check if the value is actually the fill_value rather than
254        #a value recorded at the point specified
255        print [value, fill_value]
256        if (2*fill_value) == value:
257                value = "No value found at position: "+str(point[1])+", "+str(point[0])
258        else:
259                value = "Value found at position: "+str(point[1])+", "+str(point[0])+" is: "+str(value)
260        # finally return the value
261        return value
262
263
264class CSMLDimension(object):
265    """
266    implements IDimension
267    @ivar units: The units string.
268    @ivar extent: Sequence of extent values.
269
270    """
271   
272    def __init__(self, domain, dimname, unit):
273        self.units = unit
274        self.extent = []
275        #patch to handle current limitations of multiple time dimension scanning in csml.
276        if string.lower(self.units)[:10] in ['days_since', 'seconds_si', 'minutes_si', 'hours_sinc','months _sin', 'years_sinc']:
277            if type(domain[dimname][0]) is not str   :
278                tunits=self.units.replace('_', ' ')
279                for val in domain[dimname]:
280                    csmltime= csml.csmllibs.csmltime.UDtimeToCSMLtime(cdtime.reltime(float(val), tunits).tocomp())
281                    self.extent.append(csmltime)
282            else:
283                for val in domain[dimname]:
284                    self.extent.append(str(val))
285        else:
286            for val in domain[dimname]:
287                self.extent.append(str(val))
288        #for time axis replace units with iso string
289        if dimname == 'time':
290            self.units='ISO8601'
291       
292class CSMLLayerSlab(object):
293    """
294    Implements LayerSlab
295    Represents a particular horizontal slice of a WMS layer.
296
297    ILayerSlab objects are designed to be convenient to cache.
298    They should be pickleable to enable memcached support in the future.
299
300    @ivar layer: The source ILayer instance.
301    @ivar crs: The coordinate reference system.
302    @ivar dimValues: A mapping of dimension values of this view.
303    @ivar renderOpts: The renderOpts used to create this view.
304    @ivar bbox: The bounding box as a 4-tuple.
305    """
306   
307    def __init__(self, netcdf, layer, crs, dimValues, renderOpts, bbox):
308        self._netcdf=netcdf
309        self.layer = layer
310        self.crs = crs
311        self.dimValues = dimValues
312        self.renderOpts=renderOpts
313        self.bbox=bbox
314       
315        #set colour map for ALL images from this slab
316        v=self._netcdf(layer.title)
317        tvar=v(squeeze=1)
318        #array of data
319        value=tvar.getValue()
320        self.minval=min(min(l) for l in value)
321        self.maxval=max(max(l) for l in value)
322       
323       
324    def getImage(self, bbox, width, height):
325        """
326        Create an image of a sub-bbox of a given size.
327
328        @ivar bbox: A bbox 4-tuple.
329        @ivar width: width in pixels.` 
330        @ivar height: height in pixels.
331        @return: A PIL Image object.
332
333        """
334
335        log.debug('getImage(%s, %s, %s)' % (bbox, width, height))
336       
337
338        cmap=eval(config['colourmap']) # renderOpts is hook for colourmap, for now use config
339        grid=Grid(self.layer, self._netcdf, bbox, width, height)
340        #how to handle varmin,varmax? ..read array?
341        #minval, maxval=genutil.minmax(grid.value)
342        #minval=min(min(l) for l in grid.value)
343        #maxval=max(max(l) for l in grid.value)
344        minval=self.minval
345        maxval=self.maxval
346        renderer=RGBARenderer(minval, maxval)         
347        return renderer.renderGrid(grid, bbox, width, height, cmap)
348   
349class Grid(object):
350    """A class encapsulating a simple regularly spaced, rectilinear
351    grid.  This is the only type of grid pywms is expected to
352    understand and adaptors should be provided to connect to
353    underlying implementations such as cdms or csml.
354
355    @cvar crs: Coordinate reference system
356
357    @ivar x0: coordinate of the lower left corner.
358    @ivar y0: coordinate of the lower left corner.
359    @ivar dx: The x grid spacing.
360    @ivar dy: The y grid spacing.
361    @ivar nx: The number of x grid points.
362    @ivar ny: The number of y grid points.
363    @ivar value: A masked array of the grid values.
364    @ivar ix: The dimension index of longidude in value
365    @ivar iy: The dimension index of latitude in value
366    @ivar long_name: The name of the field.
367    @ivar units: The units of the field.
368    """
369    def __init__(self, layer, netcdf, bbox, width, height):
370        #we know the axes are called latitude and longitude as the CSML code has written it:
371        #print netcdf
372        v=netcdf(layer.title)
373        #print v
374        tvar=v(latitude=(bbox[1], bbox[3]), longitude=(bbox[0],bbox[2]),squeeze=1)
375        order=tvar.getOrder()
376        #array of data
377        self.value=tvar.getValue()
378        #print self.value
379        #order of axes
380        if order == 'xy':
381            self.ix=0
382            self.iy=1
383        else:
384            self.ix=1
385            self.iy=0
386        lat = tvar.getLatitude()
387        lon = tvar.getLongitude()
388        self.x0=lon[0]
389        self.y0=lat[0]
390        self.dx=abs(lon[0]-lon[1])
391        self.dy=abs(lat[0]-lat[1])
392        #print [lon, len(lon)]
393        #print len(self.value)
394        print len(lon)
395        print len(lat)
396
397   
398
399       
400
401        self.nx=len(lon)
402        self.ny=len(lat)
403        self.long_name=tvar.id  #TODO, get long name from feature
404        self.units=tvar.units
Note: See TracBrowser for help on using the repository browser.