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

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

Moving to a different interface between the WMS controller and the
data. interfaces.py now contains ILayer et al. to represent layers
of data.

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 self.var.getTime():
126            self.dimensions = dict(time=CdmsTimeDimension(self.var.getTime()))
127        else:
128            self.dimensions = {}
129
130        if minValue is not None:
131            self.minValue = minValue
132        else:
133            self.minValue = self.var.min_value
134
135        if maxValue is not None:
136            self.maxValue = maxValue
137        else:
138            self.maxValue = self.var.max_value
139
140        try:
141            self.units = self.var.units
142        except AttributeError:
143            self.units = '?'
144
145    def selectGrid(self, bbox, dimensionSpec):
146        (lon1, lat1, lon2, lat2) = bbox
147        sel = dict(latitude=(lat1, lat2, 'cce'), longitude=(lon1, lon2, 'cce'))
148        if 'time' in self.dimensions:
149            sel['time'] = self.dimensions['time'].iso2reltime(dimensionSpec['time'])
150        sel['squeeze'] = 1
151           
152        v = self.var(**sel)
153        return self.GridClass(v, title=self.title)
154
155    def describe(self, dimensionSpec):
156        return self.var.long_name + ' at ' + self.dimensions['time']
157
158
159
160class CdmsTimeDimension(object):
161    """
162    @todo: Move to impl.py when interface migration complete.
163
164    """
165    def __init__(self, timeAxis):
166        self._axis = timeAxis
167        self.units = 'ISO 8601'
168
169    def _getExtent(self):
170        comptimes = self._axis.asComponentTime()
171        return ','.join(['%s-%s-%sT%s:%s:%sZ' % (x.year, x.month, x.day,
172                                                 x.hour, x.minute, x.second)
173                         for x in comptimes])
174    extent = property(_getExtent)
175
176    def iso2reltime(self, time, yearOverride=None):
177        mo = re.match(r'(\d+)-(\d+)-(\d+)T(\d+):(\d+):([0-9.]+)Z', time)
178        if not mo:
179            raise InvalidParameterValue('Time %s not recognised' % time, 'time')
180        (year, month, day, hour, minute, second) = mo.groups()
181        if yearOverride:
182            year = yearOverride
183        c = cdtime.comptime(int(year), int(month), int(day),
184                            int(hour), int(minute), float(second))
185        return c.torel(self._axis.units)
186
187    def iso2timeDelta(self, time):
188        return [self.iso2reltime(x) for x in time.split('D')]
189
190
191
192
193
194#-----------------------------------------------------------------------------
195
196import unittest
197
198class TestBNGGrid(unittest.TestCase):
199    def setUp(self):
200        import os
201
202        try:
203            self.data_dir = os.environ['TEST_DATA_DIR']
204            self.rainfall = cdms.open(self.data_dir+'/ukcip02/rainfall_1961-2000.nc')
205        except:
206            raise RuntimeError("""
207Test data not found.  Please set the TEST_DATA_DIR environment variable.
208""")
209
210        # wire in the render_imp logger to nose
211        import render_imp
212        render_imp.logger = render_imp.logging.getLogger('nose.render_imp')
213 
214
215    def _makeGrid(self):
216        v = self.rainfall['rainfall'][0]
217        self.assertEquals(v.id, 'rainfall')
218        return CdmsBNGGrid(v, 'UKCIP02 rainfall data')
219
220    def testGridAttributes(self):
221        g = self._makeGrid()
222
223        self.assertEquals(g.dx, 5000)
224        assert g.dy == -5000
225        self.assertEquals(g.nx, 180)
226        assert g.ny == 290
227        assert g.x0 == -200000
228        assert g.y0 == -200000
229        self.assertNotEquals(g.ix, g.iy)
230
231    def testGridValue(self):
232        g = self._makeGrid()
233        v = g.value
234
235        # Value should be masked
236        assert v.mask()
237
238    def testRender(self):
239        from render_imp import RGBARenderer
240        from matplotlib.cm import get_cmap
241
242        g = self._makeGrid()
243        # Set arbitary min/max values for now
244        r = RGBARenderer(MA.minimum(g.value), MA.maximum(g.value))
245
246        xn = g.x0+g.dx*g.nx
247        yn = g.y0+g.dy*g.ny
248        bbox = (min(g.x0, xn), min(g.y0, yn),
249                max(g.x0, xn), max(g.y0, yn))
250
251        img = r.renderGrid(g, bbox, 400, 400, get_cmap())
252        img.save('whole_domain.png')
253
254        assert img.size == (400, 400)
255
256        bbox2 = (bbox[0], bbox[1], bbox[0] + (bbox[2]-bbox[0])/2,
257                 bbox[1] + (bbox[3]-bbox[1])/2)
258        img = r.renderGrid(g, bbox2, 400, 400, get_cmap())
259        img.save('ll_quad.png')
260
261       
Note: See TracBrowser for help on using the repository browser.