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
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
|
'''
Plot map data in different projections supported by proj.4
Can be called either as
> python plot.py definitions.txt [out]
or
> python plotproj.py "+proj=lcc +lat_1=20" [out]
The output dir is optional. Uses current dir if not set.
Change PROJ and COASTLINE_DATA to path on your local system
before running the script.
'''
import os
import sys
import json
import subprocess
import numpy as np
import matplotlib.pyplot as plt
PROJ = 'proj'
COASTLINE_DATA = 'data/coastline.geojson'
GRATICULE_WIDTH = 15
N_POINTS = 100
LAT_MIN = -90
LAT_MAX = 90
LON_MIN = -180
LON_MAX = 180
def project(coordinates, proj_string, in_radians=False):
'''
Project geographical coordinates
Input:
------
coordinates: numpy ndarray of size (N,2) and type double.
longitude, latitude
proj_string: Definition of output projection
Out:
----
numpy ndarray with shape (N,2) with N pairs of 2D
coordinates (x, y).
'''
# proj expects binary input to be in radians
if not in_radians:
coordinates = np.deg2rad(coordinates)
# set up cmd call. -b for binary in/out
args = [PROJ, '-b']
args.extend(proj_string.split(' '))
proc = subprocess.Popen(args, stdin=subprocess.PIPE, stdout=subprocess.PIPE)
stdout, _ = proc.communicate(coordinates.tobytes())
out = np.frombuffer(stdout, dtype=np.double)
return np.reshape(out, coordinates.shape)
def meridian(lon, lat_min, lat_max):
'''
Calculate meridian coordinates.
'''
coords = np.zeros((N_POINTS, 2))
coords[:, 0] = lon
coords[:, 1] = np.linspace(lat_min, lat_max, N_POINTS)
return coords
def parallel(lat, lon_min, lon_max):
'''
Calculate parallel coordinates.
'''
coords = np.zeros((N_POINTS, 2))
coords[:, 0] = np.linspace(lon_min, lon_max, N_POINTS)
coords[:, 1] = lat
return coords
def plotproj(proj_str, coastline, frame, graticule, outdir):
'''
Plot map.
'''
fig = plt.figure()
axes = fig.add_subplot(111)
# Plot frame
for line in frame:
line = project(line, proj_str)
x, y = zip(*line)
plt.plot(x, y, '-k')
# Plot coastline
for feature in coastline['features']:
C = np.array(feature['geometry']['coordinates'])
if np.any(C[:, 0] > 180.0):
C[C[:, 0] > 180.0, 0] = 180.0
if np.any(C[:, 0] < -180.0):
C[C[:, 0] < -180.0, 0] = -180.0
if np.any(C[:, 1] > 90.0):
C[C[:, 1] > 90.0, 1] = 90.0
if np.any(C[:, 1] < -90.0):
C[C[:, 1] < -90.0, 1] = -90.0
coords = project(C, proj_str)
x, y = zip(*coords)
plt.plot(x, y, '-k', linewidth=0.5)
# Plot graticule
graticule = project(graticule, proj_str)
for feature in graticule:
x, y = zip(*feature)
plt.plot(x, y, '-k', linewidth=0.2)
axes.axis('off')
font = {'family': 'serif',
'color': 'black',
'style': 'italic',
'size': 12}
plt.suptitle(proj_str, fontdict=font)
plt.autoscale(tight=True)
name = proj_str.split()[0].split('=')[1]
plt.savefig(outdir + '/' + name + '.png', dpi=300)
# Clean up
fig = None
del fig
axes = None
del axes
plt.close()
if __name__ == "__main__":
'''
'''
with open(COASTLINE_DATA) as data:
coastline = json.load(data)
if os.path.exists(sys.argv[1]):
# it must be a file with proj-strings
with open(sys.argv[1]) as projs:
projs = projs.readlines()
else:
#assumes it is a single roj string
projs = [sys.argv[1]]
outdir = '.'
if len(sys.argv) > 2:
outdir = sys.argv[2]
frame = []
frame.append(parallel(LAT_MIN, LON_MIN, LON_MAX))
frame.append(parallel(LAT_MAX, LON_MIN, LON_MAX))
# build graticule
graticule = []
for lon in range(LON_MIN, LON_MAX+1, GRATICULE_WIDTH):
graticule.append(meridian(lon, LAT_MIN, LAT_MAX))
for lat in range(LAT_MIN, LAT_MAX+1, GRATICULE_WIDTH):
graticule.append(parallel(lat, LON_MIN, LON_MAX))
for i, proj in enumerate(projs):
proj_str = proj.strip()
print(i, proj_str)
plotproj(proj_str, coastline, frame, graticule, outdir)
|
