source: qesdi/geoplot/trunk/lib/geoplot/filtered_basemap.py @ 6103

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

Trying to get the different colour bars and plot options to work with the different legend types.

Line 
1'''
2Created on 18 Jun 2009
3
4@author: pnorton
5
6The filtered basemap inherits from Basemap and changes it in 3 main ways:
7
8 * The drawcoastlines and drawrivers methods now accept xLimits and yLimits and
9   only the coast or river lines within these limits should be drawn on the map.
10   This is important as it allows a single basemap to be used in multiple plots.
11   
12 * The self.coastsegs and self.riversegs arrays are parsed into Line objects, these
13   are cached after the first call. The Line objects makes calculating which
14   lines to draw faster. 
15   
16 * Implements a fix for the CAIRO backend that splits up lines which are too long
17   into mutliple shorter lines.
18'''
19
20import time
21import logging
22
23from mpl_toolkits.basemap import Basemap
24
25from matplotlib.collections import LineCollection
26import matplotlib.pyplot as plt
27
28from geoplot.utils import ispointInLimits, isRangeInLimits
29
30log = logging.getLogger(__name__)
31
32class Line(object):
33    """
34    A simple object that stores a number of x,y points along with a latitude and longitude range.
35    The latitude and longitude range correspond to the total range of the line.
36    (i.e the maximum and minimum of all the points in the line)
37    """
38   
39    def __init__(self, points, lonRange, latRange):
40        self.points = points
41        self.lonRange = lonRange
42        self.latRange = latRange
43
44class FilteredBasemap(Basemap):
45    """
46    A Filtered version of the Basemap object that should only draw coast lines
47    and river lines within the ranges specified.
48    """
49   
50
51    def __init__(self, *args, **kwargs):
52        self.__lines = None
53        self.__riverLines = None
54        Basemap.__init__(self, *args, **kwargs)
55
56    def drawcoastlines(self,linewidth=1.,color='k',antialiased=1,ax=None,zorder=None,
57                       xLimits=None, yLimits=None):
58       
59        if self.resolution is None:
60            raise AttributeError, 'there are no boundary datasets associated with this Basemap instance'
61       
62        # get current axes instance (if none specified).
63        if ax is None and self.ax is None:
64            try:
65                ax = plt.gca()
66            except:
67                import matplotlib.pyplot as plt
68                ax = plt.gca()
69               
70        elif ax is None and self.ax is not None:
71            ax = self.ax
72       
73        #log.debug("self.coastsegs has %s lines and %s points" % (len(self.coastsegs), sum([len(x) for x in self.coastsegs])))
74       
75        xLimitsMapUnits, yLimitsMapUnits = self.__call__(xLimits, yLimits)
76       
77        segs = self.__getSegments(xLimitsMapUnits, yLimitsMapUnits)
78
79        #log.debug("segs has %s lines and %s points" % (len(segs),  sum([len(x) for x in segs])))
80       
81        self._cairoLengthFix(segs)
82
83        coastlines = LineCollection(segs, antialiaseds=(antialiased,))
84       
85        coastlines.set_color(color)
86        coastlines.set_linewidth(linewidth)
87        coastlines.set_label('_nolabel_')
88       
89        if zorder is not None:
90            coastlines.set_zorder(zorder)
91           
92        ax.add_collection(coastlines)
93        # set axes limits to fit map region.
94        self.set_axes_limits(ax=ax)
95        return coastlines
96
97
98    def drawrivers(self,linewidth=0.5,color='k',antialiased=1,ax=None,zorder=None,
99                   xLimits=None, yLimits=None):
100
101        st = time.time()
102       
103        if self.resolution is None:
104            raise AttributeError, 'there are no boundary datasets associated with this Basemap instance'
105       
106        # read in river line segments, only keeping those that
107        # intersect map boundary polygon.
108        if not hasattr(self,'riversegs'):
109            self.riversegs, types = self._readboundarydata('rivers')
110           
111        # get current axes instance (if none specified).
112        if ax is None and self.ax is None:
113            try:
114                ax = plt.gca()
115            except:
116                import matplotlib.pyplot as plt
117                ax = plt.gca()
118        elif ax is None and self.ax is not None:
119            ax = self.ax
120       
121        xLimitsMapUnits, yLimitsMapUnits = self.__call__(xLimits, yLimits)
122       
123        segs = self.__getRiverSegments(xLimitsMapUnits, yLimitsMapUnits)
124
125        #log.debug("segs has %s lines and %s points" % (len(segs), sum([len(x) for x in segs])))
126 
127        self._cairoLengthFix(segs)
128 
129        rivers = LineCollection(segs, antialiaseds=(antialiased,))
130        rivers.set_color(color)
131        rivers.set_linewidth(linewidth)
132        rivers.set_label('_nolabel_')
133       
134        if zorder is not None:
135            rivers.set_zorder(zorder)
136           
137        ax.add_collection(rivers)
138        # set axes limits to fit map region.
139        self.set_axes_limits(ax=ax)
140       
141        log.debug("added rivers linecollection in %ss" % (time.time() - st,))
142        return rivers
143   
144   
145    def __getSegments(self, xLimits, yLimits):
146        """
147        Returns the coast segments within the limits xLimits and yLimits.
148        xLimits and yLimits are given in map units. If xLimits and yLimits are both
149        None then all of self.coastsegs is returned.
150        """
151       
152        if xLimits == None and yLimits == None:
153            return self.coastsegs
154       
155        #if the lines cache has't been build yet then build it using the coastsegs
156        if self.__lines == None:
157            self.__lines = self.__buildLines(self.coastsegs)
158       
159        #filter the lines using the limits
160        segs = self.__filterLines(xLimits, yLimits, self.__lines) 
161
162        return segs
163   
164    def __getRiverSegments(self, xLimits, yLimits):
165        """
166        Returns the river segments within the limits xLimits and yLimits.
167        xLimits and yLimits are given in map units. If xLimits and yLimits are both
168        None then all of self.riversegs is returned.
169        """
170               
171        if xLimits == None and yLimits == None:
172            return self.riversegs
173       
174        if self.__riverLines == None:
175            self.__riverLines = self.__buildLines(self.riversegs)
176       
177        segs = self.__filterLines(xLimits, yLimits, self.__riverLines)
178       
179        return segs
180   
181    def __buildLines(self, segs):
182        """
183        Builds a list of Line object for the segments given
184        """
185       
186        lines = []
187       
188        for points in segs:
189           
190            lons = [i[0] for i in points]
191            lats = [i[1] for i in points]
192           
193            lonRange = (min(lons), max(lons))
194            latRange = (min(lats), max(lats))
195           
196            lines.append(Line(points, lonRange, latRange))
197       
198        return lines
199
200    def __filterLines(self, xLimits, yLimits, lines):
201        """
202        Work out which of the lines fits between the limits of the plot
203        """
204       
205        def isLineInSelectionLimits(line):
206           
207            if not isRangeInLimits(line.lonRange, xLimits):
208                return False
209           
210            if not isRangeInLimits(line.latRange, yLimits):
211                return False
212           
213            return True
214       
215        # filter out the lines which are competly outside the limits
216        linesInLimits = filter(isLineInSelectionLimits, lines)
217       
218        segs = [l.points for l in linesInLimits]
219       
220        # now trim down the segments to only include bits of lines which are
221        # inside the limits
222        segs = self.__filterSegments(segs, xLimits, yLimits)
223       
224        return segs           
225       
226    def __filterSegments(self, segments, xLimits, yLimits):
227        """
228        Attempt to trim the segments in each line, removing any segments
229        at the start and end of the line that aren't within the limit.
230       
231        Not removing points in the middle of the line to avoid joining lines
232        cutting acros the image.
233       
234        Could split the lines in to multiple bits tho...
235        """
236       
237        segs = []
238        for line in segments:
239                       
240            #trim the start of the line
241            startIndex = 0
242            for i in range(len(line)-1):
243                if ispointInLimits(line[i], xLimits, yLimits):
244                    break
245                else:
246                    startIndex = i
247           
248            if startIndex == len(line) -1:
249                continue
250           
251            #trim the end of the line
252            endIndex = len(line)-1
253            for i in range(len(line)-1, startIndex, -1):
254                if ispointInLimits(line[i], xLimits, yLimits):
255                    break
256                else:
257                    endIndex = i
258           
259            #only one point of the line visible
260            if endIndex == startIndex:
261                continue
262               
263           
264            segs.append(line[startIndex:endIndex+1])
265       
266       
267        return segs
268   
269    def _cairoLengthFix(self, segs):
270        """
271        There is a limit on the length of lines that can be handled by the
272        cairo renderer, some of the segments lines are longer than this so
273        they need to be split up if they are to work with the cairo renderer.
274        """
275       
276        #!TODO: tidy this code up
277       
278        newLines = []
279        removeLines = []
280       
281        for i, x in enumerate(segs):
282
283            if len(x) > 18000:
284                log.debug("found i =%s len(x) = %s" % (i, len(x)))
285                removeLines.append(i)
286               
287                step = 18000
288                for i in range(100):
289                    f = step*i
290                    to = step*(i+1) -1
291                    if to > len(x):
292                        to = len(x)
293
294                    # can't draw a line with one point, make sure there
295                    # are at least 2 included
296                    if f - to == 1:
297                        f -= 1
298                       
299                    newLines.append(x[f:to])
300                   
301                    if to == len(x):
302                        break
303       
304        #remove the long lines
305        removeLines.reverse()
306        for i in removeLines:
307            segs.pop(i)
308       
309        #add the new segments back in
310        [segs.append(x) for x in newLines]
311       
Note: See TracBrowser for help on using the repository browser.