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

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

Modified the geoplot grid builder so that it automatically masks infinite values.

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
57        varMin, varMax = self._getVarMinAndMax()
58       
59        infValueFound = False
60       
61        if varMax == N.inf:
62            infValueFound = True
63            varMax =  self._maskInfValsInVar(self.cdmsVar).max()
64
65        reducedVar = self._getResizedVar(xLimits, yLimits)
66       
67        factor = 1.0
68        targetShape = (targetHeight* factor, targetWidth * factor)
69       
70        log.debug("xLimits = %s, yLimits = %s" % (xLimits, yLimits,))
71       
72        reducedLonMidpoints = reducedVar.getLongitude()[:]
73        reducedLatMidpoints = reducedVar.getLatitude()[:]
74       
75        subsetVar = self._subsetVariable(reducedVar, targetShape)
76       
77        # the subsetVar now contains fewer points than the reducedVar but we
78        # still need the bounds and the midpoints to fill the original area.
79       
80        # Assuming that the grid is regular (which it must be otherwise the
81        # sub-setting would drastically distort the results) we can just calculate
82        # the bounds + midpoints for the new number of points.
83       
84       
85        if GridBuilderLatLon._areAxisInOrderYX(subsetVar):
86            lonbounds = self._buildReducedBounds(reducedLonMidpoints, subsetVar.shape[1] + 1)
87            latbounds = self._buildReducedBounds(reducedLatMidpoints, subsetVar.shape[0] + 1)
88        else:
89            lonbounds = self._buildReducedBounds(reducedLonMidpoints, subsetVar.shape[0] + 1)
90            latbounds = self._buildReducedBounds(reducedLatMidpoints, subsetVar.shape[1] + 1)
91               
92        lonmidpoints = geoplot.utils.getMidpoints(lonbounds)
93        latmidpoints = geoplot.utils.getMidpoints(latbounds)
94       
95        gridBoundsX, gridBoundsY     = N.meshgrid(lonbounds, latbounds,)
96        gridMidpointX, gridMidpointY = N.meshgrid(lonmidpoints, latmidpoints)
97               
98        gridValues = self._buildGridValues(subsetVar)
99       
100        if infValueFound:
101            gridValues = self._maskInfValsInVar(gridValues)
102
103#        st = time.time()
104#        uniqueValues = N.unique(self.cdmsVar.getValue())
105#        log.debug("got unique values in = %s" % (time.time() - st,))
106       
107        return Grid(gridBoundsX, gridBoundsY, gridMidpointX, gridMidpointY, gridValues,
108                    varMax, varMin)
109   
110    def _subsetVariable(self, var, targetShape):
111       
112        ratio = [int(math.floor(var.shape[x] / float(targetShape[x]))) for x in range(2)]
113       
114        for i in range(2):
115            if ratio[i] == 0:
116                ratio[i] = 1
117       
118       
119        if GridBuilderLatLon._areAxisInOrderYX(var):
120            subsetVar = var[::ratio[0], ::ratio[1]]
121        else:
122            subsetVar = var[::ratio[1], ::ratio[0]]
123       
124        #subsetVar = var
125        log.debug("original shape = %s, targetShape = %s, ratio = %s, subsetVar.shape = %s " % \
126                  (var.shape, targetShape, ratio, subsetVar.shape))
127       
128        return subsetVar
129   
130    def _buildReducedBounds(self, originalMidpoints, size):
131        """
132        Builds a bounds array with a reduced size from a set of original
133        midpoints.
134        """
135       
136        originalBoundsRange = self._calculateBoundsRangeFromMidpoints(originalMidpoints)
137       
138        # need the bounds to be in the same order as the data so check the
139        # original lat/lon values to get the order     
140        valsLowToHigh = originalMidpoints[0] < originalMidpoints[-1]
141       
142        if valsLowToHigh:
143            bounds = N.linspace(originalBoundsRange.minimum, originalBoundsRange.maximum, size)
144        else:
145            bounds = N.linspace(originalBoundsRange.maximum, originalBoundsRange.minimum, size)
146         
147        return bounds
148   
149    def _calculateBoundsRangeFromMidpoints(self, midpoints):
150       
151        valsLowToHigh = midpoints[0] < midpoints[-1]
152       
153        # assuming the midpoints are equally spaced and the
154        # bounds are mid-way between the midpoints
155        if valsLowToHigh:
156            spacing = midpoints[1] - midpoints[0]
157            rangeMin = midpoints[0] - spacing/2.0
158            rangeMax = midpoints[-1] + spacing/2.0
159        else:
160            spacing = midpoints[0] - midpoints[1]           
161            rangeMin = midpoints[-1] - spacing/2.0
162            rangeMax = midpoints[0] + spacing/2.0
163           
164        return Range(rangeMin, rangeMax)
Note: See TracBrowser for help on using the repository browser.