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

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

Dimensions reported properly in capabilities. GetMap? confirmed as working with cdms2

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