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

Subversion URL: http://proj.badc.rl.ac.uk/svn/ndg/qesdi/geoplot/trunk/lib/geoplot/grid_builder_national.py@7577
Revision 7577, 7.9 KB checked in by astephen, 10 years ago (diff)

Committed changes to be compatible with Matplotlib and Basemap versions
1.0.0. Note that these changes "should" be backwards compatible.
Egg created from this version.

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