source: qesdi/geoplot/trunk/lib/geoplot/grid_builder_lat_lon.py @ 5890

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

Fixed a problem with the subsetting in the grid_builder_lat_lon_subset.py file.

Line 
1"""
2grid_builder_lat_lon.py
3============
4
5A GridBuilder object that can extract a lat-lon grid from a cdmsVariable. The
6variable must have a latitude and a longitude axis and no level or time data.
7
8"""
9
10#python modules
11import logging
12
13#third party modules
14import numpy as N
15
16#internal modules
17from geoplot.grid_builder_base import GridBuilderBase
18
19#set the log
20log = logging.getLogger(__name__)
21
22class GridBuilderLatLon(GridBuilderBase):
23    """
24    A grid builder that knows how to extract a lat-long grid from a cdms file.
25    """
26
27    def __init__(self, cdmsVar):
28        """
29        Constructs a GridBuilderLatLon object
30        """
31        GridBuilderBase.__init__(self, cdmsVar)
32
33    def _resizeVar(self, xLimits, yLimits):
34        """
35        returns a cdms variable that contains the data from self.cdmsVar that
36        falls within the lat-lon limits provided.
37       
38        xLimits and yLimits must be of the form (low, high)
39        """
40       
41        #check that the longitude is circular
42       
43        if not self.cdmsVar.getLongitude().isCircular():
44            log.warning("Longitude axis of variable %s isn't circular, this may cause problems subsetting." % (self.cdmsVar.id,))
45            log.warning("self.cdmsVar.getLongitude().isCircular() = %s" % (self.cdmsVar.getLongitude().isCircular(),))
46            log.warning("self.cdmsVar.getLongitude().isCircularAxis() = %s" % (self.cdmsVar.getLongitude().isCircularAxis(),))
47       
48#        log.debug('Resizing grid, xLimits: ' + repr(xLimits) + ' yLimits: ' + repr(yLimits))
49
50        lonVals =  self.cdmsVar.getLongitude()[:]
51        latVals = self.cdmsVar.getLatitude()[:]
52       
53#        log.debug("Starting longitude from %s to %s" % ( self.cdmsVar.getLongitude().getValue()[0], self.cdmsVar.getLongitude().getValue()[-1]))
54#        log.debug("Starting latitude  from %s to %s" % ( self.cdmsVar.getLatitude().getValue()[0] , self.cdmsVar.getLatitude().getValue()[-1]))
55       
56        lonLowHigh = lonVals[0] <= lonVals[-1]
57        latLowHigh = latVals[0] <= latVals[-1]
58       
59        # want to keep some extra data around the edges to make sure data
60        # is drawn right to the edge.
61        indicator = 'cce'
62               
63        if lonLowHigh:
64            lonSelection = (xLimits[0], xLimits[1], indicator)
65        else:
66            lonSelection = (xLimits[1], xLimits[0], indicator)           
67       
68        if latLowHigh:
69            latSelection = (yLimits[0], yLimits[1], indicator)
70        else:
71            latSelection = (yLimits[1], yLimits[0], indicator)
72       
73        #!TODO: should implement a fix like this that can be turned on as needed
74        #lonSelection = (lonSelection[0] - 5, lonSelection[1] + 5)
75        #latSelection = (latSelection[0] - 5, latSelection[1] + 5)
76       
77        log.debug("selecting longitude=%s latitude=%s" % (lonSelection, latSelection))
78                           
79        try:
80            newVar = self.cdmsVar(longitude=lonSelection,
81                                  latitude =latSelection)
82        except: 
83            log.exception("Exception occurred while subsetting variable")
84            newVar = self.cdmsVar
85       
86#        log.debug("New longitude from %s to %s" % ( newVar.getLongitude().getValue()[0], newVar.getLongitude().getValue()[-1]))
87#        log.debug("New latitude  from %s to %s" % ( newVar.getLatitude().getValue()[0] , newVar.getLatitude().getValue()[-1]))
88
89        return newVar
90
91
92    def _buildGridBounds(self, cdmsVar):
93        """
94        Builds two 2d Numpy array of the x and y positions of the grid
95        boundaries.
96       
97        Returns one array for the x bounds and one for the y bound. These
98        arrays contain the lower boundary for a given grid box.
99       
100        e.g. xBounds[x,y] will give the lower x boundary for the gridbox x,y.
101        """
102
103        lon = cdmsVar.getLongitude()
104        lat = cdmsVar.getLatitude()
105       
106        latBoundsMerged = GridBuilderBase._getBoundsFromAxis(lat)
107        lonBoundsMerged = GridBuilderBase._getBoundsFromAxis(lon)
108       
109        gridBoundsX, gridBoundsY = N.meshgrid(lonBoundsMerged, latBoundsMerged,)
110       
111        return (gridBoundsX, gridBoundsY)
112   
113    def _buildGridMidpoints(self, cdmsVar):
114        """
115        Builds two 2d Numpy arrays for the x and y midpoints of a given grid
116        box, this position corresponds to the location of the grid box value.
117       
118        Returns one array for the midpoint x position and another for the y
119        postion.
120       
121        e.g. xMidpoints[x,y] will give the midpoint of the gridbox x,y and the
122        position of the ascociated measurment (value[x,y]).
123        """
124
125        lon = cdmsVar.getLongitude()
126        lat = cdmsVar.getLatitude()
127
128        #need to create two arays of the one dimensional array repeated that have the same size
129        gridMidpointX, gridMidpointY = N.meshgrid(lon.getValue(), lat.getValue())
130       
131        return (gridMidpointX, gridMidpointY)
132
133    @staticmethod
134    def _areAxisInOrderYX(cdmsVar):
135        "Check if the axis on a given variable are longitude/latitude or not"
136
137        lon = cdmsVar.getLongitude()
138        lat = cdmsVar.getLatitude()
139       
140        if cdmsVar.getAxisIndex(lon) < cdmsVar.getAxisIndex(lat):
141            orderIsLatLon = False
142        else:
143            orderIsLatLon = True
144       
145        return orderIsLatLon
146
147if __name__ == '__main__':
148
149    import os
150    import pkg_resources
151    import cdms2 as cdms
152    from geoplot.mpl_imports import basemap
153    from matplotlib.backends.backend_agg import FigureCanvasAgg
154    from matplotlib.figure import Figure
155    from matplotlib import cm
156
157    outputsDir = os.path.abspath(pkg_resources.resource_filename('geoplot', '../../outputs'))
158    testDataDir = pkg_resources.resource_filename('geoplot', 'tests/data')
159
160    import log_util
161    log_util.setupGeoplotConsoleHandler(log)
162
163    f = cdms.open(testDataDir + '/ammr_example2.nc')
164    var = f('temp', squeeze = 1) # remove the additional axis
165
166    builder = GridBuilderLatLon(var)
167
168    xLimits=(-13.0, 6.2) ; yLimits=(47.0, 61.0)
169    xLimits=(None, None) ; yLimits=(None, 60.0)
170
171    grid = builder.buildGrid(xLimits, yLimits)
172
173    xLimits=(-20.0, 15.0) ; yLimits=(42.0, 65.0)
174
175
176    figsize=(600 / 80, 800 / 80)
177    fig = Figure(figsize=figsize, dpi=80)
178    axes = fig.add_axes([0.1, 0.1, 0.8, 0.8])
179
180
181    bm = basemap.Basemap(llcrnrlon=xLimits[0],
182                         llcrnrlat=yLimits[0],
183                         urcrnrlon=xLimits[1],
184                         urcrnrlat=yLimits[1],
185                         resolution='h',
186                         suppress_ticks=False)
187    cmap = cm.jet
188    cmap.set_bad("w")
189
190    sm = bm.pcolormesh(grid.boundsX,grid.boundsY,grid.values,
191                   cmap=cmap, ax=axes, edgecolors='k')
192
193    bm.drawcoastlines(ax = axes)
194    bm.drawrivers(ax = axes, color='b')
195
196    canvas = FigureCanvasAgg(fig)
197    filename = "out.png"
198    canvas.print_figure(filename)
199    log.debug("Wrote " + filename)
Note: See TracBrowser for help on using the repository browser.