source: DPPP/ukcip02_server/trunk/ukcip_server/ukcip_server/model/pywms/wms_cdms.py @ 3533

Subversion URL: http://proj.badc.rl.ac.uk/svn/ndg/DPPP/ukcip02_server/trunk/ukcip_server/ukcip_server/model/pywms/wms_cdms.py@3533
Revision 3533, 7.2 KB checked in by spascoe, 13 years ago (diff)

New interface works in view.html

Line 
1"""
2An implementation of model.WebMapService drived from a single CDMS
3file, potentially a CDML aggregation file.
4"""
5
6from model import Layer, Grid
7import cdtime, cdutil, re
8import cdms2 as cdms
9from ows_common.exceptions import InvalidParameterValue
10import numpy.oldnumeric.ma as MA
11import numpy.oldnumeric as N
12
13class GridError(Exception):
14    pass
15
16class CdmsGrid(Grid):
17    """An adaptor for a cdms variable.
18    """
19
20    def __init__(self, var, title=None):
21        """
22        @param var: The cdms variable.  This should be of shape (x, lat).
23        """
24
25        try:
26            self._var = var
27            self._setGrid()
28        except GridError, e:
29            # This isn't a simple grid
30            #!TODO: we could regrid here.
31            raise e
32        self._setMetadata(title=title)
33
34
35    def _setGrid(self):
36        """Check the grid is simple and initialise.
37        """
38
39        y = self.getYAxis()
40        dy_a = y[1:] - y[:-1]
41        if not (dy_a == (N.ones(len(dy_a)) * dy_a[0])).all():
42            raise SimpleGridError, "Y Axis not equally spaced"
43        self.y0 = y[0]
44        self.dy = dy_a[0]
45        self.ny = len(y)
46        self.iy = self._var.getAxisList().index(y)
47
48        x = self.getXAxis()
49        dx_a = x[1:] - x[:-1]
50        if not (dx_a == (N.ones(len(dx_a)) * dx_a[0])).all():
51            raise SimpleGridError, "X Axis not equally spaced"
52        self.x0 = x[0]
53        self.dx = dx_a[0]
54        self.nx = len(x)
55        self.ix = self._var.getAxisList().index(x)
56
57
58    def _setMetadata(self, title=None):
59        if title:
60            self.long_name = title
61        else:
62            try:
63                self.long_name = self._var.long_name
64            except AttributeError:
65                self.long_name = 'Unknown'
66
67        try:
68            self.units = self._var.units
69        except AttributeError:
70            self.units = ''
71
72    def _getValue(self):
73        return self._var
74    value = property(_getValue)
75
76
77class CdmsLatLonGrid(CdmsGrid):
78    """
79    Specialise CdmsGrid for EPSG:4326
80
81    """
82
83    crs = 'EPSG:4326'
84
85    def getXAxis(self):
86        return self._var.getLongitude()
87    def getYAxis(self):
88        return self._var.getLatitude()
89
90class CdmsBNGGrid(CdmsGrid):
91    """
92    Specialise CdmsGrid for British National Grid coordinate system
93
94    EPSG:27700 OSGB:36
95
96    """
97
98    crs = 'EPSG:27700'
99   
100    def getXAxis(self):
101        return self._var.getAxisList()[1]
102    def getYAxis(self):
103        return self._var.getAxisList()[0]
104
105
106
107class SimpleCdmsLayer(Layer):
108    def __init__(self, cdmsVar, minValue=None, maxValue=None, GridClass=CdmsGrid):
109        """
110        @param cdmsVar: variable object
111
112        @todo: Add crs attribute
113        """
114
115        self.GridClass = GridClass
116
117        self.var = cdmsVar
118        try:
119            long_name = self.var.long_name
120        except AttributeError:
121            long_name = self.var.id
122
123        super(SimpleCdmsLayer, self).__init__(long_name)
124
125        #if 'time' in self.var.axes:
126        #    self.dimensions = dict(time=CdmsTimeDimension(self.var.axes['time']))
127        if self.var.getTime():
128            self.dimensions = dict(time=CdmsTimeDimension(self.var.getTime()))
129        else:
130            self.dimensions = {}
131
132##         if minValue is not None:
133##             self.minValue = minValue
134##         else:
135##             self.minValue = self.var.min_value
136
137##         if maxValue is not None:
138##             self.maxValue = maxValue
139##         else:
140##             self.maxValue = self.var.max_value
141
142        self.minValue = minValue
143        self.maxValue = maxValue
144
145        try:
146            self.units = self.var.units
147        except AttributeError:
148            self.units = '?'
149
150    def selectGrid(self, bbox, dimensionSpec):
151        (lon1, lat1, lon2, lat2) = bbox
152        sel = dict(latitude=(lat1, lat2, 'cce'), longitude=(lon1, lon2, 'cce'))
153        if 'time' in self.dimensions:
154            sel['time'] = self.dimensions['time'].iso2reltime(dimensionSpec['time'])
155        sel['squeeze'] = 1
156           
157        v = self.var(**sel)
158        return self.GridClass(v, title=self.title)
159
160    def describe(self, dimensionSpec):
161        return self.var.long_name + ' at ' + self.dimensions['time']
162
163
164
165class CdmsTimeDimension(object):
166    """
167    @todo: Move to impl.py when interface migration complete.
168
169    """
170    def __init__(self, timeAxis):
171        self._axis = timeAxis
172        self.units = 'ISO 8601'
173
174    def _getExtent(self):
175        comptimes = self._axis.asComponentTime()
176        return ','.join(['%s-%s-%sT%s:%s:%sZ' % (x.year, x.month, x.day,
177                                                 x.hour, x.minute, x.second)
178                         for x in comptimes])
179    extent = property(_getExtent)
180
181    def iso2reltime(self, time, yearOverride=None):
182        mo = re.match(r'(\d+)-(\d+)-(\d+)T(\d+):(\d+):([0-9.]+)Z', time)
183        if not mo:
184            raise InvalidParameterValue('Time %s not recognised' % time, 'time')
185        (year, month, day, hour, minute, second) = mo.groups()
186        if yearOverride:
187            year = yearOverride
188        c = cdtime.comptime(int(year), int(month), int(day),
189                            int(hour), int(minute), float(second))
190        return c.torel(self._axis.units)
191
192    def iso2timeDelta(self, time):
193        return [self.iso2reltime(x) for x in time.split('D')]
194
195
196
197
198
199#-----------------------------------------------------------------------------
200
201import unittest
202
203class TestBNGGrid(unittest.TestCase):
204    def setUp(self):
205        import os
206
207        try:
208            self.data_dir = os.environ['TEST_DATA_DIR']
209            self.rainfall = cdms.open(self.data_dir+'/ukcip02/rainfall_1961-2000.nc')
210        except:
211            raise RuntimeError("""
212Test data not found.  Please set the TEST_DATA_DIR environment variable.
213""")
214
215        # wire in the render_imp logger to nose
216        import render_imp
217        render_imp.logger = render_imp.logging.getLogger('nose.render_imp')
218 
219
220    def _makeGrid(self):
221        v = self.rainfall['rainfall'][0]
222        self.assertEquals(v.id, 'rainfall')
223        return CdmsBNGGrid(v, 'UKCIP02 rainfall data')
224
225    def testGridAttributes(self):
226        g = self._makeGrid()
227
228        self.assertEquals(g.dx, 5000)
229        assert g.dy == -5000
230        self.assertEquals(g.nx, 180)
231        assert g.ny == 290
232        assert g.x0 == -200000
233        assert g.y0 == -200000
234        self.assertNotEquals(g.ix, g.iy)
235
236    def testGridValue(self):
237        g = self._makeGrid()
238        v = g.value
239
240        # Value should be masked
241        assert v.mask()
242
243    def testRender(self):
244        from render_imp import RGBARenderer
245        from matplotlib.cm import get_cmap
246
247        g = self._makeGrid()
248        # Set arbitary min/max values for now
249        r = RGBARenderer(MA.minimum(g.value), MA.maximum(g.value))
250
251        xn = g.x0+g.dx*g.nx
252        yn = g.y0+g.dy*g.ny
253        bbox = (min(g.x0, xn), min(g.y0, yn),
254                max(g.x0, xn), max(g.y0, yn))
255
256        img = r.renderGrid(g, bbox, 400, 400, get_cmap())
257        img.save('whole_domain.png')
258
259        assert img.size == (400, 400)
260
261        bbox2 = (bbox[0], bbox[1], bbox[0] + (bbox[2]-bbox[0])/2,
262                 bbox[1] + (bbox[3]-bbox[1])/2)
263        img = r.renderGrid(g, bbox2, 400, 400, get_cmap())
264        img.save('ll_quad.png')
265
266       
Note: See TracBrowser for help on using the repository browser.