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

Subversion URL: http://proj.badc.rl.ac.uk/svn/ndg/qesdi/geoplot/trunk/lib/geoplot/tests/shapefile_draw_test.py@5403
Revision 5403, 4.5 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
6import geoplot.tests.readTest as readSegments
7import pkg_resources, os
8import time
9
10# National Grid parameters
11trueOrigin = (358.0, 49.0)
12falseOrigin = (-400000, 100000)
13
14BASEMAP_RESOLUTION = 'c'
15
16
17def makeLLBasemap(xlimits, ylimits):
18    bm = basemap.Basemap(lon_0=0.0,
19                           lat_0=0.0,
20                           llcrnrlon=xlimits[0],
21                           llcrnrlat=ylimits[0],
22                           urcrnrlon=xlimits[1],
23                           urcrnrlat=ylimits[1],
24                           resolution=BASEMAP_RESOLUTION,
25                           suppress_ticks=False)
26    return bm
27
28def makeNGBasemap(xlimits, ylimits):
29    bm = basemap.Basemap(projection='tmerc',
30                           lon_0=trueOrigin[0],
31                           lat_0=trueOrigin[1],
32                           llcrnrlon=xlimits[0],
33                           llcrnrlat=ylimits[0],
34                           urcrnrlon=xlimits[1],
35                           urcrnrlat=ylimits[1],
36                           resolution=BASEMAP_RESOLUTION,
37                           suppress_ticks=False)
38    return bm   
39
40
41def plotBasemap(xlimits, ylimits, filename):
42    startTime = time.time()
43    outputsDir = os.path.abspath(pkg_resources.resource_filename('geoplot', '../../outputs'))   
44   
45    st = time.time()
46    #setup the figure
47    figsize=(600 / 80, 800 / 80)
48    fig = Figure(figsize=figsize, dpi=80)
49    axes = fig.add_axes([0.1, 0.1, 0.8, 0.8])
50       
51    eg_shape_file ="/disks/archive/real/spatial_files/region_latlon/ukcp09_region_latlon"
52#    eg_shape_file ="/disks/archive/real/spatial_files/coast_latlon/ukcp09_coast_latlon"
53#    eg_shape_file ='/disks/archive/real/spatial_files/river_latlon/ukcp09_river_latlon'
54
55    print "setup in %.3fs" % (time.time() - st,)
56
57    st = time.time()
58#    bm = makeLLBasemap(xlimits, ylimits)
59    bm = makeNGBasemap(xlimits, ylimits)
60    print "made basemap in %.3fs" % (time.time() - st,)
61       
62#    ngBm = makeNGBasemap(xlimits, ylimits)
63
64    st = time.time()
65    bm.drawrivers(ax=axes, color='b')
66    print "drawrivers in %.3fs" % (time.time() - st,)
67   
68    st = time.time()
69    segs = readSegmentsFromShapefile(eg_shape_file, bm)
70    print "read segments form shapefile in %.3fs" % (time.time() - st,)
71   
72    st = time.time()
73    writeSegments(segs,'segments.txt')
74    print "wrote segments in %.3fs" % (time.time() - st,)
75   
76    st = time.time()
77    segsFromFile = readSegments.read('segments.txt')
78    print "read segments from file in %.3fs" % (time.time() - st,)
79   
80    st = time.time()
81    drawSegments(axes, segsFromFile)
82    print "drawn segments in %.3fs" % (time.time() - st,)
83   
84    froot, ext = os.path.splitext(filename)
85
86    if ext == '.ps':
87        CanvasClass=FigureCanvasPS
88    elif ext == '.png':
89        CanvasClass=FigureCanvasAgg
90    else:
91        raise Exception("Unknown format to write to , " + ext)
92
93    canvas = CanvasClass(fig)   
94   
95   
96    fullPath = os.path.join(outputsDir, filename)
97   
98    st = time.time()
99    canvas.print_figure(fullPath)
100    print "wrote file in %.3fs" % (time.time() - st,)
101   
102    print "Wrote " + fullPath + " in %.3fs" % (time.time() - startTime,)
103
104
105def readSegmentsFromShapefile(shapefilePath, bm):
106    bm.readshapefile(shapefilePath, 'region_outlines', drawbounds=False)
107    return bm.region_outlines
108   
109def drawSegments(axes, segs):
110    col = LineCollection(segs)
111    col.set_color('black')
112    col.set_linewidth(0.3)
113    axes.add_collection(col)
114   
115def isPointInLimits(x, y, xmin, ymin, xmax, ymax):
116    if x >= xmin and x <= xmax and y >=ymin and y <= ymax:
117        return True
118    else:
119        return False
120
121
122def writeSegments(segments, outputFile):
123    fout = open(outputFile, 'w')
124    lines = []
125    print "len(segments) =" , len(segments) 
126    for line in segments:
127        linePoints = ""
128        for (x,y) in line:
129            linePoints += "%s,%s:" % (x,y)
130        linePoints = linePoints[:-1] + "\n"
131        lines.append(linePoints)
132       
133    fout.writelines(lines)
134    fout.close()
135if __name__ == '__main__':
136   
137   
138#    xlimits=(-12.0, 19.0) ; ylimits=(46.0,63.0)
139   
140#    xlimits=(-12.0, 4.0)  ; ylimits=(47.5, 62.5)
141   
142    xlimits=(-4.0, 2.0)   ; ylimits=(50.0, 54.0)
143       
144    plotBasemap(xlimits, ylimits, 'll_basemap_test.ps')
Note: See TracBrowser for help on using the repository browser.