source: lmdz_wrf/trunk/tools/geometry_tools.py @ 2411

Last change on this file since 2411 was 2411, checked in by lfita, 6 years ago

Adding:

  • `generic_tools.py': Script for geometry calculations and operations as well as definition of different standard objects and shapes
File size: 9.3 KB
Line 
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
12import numpy as np
13import matplotlib as mpl
14from mpl_toolkits.mplot3d import Axes3D
15import 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
29def 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
44def 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
61def 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
81def 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
101def 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
Note: See TracBrowser for help on using the repository browser.