1 | """ |
---|
2 | implementation of ILayerMapper, ILayer, IDimension, ILayerSlab interfaces, as defined in wms_iface.py |
---|
3 | |
---|
4 | """ |
---|
5 | import os, string |
---|
6 | import csml |
---|
7 | try: |
---|
8 | import cdms2 as cdms |
---|
9 | except: |
---|
10 | import cdms |
---|
11 | import Image |
---|
12 | from copy import copy |
---|
13 | from pywms.render_imp import RGBARenderer |
---|
14 | from matplotlib import cm |
---|
15 | import genutil |
---|
16 | from pylons import config #config must have tmpfilebuffer and csmlstore values |
---|
17 | import cdtime |
---|
18 | |
---|
19 | import logging |
---|
20 | log = logging.getLogger(__name__) |
---|
21 | |
---|
22 | class 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 | |
---|
98 | class 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 | |
---|
269 | class 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 | |
---|
297 | class 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 | |
---|
354 | class 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 |
---|