forked from MPAS-Dev/geometric_features
-
Notifications
You must be signed in to change notification settings - Fork 1
/
plot_features.py
executable file
·140 lines (121 loc) · 5.45 KB
/
plot_features.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
#!/usr/bin/env python
"""
This script plots a file containing multiple features onto a basemap using
matplotlib's basemap.
It requires basemap: http://matplotlib.org/basemap/
The -r flag is used to pass in a features file that will be plotted, and the -o
flag can optionally be used to specify the name of the image that will be
generated.
Author: Phillip J. Wolfram
Date: 08/25/2015
"""
import numpy as np
import json
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap, addcyclic
def plot_base(maptype): #{{{
if maptype == 'ortho':
map = Basemap(projection='ortho', lat_0=45, lon_0=-100, resolution='l')
map.drawmeridians(np.arange(0,360,30))
map.drawparallels(np.arange(-90,90,30))
elif maptype == 'mill':
map = Basemap(llcrnrlon=0,llcrnrlat=-90,urcrnrlon=360,urcrnrlat=90,projection='mill')
map.drawparallels(np.arange(-80,81,20),labels=[1,1,0,0])
map.drawmeridians(np.arange(0,360,60),labels=[0,0,0,1])
elif maptype == 'mill2':
map = Basemap(llcrnrlon=-180,llcrnrlat=-90,urcrnrlon=180,urcrnrlat=90,projection='mill')
map.drawparallels(np.arange(-80,81,20),labels=[1,1,0,0])
map.drawmeridians(np.arange(0,360,60),labels=[0,0,0,1])
elif maptype == 'robin':
map = Basemap(projection='robin',lon_0=0,lat_0=0)
map.drawparallels(np.arange(-80,81,20),labels=[1,1,0,0])
map.drawmeridians(np.arange(0,360,60),labels=[0,0,0,1])
elif maptype == 'robin2':
map = Basemap(projection='robin',lon_0=180,lat_0=0)
map.drawparallels(np.arange(-80,81,20),labels=[1,1,0,0])
map.drawmeridians(np.arange(0,360,60),labels=[0,0,0,1])
elif maptype == 'hammer':
map = Basemap(projection='hammer',lon_0=180)
map.drawmeridians(np.arange(0,360,30))
map.drawparallels(np.arange(-90,90,30))
elif maptype == 'northpole':
map = Basemap(projection='ortho', lat_0=90, lon_0=-100, resolution='l')
map.drawmeridians(np.arange(0,360,30))
map.drawparallels(np.arange(-90,90,30))
elif maptype == 'southpole':
map = Basemap(projection='ortho', lat_0=-90, lon_0=-100, resolution='l')
map.drawmeridians(np.arange(0,360,30))
map.drawparallels(np.arange(-90,90,30))
elif maptype == 'atlantic':
map = Basemap(projection='ortho', lat_0=0, lon_0=0, resolution='l')
map.drawmeridians(np.arange(0,360,30))
map.drawparallels(np.arange(-90,90,30))
elif maptype == 'pacific':
map = Basemap(projection='ortho', lat_0=0, lon_0=-180, resolution='l')
map.drawmeridians(np.arange(0,360,30))
map.drawparallels(np.arange(-90,90,30))
elif maptype == 'europe':
map = Basemap(projection='ortho', lat_0=0, lon_0=-90, resolution='l')
map.drawmeridians(np.arange(0,360,30))
map.drawparallels(np.arange(-90,90,30))
elif maptype == 'northamerica':
map = Basemap(projection='ortho', lat_0=0, lon_0=90, resolution='l')
map.drawmeridians(np.arange(0,360,30))
map.drawparallels(np.arange(-90,90,30))
else:
raise NameError("Didn't select a valid maptype!")
map.drawcoastlines(linewidth=0.25)
map.drawcountries(linewidth=0.25)
map.fillcontinents(color='#e0e0e0', lake_color='white')
map.drawmapboundary(fill_color='white')
return map #}}}
def plot_features_file(featurefile, plotname):
# open up the database
with open(featurefile) as f:
featuredat = json.load(f)
colors = ['blue', 'green', 'red', 'cyan', 'magenta', 'yellow', 'black']
markers = ['o', 's', 'v', '^', '>', '<', '*', 'p', 'D', 'h']
fig = plt.figure(figsize=(16,12),dpi=100)
for anum, maptype in enumerate(['mill2','mill', 'northpole', 'southpole']):
ax = fig.add_subplot(2,2,1+anum)
feature_num = 0
for feature in featuredat['features']:
polytype = feature['geometry']['type']
coords = feature['geometry']['coordinates']
feature = feature['properties']['name']
color_num = feature_num % len(colors)
marker_num = feature_num % len(markers)
map = plot_base(maptype)
try:
if polytype == 'MultiPolygon' or polytype == 'MultiLineString':
for poly in coords:
points = np.asarray(poly)
map.plot(points[:,0], points[:,1], linewidth=2.0, color=colors[color_num], latlon=True)
elif polytype == 'Polygon':
points = np.asarray(coords)
map.plot(points[:,0], points[:,1], linewidth=2.0, color=colors[color_num], latlon=True)
elif polytype == 'LineString':
points = np.asarray(coords)
# due to bug in basemap http://stackoverflow.com/questions/31839047/valueerror-in-python-basemap/32252594#32252594
lons = [points[0,0],points[0,0],points[1,0],points[1,0]]
lats = [points[0,1],points[0,1],points[1,1],points[1,1]]
map.plot(lons, lats, linewidth=2.0, color=colors[color_num], latlon=True)
elif polytype == 'Point':
points = np.asarray(coords)
map.plot(points[0], points[1], markers[marker_num], markersize=20, color=colors[color_num], latlon=True)
else:
assert False, 'Geometry %s not known.'%(polytype)
except:
print 'Error plotting %s for map type %s'%(feature, maptype)
feature_num = feature_num + 1
print 'saving ' + plotname
plt.savefig(plotname)
if __name__ == "__main__":
import argparse
parser = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawTextHelpFormatter)
parser.add_argument("-f", "--features_file", dest="features_file", help="Feature file to plot", metavar="FILE", required=True)
parser.add_argument("-o", "--features_plot", dest="features_plotname", help="Feature plot filename", metavar="FILE")
args = parser.parse_args()
if not args.features_plotname:
args.features_plotname = args.features_file.strip('.*') + '_plot.png'
plot_features_file(args.features_file, args.features_plotname)