source: qesdi/geoplot/trunk/lib/geoplot/map_base.py @ 5403

Subversion URL: http://proj.badc.rl.ac.uk/svn/ndg/qesdi/geoplot/trunk/lib/geoplot/map_base.py@5403
Revision 5403, 8.9 KB checked in by pnorton, 12 years ago (diff)

Moved QESDI tree from DCIP repository
 http://proj.badc.rl.ac.uk/svn/dcip/qesdi@3900.

Line 
1"""
2map_base.py
3============
4
5A base map class that holds the common map functionality. A map object knows how
6to create a basemap instance and draw it onto a given axes.
7
8"""
9
10import logging
11
12from matplotlib.collections import LineCollection
13import numpy as N
14
15import time
16
17log = logging.getLogger(__name__)
18log.setLevel(logging.DEBUG)
19
20BASEMAP_RESOLUTION = 'i'
21
22
23class MapBase(object):
24    """
25    Represents a type of map that can be drawn onto a matplotlib.axes.
26    """
27
28    def __init__(self, xLimits, yLimits, drawCoast, drawRivers, addShapefile, **kwargs ):
29        """
30        Constructs the map object.
31
32        @param xLimits:a tuple representing the (longitude) limits of the map
33        @type xLimits: a tuple of 2 floats, corresponding to (min, max)
34        @param yLimits:a tuple representing the (latitude) limits of the map
35        @type yLimits: a tuple of 2 floats , corresponding to (min, max)
36        @param drawCoast: indicates if coaslines should be drawn
37        @type drawCoast: boolean
38        @param drawRivers: indicates if rivers should be drawn
39        @type drawRivers: boolean
40        @param addShapefile: location of shapefile to overlay
41        @type addShapefile: string
42        """
43
44        self.addShapefile = addShapefile
45       
46        self._drawCoast = drawCoast
47        self._drawRivers = drawRivers
48       
49        self._checkLimits(xLimits, "xLimits")
50        self._xLimits = xLimits
51        self._checkLimits(yLimits, "yLimits")
52        self._yLimits = yLimits
53
54        self._setBasemap()
55
56    def drawMap(self, axes):
57        """
58        Draws the map onto an axes object.
59
60        @param axes: the axes on which to draw the map and grid.
61        @type axes: an instance of the matplotlib.axes class
62        """
63
64        startTime = time.time()
65        log.info("Drawing map")
66
67        log.debug("limits:" + str(self.xLimits) + ", " + str(self.yLimits))
68
69        log.debug("Values of drawCoast, drawRivers, addShapefile are: " + str(self.drawCoast) + "," + str(self.drawRivers) + "," + str(self.addShapefile))
70       
71#        self.basemap.drawmeridians(N.arange(0.,420.,10.),
72#                                   ax = axes,
73#                                   labels=[0,0,0,1])
74#       
75#        self.basemap.drawparallels(N.arange(-90.,120.,10.),
76#                                   ax = axes,
77#                                   labels=[1,0,0,0])       
78       
79        if self.drawCoast:
80
81            st = time.time()
82            self.basemap.drawcoastlines(ax = axes, linewidth=0.5, color="black",
83                                        antialiased=0,
84                                        xLimits=self.xLimits, yLimits=self.yLimits)
85           
86            log.debug(" drawn coastlines in %.2fs" % (time.time() - st,))
87
88        if self.drawRivers:
89            self.basemap.drawrivers(ax = axes, color='b', linewidth=0.3)
90           
91        if self.addShapefile != False:
92            log.debug("Adding shapefile: " + str(self.addShapefile)) 
93            self._drawShapefile(axes)
94
95        self._configAxes(axes)
96
97        self._resizePlot(axes)
98       
99        log.debug(" drawn map in %.2fs" % (time.time() - startTime,))
100
101    def _setBasemap(self):
102        raise NotImplementedError
103
104    def _getBasemapResolution(self):
105        "Returns the resolution that should be used when creating a basemap"
106       
107        # if the coast + rivers arn't needed don't load the resolution data,
108        # this optimisation is usefull optimisation is basemap takes quite a
109        # while to load the coast+river data, and this is skipped if resolution = None
110       
111        if self.drawCoast == False and self.drawRivers == False:
112            resolution = None
113        else:
114            resolution = BASEMAP_RESOLUTION
115       
116        log.debug("Using basemap resolution %s" % (resolution,))
117       
118        return resolution
119   
120    def _drawShapefile(self, axes):
121        "Uses the self.basemap to draw a shapefile on the axis."
122       
123        name = 'region_outlines'
124        color = 'black'
125        linewidth = 0.3
126       
127        self.basemap.readshapefile(self.addShapefile, name, 
128                                   drawbounds=False)
129       
130        segments = getattr(self.basemap, name)
131       
132        reducedSegments = self._cropShapefileSegments(segments)
133       
134        #add the reduced segments as a new line collection
135        col = LineCollection(reducedSegments)
136        col.set_color(color)
137        col.set_linewidth(linewidth)
138        axes.add_collection(col)
139        self.basemap.set_axes_limits(ax=axes)
140   
141    def _getShapefileLimits(self):
142        """
143        Returns the lat/lon limits for drawing a shapefile, will not return
144        limits larger than SHAPEFILE_LON_LIMITS and SHAPEFILE_LAT_LIMITS but
145        if self.xLimits and self.yLimits are smaller than these limits a
146        reduced limit will be returned.
147        """
148       
149        return self.xLimits, self.yLimits
150       
151   
152    def _cropShapefileSegments(self, segments):
153        """
154        reduced the shapefile segments to those which have at least 2
155        points within the shapefile limits.
156        """
157        lonLimits, latLimits = self._getShapefileLimits()
158       
159        log.debug("shapefile lonLimits=%s , latLimits = %s" % (lonLimits, latLimits,))
160       
161        # get the limits in map projection units
162        xLimits, yLimits = self.basemap(lonLimits, latLimits)
163       
164        log.debug("shapefile xLimits=%s , yLimits = %s" % (xLimits, yLimits,))
165       
166        # build a new list of reduced segments
167        remainingSegments = []
168        for linePoints in segments:
169           
170            remainingPoints = []
171           
172            for (x, y) in linePoints:
173                if self._isPointInLimits(x, y, xLimits, yLimits):
174                    remainingPoints.append( (x,y) )
175           
176            if len(remainingPoints) > 1:
177                # To avoid the ends of a line being connected across the map add
178                # all the line points not just thoes visible.
179                remainingSegments.append(linePoints)
180       
181        log.debug("number of points remaining = %s" % (sum([len(x) for x in remainingSegments]),))
182        return remainingSegments
183       
184    def _isPointInLimits(self, x, y, xLimits, yLimits):
185        """
186        Checks if a point (x,y) are within the limits xLimits and yLimits, all
187        values should be given in map projection units.
188        """
189        if x >= xLimits[0] and x <= xLimits[1] and y >=yLimits[0] and y <= yLimits[1]:
190            return True
191        else:
192            return False
193
194    def _resizePlot(self, axes):
195        """
196        during the drawing process the axes limits may be changed (differently depending
197        on what is plotted) this method re-applies the limits of the map so that the expected
198        area is displayed in the finished map.
199
200        @param axes: the axes on which to draw the map and grid.
201        @type axes: an instance of the matplotlib.axes class
202        """
203
204        (xMin, yMin) = self.basemap(self.xLimits[0], self.yLimits[0])
205        (xMax, yMax) = self.basemap(self.xLimits[1], self.yLimits[1])
206
207        #log.debug("xMin, xMax:" + str(xMin) + "," + str(xMax))
208        #log.debug("yMin, yMax:" + str(yMin) + "," + str(yMax))
209
210        axes.set_xlim(float(xMin), float(xMax))
211        axes.set_ylim(float(yMin), float(yMax))
212       
213    def _checkLimits(self, limits, name):
214       
215       
216        if limits[0] > limits[1]:
217            raise Exception("Can't set value of %s for %s, limits should be from low to high"\
218                            % (limits, name))
219
220    #getters and setters for the xLimit and yLimit properties, these are needed
221    #to try and ensure that the basemap which is dependent on the limits is
222    #rebuilt whenever a limit changes.
223
224    def __set_xLimits(self,value): 
225        self._checkLimits(value, "xLimits")
226        self._xLimits = value
227        self._setBasemap()
228    def __get_xLimits(self): return self._xLimits
229    xLimits = property(__get_xLimits, __set_xLimits)   
230
231    def __set_yLimits(self,value): 
232        self._checkLimits(value, "yLimits")
233        self._yLimits = value
234        self._setBasemap()
235    def __get_yLimits(self): return self._yLimits
236    yLimits = property(__get_yLimits, __set_yLimits)       
237
238    def __set_drawRivers(self, value):
239        if value != self._drawRivers:
240            self._drawRivers = value
241            self._setBasemap()
242           
243    def __get_drawRivers(self): return self._drawRivers
244    drawRivers = property(__get_drawRivers, __set_drawRivers)
245   
246    def __set_drawCoast(self, value):
247        if value != self._drawCoast:
248            self._drawCoast = value
249            self._setBasemap()
250           
251    def __get_drawCoast(self): return self._drawCoast
252    drawCoast = property(__get_drawCoast, __set_drawCoast)   
253
254if __name__ == '__main__':
255    #test a simple basemap
256    import pylab
257    import basemap
258   
259    bm = basemap.Basemap()
260    bm.drawcoastlines()
261    pylab.savefig("out.png")
262   
Note: See TracBrowser for help on using the repository browser.