#
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) |
---|

Line | |
---|---|

1 | """ |

2 | grid_builder_lat_lon.py |

3 | ============ |

4 | |

5 | A GridBuilder object that can extract a lat-lon grid from a cdmsVariable. The |

6 | variable must have a latitude and a longitude axis and no level or time data. |

7 | |

8 | """ |

9 | |

10 | #python modules |

11 | import logging |

12 | import time |

13 | |

14 | #third party modules |

15 | import numpy as N |

16 | |

17 | #internal modules |

18 | from geoplot.grid_builder_lat_lon import GridBuilderLatLon |

19 | from geoplot.grid_builder_base import GridBuilderBase |

20 | from geoplot.grid import Grid |

21 | |

22 | from geoplot.range import Range |

23 | import geoplot.utils |

24 | import math |

25 | |

26 | #set the log |

27 | log = logging.getLogger(__name__) |

28 | |

29 | class 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.