source: qesdi/geoplot/trunk/lib/geoplot/tests/basemap_test.py @ 5403

Subversion URL: http://proj.badc.rl.ac.uk/svn/ndg/qesdi/geoplot/trunk/lib/geoplot/tests/basemap_test.py@5403
Revision 5403, 6.3 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 
1from matplotlib.backends.backend_agg import FigureCanvasAgg
2from matplotlib.backends.backend_ps import FigureCanvasPS
3from matplotlib.collections import LineCollection
4from matplotlib.figure import Figure     
5from geoplot.mpl_imports import basemap
6
7import pkg_resources, os
8import time
9
10# National Grid parameters
11trueOrigin = (358.0, 49.0)
12falseOrigin = (-400000, 100000)
13
14BASEMAP_RESOLUTION = 'h'
15
16
17maxXLimits=(-20.0, 20.0);maxYLimits=(40.0,65.0)
18
19llcrnrlon = maxXLimits[0]
20llcrnrlat = maxYLimits[0]
21urcrnrlon = maxXLimits[1]
22urcrnrlat = maxYLimits[1]
23
24# Create a basemap with the origin at the true origin.
25nationalGrid = basemap.Basemap(projection='tmerc',
26                               lon_0=trueOrigin[0],
27                               lat_0=trueOrigin[1],
28                               llcrnrlon=llcrnrlon, llcrnrlat=llcrnrlat,
29                               urcrnrlon=urcrnrlon, urcrnrlat=urcrnrlat,
30                               resolution=BASEMAP_RESOLUTION
31                               )
32 
33
34def plotBasemap(xlimits, ylimits, filename):
35    outputsDir = os.path.abspath(pkg_resources.resource_filename('geoplot', '../../outputs'))   
36
37    #setup the figure
38    figsize=(600 / 80, 800 / 80)
39    fig = Figure(figsize=figsize, dpi=80)
40    axes = fig.add_axes([0.1, 0.1, 0.8, 0.8])
41       
42    st = time.time()
43    bm = basemap.Basemap(projection='tmerc',
44                               lon_0=trueOrigin[0],
45                               lat_0=trueOrigin[1],
46                               llcrnrlon=xlimits[0], llcrnrlat=ylimits[0],
47                               urcrnrlon=xlimits[1], urcrnrlat=ylimits[1],
48                               resolution=BASEMAP_RESOLUTION
49                               )
50   
51#    bm.drawcoastlines(ax = axes, linewidth=0.3, color="black")#
52#    bm.drawrivers(ax = axes, color='b', linewidth=0.3)
53
54#    eg_shape_file = '/disks/westerly1/pnorton/new_shapefiles/Admin_2D_SimplifyPolygon1'
55#    eg_shape_file = '/disks/westerly1/pnorton/new_shapefiles/river_2D_SimplifyPolygon1'
56#    eg_shape_file = '/disks/westerly1/pnorton/new_shapefiles/coast_outline_SimplifyPolygo'           
57#    eg_shape_file ="/disks/westerly1/pnorton/old_shapefiles/ukcp09_region_latlon"
58#    eg_shape_file ="/disks/westerly1/pnorton/old_shapefiles/ukcp09_coast_latlon"
59#    eg_shape_file ='/disks/westerly1/pnorton/old_shapefiles/ukcp09_river_latlon'
60
61#    eg_shape_file ="/disks/archive/real/spatial_files/region_latlon/ukcp09_region_latlon"
62#    eg_shape_file ="/disks/archive/real/spatial_files/coast_latlon/ukcp09_coast_latlon"
63#    eg_shape_file ='/disks/archive/real/spatial_files/river_latlon/ukcp09_river_latlon'
64
65#    eg_shape_file ="/disks/westerly1/pnorton/shapefiles_3/2DAdmin"
66#    eg_shape_file ="/disks/westerly1/pnorton/shapefiles_3/2DRiver"
67
68    eg_shape_file ="/disks/westerly1/pnorton/sf/shapefiles_4/2DAdmin"
69#    eg_shape_file ="/disks/westerly1/pnorton/sf/shapefiles_4/2DRiver"
70#   eg_shape_file ="/disks/westerly1/pnorton/sf/shapefiles_4/2DCoastb"
71
72    print "eg_shape_file =" , eg_shape_file
73    st = time.time()
74    bm.readshapefile(eg_shape_file, "region_outlines", color='black', 
75                     ax=axes, linewidth=0.3)
76   
77    print "read shape file in %.2fs" % (time.time() - st,)
78   
79#    print bm.region_outlines
80   
81    print "drawn in %fs" % (time.time() - st,)
82   
83#    import geoplot.tests.readTest
84#    st = time.time()
85#    segs = geoplot.tests.readTest.read()
86#    print "read segments in %.2fs" % (time.time() - st,)
87
88    xlims = axes.get_xlim()
89    ylims = axes.get_ylim()
90   
91    xmid = (xlims[0] + xlims[1])/2.0
92    ymid = (ylims[0] + ylims[1])/2.0
93   
94
95    segs = bm.region_outlines
96   
97    print "plot width = %s" % (xlims[1] - xlims[0], )
98   
99    print "total points in segments = ", sum([len(x) for x in segs])
100       
101    visibleSegments = []
102   
103    for linePoints in segs:
104       
105        visiblePoints = []
106       
107        for (x, y) in linePoints:
108            if isPointInLimits(x, y, xlims[0], ylims[0], xlims[1], ylims[1]):
109                visiblePoints.append( (x,y) )
110       
111        if len(visiblePoints) > 1:
112            visibleSegments.append(visiblePoints)
113               
114           
115    segs = visibleSegments
116   
117    print "total points in cropped segments = ", sum([len(x) for x in segs])
118   
119   
120    segs2 = []
121   
122    useEvery = 4
123   
124    for x in segs:
125        l = []
126       
127        l.append(x[0])
128       
129        i = 0
130        while i < len(x):
131            l.append(x[i])           
132            i += useEvery
133       
134        if len(x) < useEvery:
135            mid = int(float(len(x))/2.0)
136            l.append(x[mid])           
137            l.append(x[-1])
138       
139        segs2.append(l)
140   
141    print "total points in reduced segments = ", sum([len(x) for x in segs2])   
142   
143   
144#    for i in range(10000):
145#        s = i * 1000
146#        segments.append (
147#            [(xlims[0], ymid + s),(xmid, ymid),(xlims[1], ymid - s)]
148#            )
149    axes.collections = []
150#    col = LineCollection(segs); print "using cropped segments"
151    col = LineCollection(segs2) ; print "using reduced segments"
152    col.set_color('black')
153    col.set_linewidth(0.3)
154    axes.add_collection(col)
155   
156#    (xMin, yMin) = bm(xlimits[0], ylimits[0])
157#    (xMax, yMax) = bm(xlimits[1], ylimits[1])
158##
159#    axes.set_xlim(float(xMin), float(xMax))
160#    axes.set_ylim(float(yMin), float(yMax))   
161#   
162    print "axes.get_xlim() =" , axes.get_xlim()
163    print "axes.get_ylim() =" , axes.get_ylim()
164       
165   
166    froot, ext = os.path.splitext(filename)
167
168    if ext == '.ps':
169        CanvasClass=FigureCanvasPS
170    elif ext == '.png':
171        CanvasClass=FigureCanvasAgg
172    else:
173        raise Exception("Unknown format to write to , " + ext)
174
175    canvas = CanvasClass(fig)   
176   
177   
178    fullPath = os.path.join(outputsDir, filename)
179    canvas.print_figure(fullPath)
180    print "Wrote " + fullPath
181
182def isPointInLimits(x, y, xmin, ymin, xmax, ymax):
183    if x >= xmin and x <= xmax and y >=ymin and y <= ymax:
184        return True
185    else:
186        return False
187
188
189if __name__ == '__main__':
190   
191   
192    xlimits=(-10.0, 19.0);ylimits=(46.0,63.0)
193     
194    xlimits=(-4.0, 2.0) ; ylimits=(50.0, 54.0)
195#    xlimits=(-10.0, 4.0) ; ylimits=(47.5, 62.5)
196
197#    xlimits = (-10.0, 19.0); ylimits = (46.0, 63.0)
198   
199#    xlimits=(-10.0, 6.0) ; ylimits=(47.5, 62.5)
200
201       
202    plotBasemap(xlimits, ylimits, 'basemap_test.ps')
Note: See TracBrowser for help on using the repository browser.