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

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

Woring on 'ang_circ'
File size: 17.1 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
17errormsg = 'ERROR -- error -- ERROR -- error'
18
19####### Contents:
20# deg_deci: Function to pass from degrees [deg, minute, sec] to decimal angles [rad]
21# dist_points: Function to provide the distance between two points
22# multi_rotate_2D: Function to rotate multiple vectors by a certain angle in the plain
23# position_sphere: Function to tranform fom a point in lon, lat deg coordinates to
24#   cartesian coordinates over an sphere
25# rotate_2D: Function to rotate a vector by a certain angle in the plain
26# rotate_line2D: Function to rotate a line given by 2 pairs of x,y coordinates by a
27#   certain angle in the plain
28# rotate_lines2D: Function to rotate multiple lines given by mulitple pars of x,y
29#   coordinates by a certain angle in the plain
30# spheric_line: Function to transform a series of locations in lon, lat coordinates
31#   to x,y,z over an 3D spaceFunction to provide coordinates of a line  on a 3D space
32
33## Shapes/objects
34# ellipse_polar: Function to determine an ellipse from its center and polar coordinates
35# surface_sphere: Function to provide an sphere as matrix of x,y,z coordinates
36
37## Plotting
38# plot_sphere: Function to plot an sphere and determine which standard lines will be
39#   also drawn
40
41def deg_deci(angle):
42    """ Function to pass from degrees [deg, minute, sec] to decimal angles [rad]
43      angle: list of [deg, minute, sec] to pass
44    >>> deg_deci([41., 58., 34.])
45    0.732621346072
46    """
47    fname = 'deg_deci'
48
49    deg = np.abs(angle[0]) + np.abs(angle[1])/60. + np.abs(angle[2])/3600.
50
51    if angle[0] < 0.: deg = -deg*np.pi/180.
52    else: deg = deg*np.pi/180.
53
54    return deg
55
56def position_sphere(radii, alpha, beta):
57    """ Function to tranform fom a point in lon, lat deg coordinates to cartesian 
58          coordinates over an sphere
59      radii: radii of the sphere
60      alpha: longitude of the point
61      beta: latitude of the point
62    >>> position_sphere(10., 30., 45.)
63    (0.81031678432964027, -5.1903473778327376, 8.5090352453411846
64    """
65    fname = 'position_sphere'
66
67    xpt = radii*np.cos(beta)*np.cos(alpha)
68    ypt = radii*np.cos(beta)*np.sin(alpha)
69    zpt = radii*np.sin(beta)
70
71    return xpt, ypt, zpt
72
73def spheric_line(radii,lon,lat):
74    """ Function to transform a series of locations in lon, lat coordinates to x,y,z
75          over an 3D space
76      radii: radius of the sphere
77      lon: array of angles along longitudes
78      lat: array of angles along latitudes
79    """
80    fname = 'spheric_line'
81
82    Lint = lon.shape[0]
83    coords = np.zeros((Lint,3), dtype=np.float)
84
85    for iv in range(Lint):
86        coords[iv,:] = position_sphere(radii, lon[iv], lat[iv])
87
88    return coords
89
90def rotate_2D(vector, angle):
91    """ Function to rotate a vector by a certain angle in the plain
92      vector= vector to rotate [y, x]
93      angle= angle to rotate [rad]
94    >>> rotate_2D(np.array([1.,0.]), np.pi/4.)
95    [ 0.70710678 -0.70710678]
96    """
97    fname = 'rotate_2D'
98
99    rotmat = np.zeros((2,2), dtype=np.float)
100
101    rotmat[0,0] = np.cos(angle)
102    rotmat[0,1] = -np.sin(angle)
103    rotmat[1,0] = np.sin(angle)
104    rotmat[1,1] = np.cos(angle)
105
106    rotvector = np.zeros((2), dtype=np.float)
107
108    vecv = np.zeros((2), dtype=np.float)
109
110    # Unifying vector
111    modvec = vector[0]**2+vector[1]**2
112    if modvec != 0: 
113        vecv[0] = vector[1]/modvec
114        vecv[1] = vector[0]/modvec
115
116        rotvec = np.matmul(rotmat, vecv)
117        rotvec = np.where(np.abs(rotvec) < 1.e-7, 0., rotvec)
118
119        rotvector[0] = rotvec[1]*modvec
120        rotvector[1] = rotvec[0]*modvec
121
122    return rotvector
123
124def multi_rotate_2D(vectors, angle):
125    """ Function to rotate multiple vectors by a certain angle in the plain
126      line= matrix of vectors to rotate
127      angle= angle to rotate [rad]
128    >>> square = np.zeros((4,2), dtype=np.float)
129    >>> square[0,:] = [-0.5,-0.5]
130    >>> square[1,:] = [0.5,-0.5]
131    >>> square[2,:] = [0.5,0.5]
132    >>> square[3,:] = [-0.5,0.5]
133    >>> multi_rotate_2D(square, np.pi/4.)
134    [[-0.70710678  0.        ]
135     [ 0.         -0.70710678]
136     [ 0.70710678  0.        ]
137     [ 0.          0.70710678]]
138    """
139    fname = 'multi_rotate_2D'
140
141    rotvecs = np.zeros(vectors.shape, dtype=np.float)
142
143    Nvecs = vectors.shape[0]
144    for iv in range(Nvecs):
145        rotvecs[iv,:] = rotate_2D(vectors[iv,:], angle)
146
147    return rotvecs
148
149def rotate_line2D(line, angle):
150    """ Function to rotate a line given by 2 pairs of x,y coordinates by a certain
151          angle in the plain
152      line= line to rotate as couple of points [[y0,x0], [y1,x1]]
153      angle= angle to rotate [rad]
154    >>> rotate_line2D(np.array([[0.,0.], [1.,0.]]), np.pi/4.)
155    [[ 0.          0.        ]
156     [0.70710678  -0.70710678]]
157    """
158    fname = 'rotate_2D'
159
160    rotline = np.zeros((2,2), dtype=np.float)
161    rotline[0,:] = rotate_2D(line[0,:], angle)
162    rotline[1,:] = rotate_2D(line[1,:], angle)
163
164    return rotline
165
166def rotate_lines2D(lines, angle):
167    """ Function to rotate multiple lines given by mulitple pars of x,y coordinates 
168          by a certain angle in the plain
169      line= matrix of N couples of points [N, [y0,x0], [y1,x1]]
170      angle= angle to rotate [rad]
171    >>> square = np.zeros((4,2,2), dtype=np.float)
172    >>> square[0,0,:] = [-0.5,-0.5]
173    >>> square[0,1,:] = [0.5,-0.5]
174    >>> square[1,0,:] = [0.5,-0.5]
175    >>> square[1,1,:] = [0.5,0.5]
176    >>> square[2,0,:] = [0.5,0.5]
177    >>> square[2,1,:] = [-0.5,0.5]
178    >>> square[3,0,:] = [-0.5,0.5]
179    >>> square[3,1,:] = [-0.5,-0.5]
180    >>> rotate_lines2D(square, np.pi/4.)
181    [[[-0.70710678  0.        ]
182      [ 0.         -0.70710678]]
183
184     [[ 0.         -0.70710678]
185      [ 0.70710678  0.        ]]
186
187     [[ 0.70710678  0.        ]
188      [ 0.          0.70710678]]
189
190     [[ 0.          0.70710678]
191      [-0.70710678  0.        ]]]
192    """
193    fname = 'rotate_lines2D'
194
195    rotlines = np.zeros(lines.shape, dtype=np.float)
196
197    Nlines = lines.shape[0]
198    for il in range(Nlines):
199        line = np.zeros((2,2), dtype=np.float)
200        line[0,:] = lines[il,0,:]
201        line[1,:] = lines[il,1,:]
202
203        rotlines[il,:,:] = rotate_line2D(line, angle)
204
205    return rotlines
206
207def dist_points(ptA, ptB):
208    """ Function to provide the distance between two points
209      ptA: coordinates of the point A [yA, xA]
210      ptB: coordinates of the point B [yB, xB]
211    >>> dist_points([1.,1.], [-1.,-1.])
212    2.82842712475
213    """
214    fname = 'dist_points'
215
216    dist = np.sqrt( (ptA[0]-ptB[0])**2 + (ptA[1]-ptB[1])**2)
217
218    return dist
219
220####### ###### ##### #### ### ## #
221# Shapes/objects
222
223def surface_sphere(radii,Npts):
224    """ Function to provide an sphere as matrix of x,y,z coordinates
225      radii: radii of the sphere
226      Npts: number of points to discretisize longitues (half for latitudes)
227    """
228    fname = 'surface_sphere'
229
230    sphereup = np.zeros((3,Npts/2,Npts), dtype=np.float)
231    spheredown = np.zeros((3,Npts/2,Npts), dtype=np.float)
232    for ia in range(Npts):
233        alpha = ia*2*np.pi/(Npts-1)
234        for ib in range(Npts/2):
235            beta = ib*np.pi/(2.*(Npts/2-1))
236            sphereup[:,ib,ia] = position_sphere(radii, alpha, beta)
237        for ib in range(Npts/2):
238            beta = -ib*np.pi/(2.*(Npts/2-1))
239            spheredown[:,ib,ia] = position_sphere(radii, alpha, beta)
240
241    return sphereup, spheredown
242
243def ellipse_polar(c, a, b, Nang=100):
244    """ Function to determine an ellipse from its center and polar coordinates
245        FROM: https://en.wikipedia.org/wiki/Ellipse
246      c= coordinates of the center
247      a= distance major axis
248      b= distance minor axis
249      Nang= number of angles to use
250    """
251    fname = 'ellipse_polar'
252
253    if np.mod(Nang,2) == 0: Nang=Nang+1
254 
255    dtheta = 2*np.pi/(Nang-1)
256
257    ellipse = np.zeros((Nang,2), dtype=np.float)
258    for ia in range(Nang):
259        theta = dtheta*ia
260        rad = a*b/np.sqrt( (b*np.cos(theta))**2 + (a*np.sin(theta))**2 )
261        x = rad*np.cos(theta)
262        y = rad*np.sin(theta)
263        ellipse[ia,:] = [y+c[0],x+c[1]]
264
265    return ellipse
266
267def hyperbola_polar(a, b, Nang=100):
268    """ Fcuntion to determine an hyperbola in polar coordinates
269        FROM: https://en.wikipedia.org/wiki/Hyperbola#Polar_coordinates
270          x^2/a^2 - y^2/b^2 = 1
271      a= x-parameter
272      y= y-parameter
273      Nang= number of angles to use
274      DOES NOT WORK!!!!
275    """
276    fname = 'hyperbola_polar'
277
278    dtheta = 2.*np.pi/(Nang-1)
279
280    # Positive branch
281    hyperbola_p = np.zeros((Nang,2), dtype=np.float)
282    for ia in range(Nang):
283        theta = dtheta*ia
284        x = a*np.cosh(theta)
285        y = b*np.sinh(theta)
286        hyperbola_p[ia,:] = [y,x]
287
288    # Negative branch
289    hyperbola_n = np.zeros((Nang,2), dtype=np.float)
290    for ia in range(Nang):
291        theta = dtheta*ia
292        x = -a*np.cosh(theta)
293        y = b*np.sinh(theta)
294        hyperbola_n[ia,:] = [y,x]
295
296    return hyperbola_p, hyperbola_n
297
298def circ_sec(ptA, ptB, radii, Nang=100):
299    """ Function union of point A and B by a section of a circle
300      ptA= coordinates od the point A [yA, xA]
301      ptB= coordinates od the point B [yB, xB]
302      radii= radi of the circle to use to unite the points
303      Nang= amount of angles to use
304    """
305    fname = 'circ_sec'
306
307    distAB = dist_points(ptA,ptB)
308
309    if distAB > radii:
310        print errormsg
311        print '  ' + fname + ': radii=', radii, " too small for the distance " +     \
312          "between points !!"
313        print '    distance between points:', distAB
314        quit(-1)
315
316    # Coordinate increments
317    dAB = np.abs(ptA-ptB)
318
319    # angle of the circular section joining points
320    tottheta = 2.*np.arcsin(distAB/2./radii)
321
322    # center along coincident bisection of the union
323    xcc = -radii
324    ycc = 0.
325
326    # Angle of the points
327    pttheta = np.arctan2(ptB[0]-ptA[0],ptB[1]-ptA[1])
328
329    # rotating the position of the center
330    rotang = np.pi/2.-pttheta
331
332    newcc = rotate_line2D(np.array([[ycc,xcc], [0.,0.]]), rotang)
333    print newcc
334
335    quit()
336
337    dtheta = np.abs(tottheta)/(Nang-1)
338    if sign == 1:
339        dtheta = dtheta*(-1.)
340
341    print 'Lluis tottheta:', tottheta*180./np.pi, 'dtheta:', dtheta*180./np.pi, 'c:', xc,yc
342
343    circ_sec = np.zeros((Nang,2), dtype=np.float)
344    for ia in range(Nang):
345        theta = dtheta*ia
346        x = radii*np.cos(theta)
347        y = radii*np.sin(theta)
348
349        circ_sec[ia,:] = [y+yc,x+xc]
350        print ia, 'Lluis xy:', x,y, 'circ', circ_sec[ia,:]
351   
352    return circ_sec
353
354# FROM: http://www.photographers1.com/Sailing/NauticalTerms&Nomenclature.html
355def zsailing_boat(length=10., beam=3., sternbp=0.5):
356    """ Function to define an schematic boat from the z-plane
357      length: length of the boat
358      beam: beam of the boat
359      sternbp: beam at stern as percentage of beam
360    """
361    fname = 'zsailing_boat'
362
363   
364
365    return boat
366
367####### ####### ##### #### ### ## #
368# Plotting
369
370def plot_sphere(iazm=-60., iele=30., dist=10., Npts=100, radii=10,                   \
371  drwsfc=[True,True], colsfc=['#AAAAAA','#646464'],                                  \
372  drwxline = True, linex=[':','b',2.], drwyline = True, liney=[':','r',2.],          \
373  drwzline = True, linez=['-.','g',2.], drwxcline=[True,True],                       \
374  linexc=[['-','#646400',1.],['--','#646400',1.]],                                   \
375  drwequator=[True,True], lineeq=[['-','#AA00AA',1.],['--','#AA00AA',1.]],           \
376  drwgreeenwhich=[True,True], linegw=[['-','k',1.],['--','k',1.]]):
377    """ Function to plot an sphere and determine which standard lines will be also
378        drawn
379      iazm: azimut of the camera form the sphere
380      iele: elevation of the camera form the sphere
381      dist: distance of the camera form the sphere
382      Npts: Resolution for the sphere
383      radii: radius of the sphere
384      drwsfc: whether 'up' and 'down' portions of the sphere should be drawn
385      colsfc: colors of the surface of the sphere portions ['up', 'down']
386      drwxline: whether x-axis line should be drawn
387      linex: properties of the x-axis line ['type', 'color', 'wdith']
388      drwyline: whether y-axis line should be drawn
389      liney: properties of the y-axis line ['type', 'color', 'wdith']
390      drwzline: whether z-axis line should be drawn
391      linez: properties of the z-axis line ['type', 'color', 'wdith']
392      drwequator: whether 'front' and 'back' portions of the Equator should be drawn
393      lineeq: properties of the lines 'front' and 'back' of the Equator
394      drwgreeenwhich: whether 'front', 'back' portions of Greenqhich should be drawn
395      linegw: properties of the lines 'front' and 'back' Greenwhich
396      drwxcline: whether 'front', 'back' 90 line (lon=90., lon=270.) should be drawn
397      linexc: properties of the lines 'front' and 'back' for the 90 line
398    """
399    fname = 'plot_sphere'
400
401    iazmrad = iazm*np.pi/180.
402    ielerad = iele*np.pi/180.
403
404    # 3D surface Sphere
405    sfcsphereu, sfcsphered = surface_sphere(radii,Npts)
406   
407    # greenwhich
408    if iazmrad > np.pi/2. and iazmrad < 3.*np.pi/2.:
409        ia=np.pi-ielerad
410    else:
411        ia=0.-ielerad
412    ea=ia+np.pi
413    da = (ea-ia)/(Npts-1)
414    beta = np.arange(ia,ea+da,da)[0:Npts]
415    alpha = np.zeros((Npts), dtype=np.float)
416    greenwhichc = spheric_line(radii,alpha,beta)
417    ia=ea+0.
418    ea=ia+np.pi
419    da = (ea-ia)/(Npts-1)
420    beta = np.arange(ia,ea+da,da)[0:Npts]
421    greenwhichd = spheric_line(radii,alpha,beta)
422
423    # Equator
424    ia=np.pi-iazmrad/2.
425    ea=ia+np.pi
426    da = (ea-ia)/(Npts-1)
427    alpha = np.arange(ia,ea+da,da)[0:Npts]
428    beta = np.zeros((Npts), dtype=np.float)
429    equatorc = spheric_line(radii,alpha,beta)
430    ia=ea+0.
431    ea=ia+np.pi
432    da = (ea-ia)/(Npts-1)
433    alpha = np.arange(ia,ea+da,da)[0:Npts]
434    equatord = spheric_line(radii,alpha,beta)
435
436    # 90 line
437    if iazmrad > np.pi and iazmrad < 2.*np.pi:
438        ia=3.*np.pi/2. + ielerad
439    else:
440        ia=np.pi/2. - ielerad
441    if ielerad < 0.:
442        ia = ia + np.pi
443    ea=ia+np.pi
444    da = (ea-ia)/(Npts-1)
445    beta = np.arange(ia,ea+da,da)[0:Npts]
446    alpha = np.ones((Npts), dtype=np.float)*np.pi/2.
447    xclinec = spheric_line(radii,alpha,beta)
448    ia=ea+0.
449    ea=ia+np.pi
450    da = (ea-ia)/(Npts-1)
451    beta = np.arange(ia,ea+da,da)[0:Npts]
452    xclined = spheric_line(radii,alpha,beta)
453
454    # x line
455    xline = np.zeros((2,3), dtype=np.float)
456    xline[0,:] = position_sphere(radii, 0., 0.)
457    xline[1,:] = position_sphere(radii, np.pi, 0.)
458
459    # y line
460    yline = np.zeros((2,3), dtype=np.float)
461    yline[0,:] = position_sphere(radii, np.pi/2., 0.)
462    yline[1,:] = position_sphere(radii, 3*np.pi/2., 0.)
463
464    # z line
465    zline = np.zeros((2,3), dtype=np.float)
466    zline[0,:] = position_sphere(radii, 0., np.pi/2.)
467    zline[1,:] = position_sphere(radii, 0., -np.pi/2.)
468
469    fig = plt.figure()
470    ax = fig.gca(projection='3d')
471
472    # Sphere surface
473    if drwsfc[0]:
474        ax.plot_surface(sfcsphereu[0,:,:], sfcsphereu[1,:,:], sfcsphereu[2,:,:],     \
475          color=colsfc[0])
476    if drwsfc[1]:
477        ax.plot_surface(sfcsphered[0,:,:], sfcsphered[1,:,:], sfcsphered[2,:,:],     \
478          color=colsfc[1])
479
480    # greenwhich
481    linev = linegw[0]
482    if drwgreeenwhich[0]:
483        ax.plot(greenwhichc[:,0], greenwhichc[:,1], greenwhichc[:,2], linev[0],      \
484          color=linev[1], linewidth=linev[2],  label='Greenwich')
485    linev = linegw[1]
486    if drwgreeenwhich[1]:
487        ax.plot(greenwhichd[:,0], greenwhichd[:,1], greenwhichd[:,2], linev[0],      \
488          color=linev[1], linewidth=linev[2])
489
490    # Equator
491    linev = lineeq[0]
492    if drwequator[0]:
493        ax.plot(equatorc[:,0], equatorc[:,1], equatorc[:,2], linev[0],               \
494          color=linev[1], linewidth=linev[2], label='Equator')
495    linev = lineeq[1]
496    if drwequator[1]:
497        ax.plot(equatord[:,0], equatord[:,1], equatord[:,2], linev[0],               \
498          color=linev[1], linewidth=linev[2])
499
500    # 90line
501    linev = linexc[0]
502    if drwxcline[0]:
503        ax.plot(xclinec[:,0], xclinec[:,1], xclinec[:,2], linev[0], color=linev[1],  \
504          linewidth=linev[2], label='90-line')
505    linev = linexc[1]
506    if drwxcline[1]:
507        ax.plot(xclined[:,0], xclined[:,1], xclined[:,2], linev[0], color=linev[1],  \
508          linewidth=linev[2])
509
510    # x line
511    linev = linex
512    if drwxline:
513        ax.plot([xline[0,0],xline[1,0]], [xline[0,1],xline[1,1]],                    \
514          [xline[0,2],xline[1,2]], linev[0], color=linev[1], linewidth=linev[2],  label='xline')
515
516    # y line
517    linev = liney
518    if drwyline:
519        ax.plot([yline[0,0],yline[1,0]], [yline[0,1],yline[1,1]],                    \
520          [yline[0,2],yline[1,2]], linev[0], color=linev[1], linewidth=linev[2],  label='yline')
521
522    # z line
523    linev = linez
524    if drwzline:
525        ax.plot([zline[0,0],zline[1,0]], [zline[0,1],zline[1,1]],                    \
526          [zline[0,2],zline[1,2]], linev[0], color=linev[1], linewidth=linev[2],  label='zline')
527
528    plt.legend()
529
530    return fig, ax
531
Note: See TracBrowser for help on using the repository browser.