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

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

Cleaned up a the import statements for the geoplot modules, also fixed a couple of the tests.

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