source: qesdi/geoplot/trunk/lib/geoplot/grid_builder_national.py @ 5826

Subversion URL: http://proj.badc.rl.ac.uk/svn/ndg/qesdi/geoplot/trunk/lib/geoplot/grid_builder_national.py@5826
Revision 5826, 7.5 KB checked in by pnorton, 11 years ago (diff)

Improved geoplot's behaviour when dealing with variables with axis in the order of lon/lat instead of lat/lon.

Line 
1"""
2grid_builder_national.py
3============
4
5A GridBuilder object that can extract a lat-lon grid form a cdms vairable
6conating a grid in terms of british national grid coordinates.
7
8"""
9
10#python modules
11import logging
12from math import *
13
14#third party modules
15from geoplot.mpl_imports import basemap
16
17#internal modules
18import numpy as N
19
20from geoplot.grid_builder_base import GridBuilderBase
21
22#setup the logging
23log = logging.getLogger(__name__)
24
25# National Grid parameters
26trueOrigin = (358.0, 49.0)
27falseOrigin = (-400000, 100000)
28
29maxXBounds=(-20.0,  15.0); maxYBounds=( 40.0, 70.0)
30
31# Create a basemap with the origin at the true origin.
32nationalGrid = basemap.Basemap(projection='tmerc',
33                               lon_0=trueOrigin[0],
34                               lat_0=trueOrigin[1],
35                               llcrnrlon=trueOrigin[0], llcrnrlat=trueOrigin[1],
36                               urcrnrlon=10, urcrnrlat=10,
37                               resolution=None)
38
39class GridBuilderNational(GridBuilderBase):
40
41    def __init__(self, cdmsVar):
42        GridBuilderBase.__init__(self, cdmsVar)
43
44    def _resizeVar(self, xLimits, yLimits):
45        """
46        returns a cdms variable that contains the data from self.cdmsVar that
47        falls within the lat-lon limits provided.
48        """
49        log.debug('Resizing grid, xLimits: ' + repr(xLimits) + ' yLimits: ' + repr(yLimits))
50
51        if xLimits[0] < maxXBounds[0] or xLimits[1] > maxXBounds[1]:
52            raise Exception("xLimits (=%s) is outside maximum %s" % (xLimits, maxXBounds))
53       
54        if yLimits[0] < maxYBounds[0] or yLimits[1] > maxYBounds[1]:
55            raise Exception("yLimits (=%s) is outside maximum %s" % (yLimits, maxYBounds))
56       
57
58        eastLow, northLow = nationalGrid(xLimits[0], yLimits[0])
59        eastHigh, northHigh = nationalGrid(xLimits[1], yLimits[1])
60
61        #need to add the false origin to shift from map units to British national grid
62        eastLow = eastLow - falseOrigin[0]
63        eastHigh = eastHigh - falseOrigin[0]
64        northLow = northLow - falseOrigin[1]
65        northHigh = northHigh - falseOrigin[1]
66
67        log.debug("eastings selection = %s, %s" % (eastLow, eastHigh,))
68        log.debug("northings selection = %s, %s" % (northLow, northHigh,))
69
70        eastgingVals =  self.cdmsVar.getAxis(self.cdmsVar.getAxisIndex('eastings'))[:]
71        northingVals = self.cdmsVar.getAxis(self.cdmsVar.getAxisIndex('northings'))[:]
72       
73        eastingsLowHigh = eastgingVals[0] <= eastgingVals[-1]
74        northingsLowHigh = northingVals[0] <= northingVals[-1]
75           
76        log.debug("eastingsLowHigh = %s" % (eastingsLowHigh,))
77        log.debug("northingsLowHigh = %s" % (northingsLowHigh,))
78       
79        if eastingsLowHigh:
80            eastingsSelection = (eastLow, eastHigh)
81        else:
82            eastingsSelection = (eastHigh, eastLow)       
83       
84        if northingsLowHigh:
85            northingsSelection = (northLow, northHigh)
86        else:
87            northingsSelection = (northHigh, northLow)
88       
89        log.debug("selecting eastingsSelection %s" % (eastingsSelection,))
90        log.debug("selecting northingsSelection %s" % (northingsSelection,))
91       
92        newVar = None
93        try:
94            newVar = self.cdmsVar(northings=northingsSelection, eastings=eastingsSelection)
95        except:
96           
97            log.exception("Exception occurred while trying to subset variable "
98                          "with northings=%s and eastings=%s. Using original variable." \
99                          % (northingsSelection, eastingsSelection))
100            newVar = self.cdmsVar           
101           
102        else:
103#            log.debug("newVar = %s" % (newVar,))
104            log.debug("newVar.getAxisIds()= %s" % (newVar.getAxisIds(),))       
105
106
107        return newVar
108   
109    def _buildGridBounds(self, cdmsVar):
110        """
111        Builds two 2d Numpy array of the x and y positions of the grid
112        boundaries.
113       
114        Returns one array for the x bounds and one for the y bound. These
115        arrays contain the lower boundary for a given grid box.
116       
117        e.g. xBounds[x,y] will give the lower x boundary for the gridbox x,y.
118        """
119       
120        ids = cdmsVar.getAxisIds()
121        log.debug("ids = %s" % (ids,))
122       
123        assert ids == ['northings', 'eastings']
124       
125        yAxis = cdmsVar.getAxis(0)
126        xAxis = cdmsVar.getAxis(1)
127       
128        xBounds = GridBuilderBase._getBoundsFromAxis(xAxis)
129        yBounds = GridBuilderBase._getBoundsFromAxis(yAxis)
130
131        xMesh, yMesh = N.meshgrid(xBounds, yBounds)
132
133        #transform from national grid into lat lon
134        #log.debug("BEFORE\n" + str(self.xMesh) +"\n"+ str(self.yMesh))
135
136        #remove the false origin before undoing the transform
137        xMesh = xMesh + falseOrigin[0]
138        yMesh = yMesh + falseOrigin[1]
139
140        xMesh , yMesh  = nationalGrid(xMesh, yMesh, inverse=True)
141        #log.debug("AFTER\n" + str(self.xMesh) +"\n"+ str(self.yMesh))
142
143        return (xMesh, yMesh)
144
145    def _buildGridMidpoints(self, cdmsVar):
146        """
147        Builds two 2d Numpy arrays for the x and y midpoints of a given grid
148        box, this position corresponds to the location of the grid box value.
149       
150        Returns one array for the midpoint x position and another for the y
151        postion.
152       
153        e.g. xMidpoints[x,y] will give the midpoint of the gridbox x,y and the
154        position of the ascociated measurment (value[x,y]).
155        """
156
157        yAxis = cdmsVar.getAxis(0)
158        xAxis = cdmsVar.getAxis(1)
159
160        yPos = yAxis.getValue()
161        xPos = xAxis.getValue()
162
163        xAxis, yAxis = N.meshgrid(xPos, yPos)
164       
165        #remove the false origin before undoing the transform
166        xAxis = xAxis + falseOrigin[0]
167        yAxis = yAxis + falseOrigin[1]
168
169        xAxis , yAxis  = nationalGrid(xAxis, yAxis, inverse=True)
170       
171        return (xAxis, yAxis)
172
173    @staticmethod 
174    def _areAxisInOrderYX(cdmsVar):
175        #!TODO: implement this method
176        return True
177   
178if __name__ == '__main__':
179
180    import os
181    import pkg_resources
182
183    outputsDir = os.path.abspath(pkg_resources.resource_filename('geoplot', '../../outputs'))
184    testDataDir = pkg_resources.resource_filename('geoplot', 'tests/data')
185
186    import log_util
187    log_util.setupGeoplotConsoleHandler(log)
188   
189    dataFile = testDataDir + '/ng_test.nc'
190    dataVariable = 'mean_temp'
191    f = cdms.open(dataFile)
192    var = f(dataVariable, squeeze=1)
193   
194    builder = GridBuilderNational(var)
195
196    xLimits=(-13.0, 6.2) ; yLimits=(47.0, 61.0)
197
198    grid = builder.buildGrid(xLimits, yLimits)
199
200    from matplotlib.toolkits import basemap
201    from matplotlib.backends.backend_agg import FigureCanvasAgg
202    from matplotlib.figure import Figure
203    from matplotlib import cm
204
205    width=600;height=800;dpi=80
206    figsize=(width / dpi, height / dpi)
207    fig = Figure(figsize=figsize, dpi=dpi)
208    axes = fig.add_axes([0.1, 0.1, 0.8, 0.8])
209
210
211    bm = basemap.Basemap(llcrnrlon=xLimits[0],
212                         llcrnrlat=yLimits[0],
213                         urcrnrlon=xLimits[1],
214                         urcrnrlat=yLimits[1],
215                         resolution='h',
216                         suppress_ticks=False)
217    cmap = cm.jet
218    cmap.set_bad("w")
219
220    sm = bm.pcolormesh(grid.boundsX,grid.boundsY,grid.values,
221                   cmap=cmap, ax=axes, edgecolors='k')
222
223    bm.drawcoastlines(ax = axes)
224    bm.drawrivers(ax = axes, color='b')
225
226    canvas = FigureCanvasAgg(fig)
227    filename = "out.png"
228    canvas.print_figure(filename)
229    log.debug("Wrote " + filename)
Note: See TracBrowser for help on using the repository browser.