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

Subversion URL: http://proj.badc.rl.ac.uk/svn/ndg/qesdi/geoplot/trunk/lib/geoplot/grid_builder_lat_lon_subset.py@5890
Revision 5890, 5.6 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
12import time
13
14#third party modules
15import numpy as N
16
17#internal modules
18from geoplot.grid_builder_lat_lon import GridBuilderLatLon
19from geoplot.grid_builder_base import GridBuilderBase
20from geoplot.grid import Grid
21
22from geoplot.range import Range
23import geoplot.utils
24import math
25
26#set the log
27log = logging.getLogger(__name__)
28
29class GridBuilderLatLonSubset(GridBuilderLatLon):
30    """
31    A grid builder that knows how to extract a lat-long grid from a cdms file.
32    """
33
34    def __init__(self, cdmsVar):
35        """
36        Constructs a GridBuilderLatLon object
37        """
38        GridBuilderLatLon.__init__(self, cdmsVar)
39
40    def buildGrid(self, targetWidth, targetHeight, xLimits=None, yLimits=None):
41        """
42        builds a grid object from the data in the cdmsVar ascociated with the
43        grid builder object.
44       
45        If xLimits or yLimits values are given then the resulting grid will
46        contain the portion of the data in the cdms variable that falls within
47        the limits given.
48       
49        @keyword xLimits: (optional) longitude limits of the resulting grid
50        @type xLimits: a tuple of (MinLongitude, MaxLongitude)
51        @keyword yLimits: (optional) latitude limits of the resulting grid
52        @type yLimits: a tuple of (MinLatitude, MaxLatitude)
53        @return: a grid built using the cdms variable data
54        @rtype: geoplot.Grid
55        """
56#        targetWidth = 12
57#        targetHeight = 7
58        reducedVar = self._getResizedVar(xLimits, yLimits)
59       
60        factor = 3.0
61        targetShape = (targetHeight* factor, targetWidth * factor)
62       
63        log.debug("xLimits = %s, yLimits = %s" % (xLimits, yLimits,))
64       
65        reducedLonMidpoints = reducedVar.getLongitude()[:]
66        reducedLatMidpoints = reducedVar.getLatitude()[:]
67       
68        subsetVar = self._subsetVariable(reducedVar, targetShape)
69       
70        # the subsetVar now contains fewer points than the reducedVar but we
71        # still need the bounds and the midpoints to fill the original area.
72       
73        # Assuming that the grid is regular (which it must be otherwise the
74        # sub-setting would drastically distort the results) we can just calculate
75        # the bounds + midpoints for the new number of points.
76       
77        if GridBuilderLatLon._areAxisInOrderYX(subsetVar):
78            lonbounds = self._buildReducedBounds(reducedLonMidpoints, subsetVar.shape[1] + 1)
79            latbounds = self._buildReducedBounds(reducedLatMidpoints, subsetVar.shape[0] + 1)
80        else:
81            lonbounds = self._buildReducedBounds(reducedLonMidpoints, subsetVar.shape[0] + 1)
82            latbounds = self._buildReducedBounds(reducedLatMidpoints, subsetVar.shape[1] + 1)
83               
84        lonmidpoints = geoplot.utils.getMidpoints(lonbounds)
85        latmidpoints = geoplot.utils.getMidpoints(latbounds)
86       
87        gridBoundsX, gridBoundsY     = N.meshgrid(lonbounds, latbounds,)
88        gridMidpointX, gridMidpointY = N.meshgrid(lonmidpoints, latmidpoints)
89               
90        gridValues = self._buildGridValues(subsetVar)
91       
92#        st = time.time()
93#        uniqueValues = N.unique(self.cdmsVar.getValue())
94#        log.debug("got unique values in = %s" % (time.time() - st,))
95       
96        return Grid(gridBoundsX, gridBoundsY, gridMidpointX, gridMidpointY, gridValues,
97                    self.cdmsVar[:].max(), self.cdmsVar[:].min())
98   
99    def _subsetVariable(self, var, targetShape):
100       
101        ratio = [int(math.floor(var.shape[x] / float(targetShape[x]))) for x in range(2)]
102       
103        for i in range(2):
104            if ratio[i] == 0:
105                ratio[i] = 1
106       
107       
108        if GridBuilderLatLon._areAxisInOrderYX(var):
109            subsetVar = var[::ratio[0], ::ratio[1]]
110        else:
111            subsetVar = var[::ratio[1], ::ratio[0]]
112       
113        #subsetVar = var
114        log.debug("original shape = %s, targetShape = %s, ratio = %s, subsetVar.shape = %s " % \
115                  (var.shape, targetShape, ratio, subsetVar.shape))
116       
117        return subsetVar
118   
119    def _buildReducedBounds(self, originalMidpoints, size):
120        """
121        Builds a bounds array with a reduced size from a set of original
122        midpoints.
123        """
124       
125        originalBoundsRange = self._calculateBoundsRangeFromMidpoints(originalMidpoints)
126       
127        # need the bounds to be in the same order as the data so check the
128        # original lat/lon values to get the order     
129        valsLowToHigh = originalMidpoints[0] < originalMidpoints[-1]
130       
131        if valsLowToHigh:
132            bounds = N.linspace(originalBoundsRange.minimum, originalBoundsRange.maximum, size)
133        else:
134            bounds = N.linspace(originalBoundsRange.maximum, originalBoundsRange.minimum, size)
135         
136        return bounds
137   
138    def _calculateBoundsRangeFromMidpoints(self, midpoints):
139       
140        valsLowToHigh = midpoints[0] < midpoints[-1]
141       
142        # assuming the midpoints are equally spaced and the
143        # bounds are mid-way between the midpoints
144        if valsLowToHigh:
145            spacing = midpoints[1] - midpoints[0]
146            rangeMin = midpoints[0] - spacing/2.0
147            rangeMax = midpoints[-1] + spacing/2.0
148        else:
149            spacing = midpoints[0] - midpoints[1]           
150            rangeMin = midpoints[-1] - spacing/2.0
151            rangeMax = midpoints[0] + spacing/2.0
152           
153        return Range(rangeMin, rangeMax)
Note: See TracBrowser for help on using the repository browser.