1 | # Python tools to manage netCDF files. |
---|
2 | # L. Fita, CIMA. Mrch 2019 |
---|
3 | # More information at: http://www.xn--llusfb-5va.cat/python/PyNCplot |
---|
4 | # |
---|
5 | # pyNCplot and its component geometry_tools.py comes with ABSOLUTELY NO WARRANTY. |
---|
6 | # This work is licendes under a Creative Commons |
---|
7 | # Attribution-ShareAlike 4.0 International License (http://creativecommons.org/licenses/by-sa/4.0) |
---|
8 | # |
---|
9 | ## Script for geometry calculations and operations as well as definition of different |
---|
10 | ### standard objects and shapes |
---|
11 | |
---|
12 | import numpy as np |
---|
13 | import matplotlib as mpl |
---|
14 | from mpl_toolkits.mplot3d import Axes3D |
---|
15 | import matplotlib.pyplot as plt |
---|
16 | |
---|
17 | ####### Contents: |
---|
18 | # deg_deci: Function to pass from degrees [deg, minute, sec] to decimal angles [rad] |
---|
19 | # position_sphere: Function to tranform fom a point in lon, lat deg coordinates to |
---|
20 | # cartesian coordinates over an sphere |
---|
21 | # surface_sphere: Function to provide an sphere as matrix of x,y,z coordinates |
---|
22 | # spheric_line: Function to transform a series of locations in lon, lat coordinates |
---|
23 | # to x,y,z over an 3D spaceFunction to provide coordinates of a line on a 3D space |
---|
24 | |
---|
25 | ## Plotting |
---|
26 | # plot_sphere: Function to plot an sphere and determine which standard lines will be |
---|
27 | # also drawn |
---|
28 | |
---|
29 | def deg_deci(angle): |
---|
30 | """ Function to pass from degrees [deg, minute, sec] to decimal angles [rad] |
---|
31 | angle: list of [deg, minute, sec] to pass |
---|
32 | >>> deg_deci([41., 58., 34.]) |
---|
33 | 0.732621346072 |
---|
34 | """ |
---|
35 | fname = 'deg_deci' |
---|
36 | |
---|
37 | deg = np.abs(angle[0]) + np.abs(angle[1])/60. + np.abs(angle[2])/3600. |
---|
38 | |
---|
39 | if angle[0] < 0.: deg = -deg*np.pi/180. |
---|
40 | else: deg = deg*np.pi/180. |
---|
41 | |
---|
42 | return deg |
---|
43 | |
---|
44 | def position_sphere(radii, alpha, beta): |
---|
45 | """ Function to tranform fom a point in lon, lat deg coordinates to cartesian |
---|
46 | coordinates over an sphere |
---|
47 | radii: radii of the sphere |
---|
48 | alpha: longitude of the point |
---|
49 | beta: latitude of the point |
---|
50 | >>> position_sphere(10., 30., 45.) |
---|
51 | (0.81031678432964027, -5.1903473778327376, 8.5090352453411846 |
---|
52 | """ |
---|
53 | fname = 'position_sphere' |
---|
54 | |
---|
55 | xpt = radii*np.cos(beta)*np.cos(alpha) |
---|
56 | ypt = radii*np.cos(beta)*np.sin(alpha) |
---|
57 | zpt = radii*np.sin(beta) |
---|
58 | |
---|
59 | return xpt, ypt, zpt |
---|
60 | |
---|
61 | def surface_sphere(radii,Npts): |
---|
62 | """ Function to provide an sphere as matrix of x,y,z coordinates |
---|
63 | radii: radii of the sphere |
---|
64 | Npts: number of points to discretisize longitues (half for latitudes) |
---|
65 | """ |
---|
66 | fname = 'surface_sphere' |
---|
67 | |
---|
68 | sphereup = np.zeros((3,Npts/2,Npts), dtype=np.float) |
---|
69 | spheredown = np.zeros((3,Npts/2,Npts), dtype=np.float) |
---|
70 | for ia in range(Npts): |
---|
71 | alpha = ia*2*np.pi/(Npts-1) |
---|
72 | for ib in range(Npts/2): |
---|
73 | beta = ib*np.pi/(2.*(Npts/2-1)) |
---|
74 | sphereup[:,ib,ia] = position_sphere(radii, alpha, beta) |
---|
75 | for ib in range(Npts/2): |
---|
76 | beta = -ib*np.pi/(2.*(Npts/2-1)) |
---|
77 | spheredown[:,ib,ia] = position_sphere(radii, alpha, beta) |
---|
78 | |
---|
79 | return sphereup, spheredown |
---|
80 | |
---|
81 | def spheric_line(radii,lon,lat): |
---|
82 | """ Function to transform a series of locations in lon, lat coordinates to x,y,z |
---|
83 | over an 3D space |
---|
84 | radii: radius of the sphere |
---|
85 | lon: array of angles along longitudes |
---|
86 | lat: array of angles along latitudes |
---|
87 | """ |
---|
88 | fname = 'spheric_line' |
---|
89 | |
---|
90 | Lint = lon.shape[0] |
---|
91 | coords = np.zeros((Lint,3), dtype=np.float) |
---|
92 | |
---|
93 | for iv in range(Lint): |
---|
94 | coords[iv,:] = position_sphere(radii, lon[iv], lat[iv]) |
---|
95 | |
---|
96 | return coords |
---|
97 | |
---|
98 | ####### ####### ##### #### ### ## # |
---|
99 | # Plotting |
---|
100 | |
---|
101 | def plot_sphere(iazm=-60., iele=30., dist=10., Npts=100, radii=10, \ |
---|
102 | drwsfc=[True,True], colsfc=['#AAAAAA','#646464'], \ |
---|
103 | drwxline = True, linex=[':','b',2.], drwyline = True, liney=[':','r',2.], \ |
---|
104 | drwzline = True, linez=['-.','g',2.], drwxcline=[True,True], \ |
---|
105 | linexc=[['-','#646400',1.],['--','#646400',1.]], \ |
---|
106 | drwequator=[True,True], lineeq=[['-','#AA00AA',1.],['--','#AA00AA',1.]], \ |
---|
107 | drwgreeenwhich=[True,True], linegw=[['-','k',1.],['--','k',1.]]): |
---|
108 | """ Function to plot an sphere and determine which standard lines will be also |
---|
109 | drawn |
---|
110 | iazm: azimut of the camera form the sphere |
---|
111 | iele: elevation of the camera form the sphere |
---|
112 | dist: distance of the camera form the sphere |
---|
113 | Npts: Resolution for the sphere |
---|
114 | radii: radius of the sphere |
---|
115 | drwsfc: whether 'up' and 'down' portions of the sphere should be drawn |
---|
116 | colsfc: colors of the surface of the sphere portions ['up', 'down'] |
---|
117 | drwxline: whether x-axis line should be drawn |
---|
118 | linex: properties of the x-axis line ['type', 'color', 'wdith'] |
---|
119 | drwyline: whether y-axis line should be drawn |
---|
120 | liney: properties of the y-axis line ['type', 'color', 'wdith'] |
---|
121 | drwzline: whether z-axis line should be drawn |
---|
122 | linez: properties of the z-axis line ['type', 'color', 'wdith'] |
---|
123 | drwequator: whether 'front' and 'back' portions of the Equator should be drawn |
---|
124 | lineeq: properties of the lines 'front' and 'back' of the Equator |
---|
125 | drwgreeenwhich: whether 'front', 'back' portions of Greenqhich should be drawn |
---|
126 | linegw: properties of the lines 'front' and 'back' Greenwhich |
---|
127 | drwxcline: whether 'front', 'back' 90 line (lon=90., lon=270.) should be drawn |
---|
128 | linexc: properties of the lines 'front' and 'back' for the 90 line |
---|
129 | """ |
---|
130 | fname = 'plot_sphere' |
---|
131 | |
---|
132 | iazmrad = iazm*np.pi/180. |
---|
133 | ielerad = iele*np.pi/180. |
---|
134 | |
---|
135 | # 3D surface Sphere |
---|
136 | sfcsphereu, sfcsphered = surface_sphere(radii,Npts) |
---|
137 | |
---|
138 | # greenwhich |
---|
139 | if iazmrad > np.pi/2. and iazmrad < 3.*np.pi/2.: |
---|
140 | ia=np.pi-ielerad |
---|
141 | else: |
---|
142 | ia=0.-ielerad |
---|
143 | ea=ia+np.pi |
---|
144 | da = (ea-ia)/(Npts-1) |
---|
145 | beta = np.arange(ia,ea+da,da)[0:Npts] |
---|
146 | alpha = np.zeros((Npts), dtype=np.float) |
---|
147 | greenwhichc = spheric_line(radii,alpha,beta) |
---|
148 | ia=ea+0. |
---|
149 | ea=ia+np.pi |
---|
150 | da = (ea-ia)/(Npts-1) |
---|
151 | beta = np.arange(ia,ea+da,da)[0:Npts] |
---|
152 | greenwhichd = spheric_line(radii,alpha,beta) |
---|
153 | |
---|
154 | # Equator |
---|
155 | ia=np.pi-iazmrad/2. |
---|
156 | ea=ia+np.pi |
---|
157 | da = (ea-ia)/(Npts-1) |
---|
158 | alpha = np.arange(ia,ea+da,da)[0:Npts] |
---|
159 | beta = np.zeros((Npts), dtype=np.float) |
---|
160 | equatorc = spheric_line(radii,alpha,beta) |
---|
161 | ia=ea+0. |
---|
162 | ea=ia+np.pi |
---|
163 | da = (ea-ia)/(Npts-1) |
---|
164 | alpha = np.arange(ia,ea+da,da)[0:Npts] |
---|
165 | equatord = spheric_line(radii,alpha,beta) |
---|
166 | |
---|
167 | # 90 line |
---|
168 | if iazmrad > np.pi and iazmrad < 2.*np.pi: |
---|
169 | ia=3.*np.pi/2. + ielerad |
---|
170 | else: |
---|
171 | ia=np.pi/2. - ielerad |
---|
172 | if ielerad < 0.: |
---|
173 | ia = ia + np.pi |
---|
174 | ea=ia+np.pi |
---|
175 | da = (ea-ia)/(Npts-1) |
---|
176 | beta = np.arange(ia,ea+da,da)[0:Npts] |
---|
177 | alpha = np.ones((Npts), dtype=np.float)*np.pi/2. |
---|
178 | xclinec = spheric_line(radii,alpha,beta) |
---|
179 | ia=ea+0. |
---|
180 | ea=ia+np.pi |
---|
181 | da = (ea-ia)/(Npts-1) |
---|
182 | beta = np.arange(ia,ea+da,da)[0:Npts] |
---|
183 | xclined = spheric_line(radii,alpha,beta) |
---|
184 | |
---|
185 | # x line |
---|
186 | xline = np.zeros((2,3), dtype=np.float) |
---|
187 | xline[0,:] = position_sphere(radii, 0., 0.) |
---|
188 | xline[1,:] = position_sphere(radii, np.pi, 0.) |
---|
189 | |
---|
190 | # y line |
---|
191 | yline = np.zeros((2,3), dtype=np.float) |
---|
192 | yline[0,:] = position_sphere(radii, np.pi/2., 0.) |
---|
193 | yline[1,:] = position_sphere(radii, 3*np.pi/2., 0.) |
---|
194 | |
---|
195 | # z line |
---|
196 | zline = np.zeros((2,3), dtype=np.float) |
---|
197 | zline[0,:] = position_sphere(radii, 0., np.pi/2.) |
---|
198 | zline[1,:] = position_sphere(radii, 0., -np.pi/2.) |
---|
199 | |
---|
200 | fig = plt.figure() |
---|
201 | ax = fig.gca(projection='3d') |
---|
202 | |
---|
203 | # Sphere surface |
---|
204 | if drwsfc[0]: |
---|
205 | ax.plot_surface(sfcsphereu[0,:,:], sfcsphereu[1,:,:], sfcsphereu[2,:,:], \ |
---|
206 | color=colsfc[0]) |
---|
207 | if drwsfc[1]: |
---|
208 | ax.plot_surface(sfcsphered[0,:,:], sfcsphered[1,:,:], sfcsphered[2,:,:], \ |
---|
209 | color=colsfc[1]) |
---|
210 | |
---|
211 | # greenwhich |
---|
212 | linev = linegw[0] |
---|
213 | if drwgreeenwhich[0]: |
---|
214 | ax.plot(greenwhichc[:,0], greenwhichc[:,1], greenwhichc[:,2], linev[0], \ |
---|
215 | color=linev[1], linewidth=linev[2], label='Greenwich') |
---|
216 | linev = linegw[1] |
---|
217 | if drwgreeenwhich[1]: |
---|
218 | ax.plot(greenwhichd[:,0], greenwhichd[:,1], greenwhichd[:,2], linev[0], \ |
---|
219 | color=linev[1], linewidth=linev[2]) |
---|
220 | |
---|
221 | # Equator |
---|
222 | linev = lineeq[0] |
---|
223 | if drwequator[0]: |
---|
224 | ax.plot(equatorc[:,0], equatorc[:,1], equatorc[:,2], linev[0], \ |
---|
225 | color=linev[1], linewidth=linev[2], label='Equator') |
---|
226 | linev = lineeq[1] |
---|
227 | if drwequator[1]: |
---|
228 | ax.plot(equatord[:,0], equatord[:,1], equatord[:,2], linev[0], \ |
---|
229 | color=linev[1], linewidth=linev[2]) |
---|
230 | |
---|
231 | # 90line |
---|
232 | linev = linexc[0] |
---|
233 | if drwxcline[0]: |
---|
234 | ax.plot(xclinec[:,0], xclinec[:,1], xclinec[:,2], linev[0], color=linev[1], \ |
---|
235 | linewidth=linev[2], label='90-line') |
---|
236 | linev = linexc[1] |
---|
237 | if drwxcline[1]: |
---|
238 | ax.plot(xclined[:,0], xclined[:,1], xclined[:,2], linev[0], color=linev[1], \ |
---|
239 | linewidth=linev[2]) |
---|
240 | |
---|
241 | # x line |
---|
242 | linev = linex |
---|
243 | if drwxline: |
---|
244 | ax.plot([xline[0,0],xline[1,0]], [xline[0,1],xline[1,1]], \ |
---|
245 | [xline[0,2],xline[1,2]], linev[0], color=linev[1], linewidth=linev[2], label='xline') |
---|
246 | |
---|
247 | # y line |
---|
248 | linev = liney |
---|
249 | if drwyline: |
---|
250 | ax.plot([yline[0,0],yline[1,0]], [yline[0,1],yline[1,1]], \ |
---|
251 | [yline[0,2],yline[1,2]], linev[0], color=linev[1], linewidth=linev[2], label='yline') |
---|
252 | |
---|
253 | # z line |
---|
254 | linev = linez |
---|
255 | if drwzline: |
---|
256 | ax.plot([zline[0,0],zline[1,0]], [zline[0,1],zline[1,1]], \ |
---|
257 | [zline[0,2],zline[1,2]], linev[0], color=linev[1], linewidth=linev[2], label='zline') |
---|
258 | |
---|
259 | plt.legend() |
---|
260 | |
---|
261 | return fig, ax |
---|
262 | |
---|