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

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

Re-ordering documentation
File size: 36.8 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
16import os
17import generic_tools as gen
18
19errormsg = 'ERROR -- error -- ERROR -- error'
20infmsg = 'INFORMATION -- information -- INFORMATION -- information'
21
22####### Contents:
23# deg_deci: Function to pass from degrees [deg, minute, sec] to decimal angles [rad]
24# dist_points: Function to provide the distance between two points
25# max_coords_poly: Function to provide the extremes of the coordinates of a polygon
26# mirror_polygon: Function to reflex a polygon for a given axis
27# position_sphere: Function to tranform fom a point in lon, lat deg coordinates to
28#   cartesian coordinates over an sphere
29# read_join_poly: Function to read an ASCII file with the combination of polygons
30# rotate_2D: Function to rotate a vector by a certain angle in the plain
31# rotate_polygon_2D: Function to rotate 2D plain the vertices of a polygon
32# rotate_line2D: Function to rotate a line given by 2 pairs of x,y coordinates by a
33#   certain angle in the plain
34# rotate_lines2D: Function to rotate multiple lines given by mulitple pars of x,y
35#   coordinates by a certain angle in the plain
36# spheric_line: Function to transform a series of locations in lon, lat coordinates
37#   to x,y,z over an 3D spaceFunction to provide coordinates of a line  on a 3D space
38# write_join_poly: Function to write an ASCII file with the combination of polygons
39
40## Shapes/objects
41# circ_sec: Function union of point A and B by a section of a circle
42# ellipse_polar: Function to determine an ellipse from its center and polar coordinates
43# p_doubleArrow: Function to provide an arrow with double lines
44# p_circle: Function to get a polygon of a circle
45# p_reg_polygon: Function to provide a regular polygon of Nv vertices
46# p_reg_star: Function to provide a regular star of Nv vertices
47# p_sinusiode: Function to get coordinates of a sinusoidal curve
48# p_square: Function to get a polygon square
49# p_spiral: Function to provide a polygon of an Archimedean spiral
50# p_triangle: Function to provide the polygon of a triangle from its 3 vertices
51# surface_sphere: Function to provide an sphere as matrix of x,y,z coordinates
52# z_boat: Function to define an schematic boat from the z-plane
53# zsailing_boat: Function to define an schematic sailing boat from the z-plane with sails
54
55## Plotting
56# plot_sphere: Function to plot an sphere and determine which standard lines will be
57#   also drawn
58
59def deg_deci(angle):
60    """ Function to pass from degrees [deg, minute, sec] to decimal angles [rad]
61      angle: list of [deg, minute, sec] to pass
62    >>> deg_deci([41., 58., 34.])
63    0.732621346072
64    """
65    fname = 'deg_deci'
66
67    deg = np.abs(angle[0]) + np.abs(angle[1])/60. + np.abs(angle[2])/3600.
68
69    if angle[0] < 0.: deg = -deg*np.pi/180.
70    else: deg = deg*np.pi/180.
71
72    return deg
73
74def position_sphere(radii, alpha, beta):
75    """ Function to tranform fom a point in lon, lat deg coordinates to cartesian 
76          coordinates over an sphere
77      radii: radii of the sphere
78      alpha: longitude of the point
79      beta: latitude of the point
80    >>> position_sphere(10., 30., 45.)
81    (0.81031678432964027, -5.1903473778327376, 8.5090352453411846
82    """
83    fname = 'position_sphere'
84
85    xpt = radii*np.cos(beta)*np.cos(alpha)
86    ypt = radii*np.cos(beta)*np.sin(alpha)
87    zpt = radii*np.sin(beta)
88
89    return xpt, ypt, zpt
90
91def spheric_line(radii,lon,lat):
92    """ Function to transform a series of locations in lon, lat coordinates to x,y,z
93          over an 3D space
94      radii: radius of the sphere
95      lon: array of angles along longitudes
96      lat: array of angles along latitudes
97    """
98    fname = 'spheric_line'
99
100    Lint = lon.shape[0]
101    coords = np.zeros((Lint,3), dtype=np.float)
102
103    for iv in range(Lint):
104        coords[iv,:] = position_sphere(radii, lon[iv], lat[iv])
105
106    return coords
107
108def rotate_2D(vector, angle):
109    """ Function to rotate a vector by a certain angle in the plain
110      vector= vector to rotate [y, x]
111      angle= angle to rotate [rad]
112    >>> rotate_2D(np.array([1.,0.]), np.pi/4.)
113    [ 0.70710678 -0.70710678]
114    """
115    fname = 'rotate_2D'
116
117    rotmat = np.zeros((2,2), dtype=np.float)
118
119    rotmat[0,0] = np.cos(angle)
120    rotmat[0,1] = -np.sin(angle)
121    rotmat[1,0] = np.sin(angle)
122    rotmat[1,1] = np.cos(angle)
123
124    rotvector = np.zeros((2), dtype=np.float)
125
126    vecv = np.zeros((2), dtype=np.float)
127
128    # Unifying vector
129    modvec = vector[0]**2+vector[1]**2
130    if modvec != 0: 
131        vecv[0] = vector[1]/modvec
132        vecv[1] = vector[0]/modvec
133
134        rotvec = np.matmul(rotmat, vecv)
135        rotvec = np.where(np.abs(rotvec) < 1.e-7, 0., rotvec)
136
137        rotvector[0] = rotvec[1]*modvec
138        rotvector[1] = rotvec[0]*modvec
139
140    return rotvector
141
142def rotate_polygon_2D(vectors, angle):
143    """ Function to rotate 2D plain the vertices of a polygon
144      line= matrix of vectors to rotate
145      angle= angle to rotate [rad]
146    >>> square = np.zeros((4,2), dtype=np.float)
147    >>> square[0,:] = [-0.5,-0.5]
148    >>> square[1,:] = [0.5,-0.5]
149    >>> square[2,:] = [0.5,0.5]
150    >>> square[3,:] = [-0.5,0.5]
151    >>> rotate_polygon_2D(square, np.pi/4.)
152    [[-0.70710678  0.        ]
153     [ 0.         -0.70710678]
154     [ 0.70710678  0.        ]
155     [ 0.          0.70710678]]
156    """
157    fname = 'rotate_polygon_2D'
158
159    rotvecs = np.zeros(vectors.shape, dtype=np.float)
160
161    Nvecs = vectors.shape[0]
162    for iv in range(Nvecs):
163        rotvecs[iv,:] = rotate_2D(vectors[iv,:], angle)
164
165    return rotvecs
166
167def rotate_line2D(line, angle):
168    """ Function to rotate a line given by 2 pairs of x,y coordinates by a certain
169          angle in the plain
170      line= line to rotate as couple of points [[y0,x0], [y1,x1]]
171      angle= angle to rotate [rad]
172    >>> rotate_line2D(np.array([[0.,0.], [1.,0.]]), np.pi/4.)
173    [[ 0.          0.        ]
174     [0.70710678  -0.70710678]]
175    """
176    fname = 'rotate_2D'
177
178    rotline = np.zeros((2,2), dtype=np.float)
179    rotline[0,:] = rotate_2D(line[0,:], angle)
180    rotline[1,:] = rotate_2D(line[1,:], angle)
181
182    return rotline
183
184def rotate_lines2D(lines, angle):
185    """ Function to rotate multiple lines given by mulitple pars of x,y coordinates 
186          by a certain angle in the plain
187      line= matrix of N couples of points [N, [y0,x0], [y1,x1]]
188      angle= angle to rotate [rad]
189    >>> square = np.zeros((4,2,2), dtype=np.float)
190    >>> square[0,0,:] = [-0.5,-0.5]
191    >>> square[0,1,:] = [0.5,-0.5]
192    >>> square[1,0,:] = [0.5,-0.5]
193    >>> square[1,1,:] = [0.5,0.5]
194    >>> square[2,0,:] = [0.5,0.5]
195    >>> square[2,1,:] = [-0.5,0.5]
196    >>> square[3,0,:] = [-0.5,0.5]
197    >>> square[3,1,:] = [-0.5,-0.5]
198    >>> rotate_lines2D(square, np.pi/4.)
199    [[[-0.70710678  0.        ]
200      [ 0.         -0.70710678]]
201
202     [[ 0.         -0.70710678]
203      [ 0.70710678  0.        ]]
204
205     [[ 0.70710678  0.        ]
206      [ 0.          0.70710678]]
207
208     [[ 0.          0.70710678]
209      [-0.70710678  0.        ]]]
210    """
211    fname = 'rotate_lines2D'
212
213    rotlines = np.zeros(lines.shape, dtype=np.float)
214
215    Nlines = lines.shape[0]
216    for il in range(Nlines):
217        line = np.zeros((2,2), dtype=np.float)
218        line[0,:] = lines[il,0,:]
219        line[1,:] = lines[il,1,:]
220
221        rotlines[il,:,:] = rotate_line2D(line, angle)
222
223    return rotlines
224
225def dist_points(ptA, ptB):
226    """ Function to provide the distance between two points
227      ptA: coordinates of the point A [yA, xA]
228      ptB: coordinates of the point B [yB, xB]
229    >>> dist_points([1.,1.], [-1.,-1.])
230    2.82842712475
231    """
232    fname = 'dist_points'
233
234    dist = np.sqrt( (ptA[0]-ptB[0])**2 + (ptA[1]-ptB[1])**2)
235
236    return dist
237
238def max_coords_poly(polygon):
239    """ Function to provide the extremes of the coordinates of a polygon
240      polygon: coordinates [Nvertexs, 2] of a polygon
241    >>> square = np.zeros((4,2), dtype=np.float)
242    >>> square[0,:] = [-0.5,-0.5]
243    >>> square[1,:] = [0.5,-0.5]
244    >>> square[2,:] = [0.5,0.5]
245    >>> square[3,:] = [-0.5,0.5]
246    >>> max_coords_poly(square)
247    [-0.5, 0.5], [-0.5, 0.5], [0.5, 0.5], 0.5
248    """
249    fname = 'max_coords_poly'
250
251    # x-coordinate min/max
252    nx = np.min(polygon[:,1])
253    xx = np.max(polygon[:,1])
254
255    # y-coordinate min/max
256    ny = np.min(polygon[:,0])
257    xy = np.max(polygon[:,0])
258
259    # x/y-coordinate maximum of absolute values
260    axx = np.max(np.abs(polygon[:,1]))
261    ayx = np.max(np.abs(polygon[:,0]))
262
263    # absolute maximum
264    xyx = np.max([axx, ayx])
265
266    return [nx, xx], [ny, xy], [ayx, axx], xyx
267
268def mirror_polygon(polygon,axis):
269    """ Function to reflex a polygon for a given axis
270      polygon: polygon to mirror
271      axis: axis at which mirror is located ('x' or 'y')
272    """
273    fname = 'mirror_polygon'
274
275    reflex = np.zeros(polygon.shape, dtype=np.float)
276
277    N = polygon.shape[0]
278    if axis == 'x':
279        for iv in range(N):
280            reflex[iv,:] = [-polygon[iv,0], polygon[iv,1]]
281    elif axis == 'y':
282        for iv in range(N):
283            reflex[iv,:] = [polygon[iv,0], -polygon[iv,1]]
284
285    return reflex
286
287####### ###### ##### #### ### ## #
288# Shapes/objects
289
290def surface_sphere(radii,Npts):
291    """ Function to provide an sphere as matrix of x,y,z coordinates
292      radii: radii of the sphere
293      Npts: number of points to discretisize longitues (half for latitudes)
294    """
295    fname = 'surface_sphere'
296
297    sphereup = np.zeros((3,Npts/2,Npts), dtype=np.float)
298    spheredown = np.zeros((3,Npts/2,Npts), dtype=np.float)
299    for ia in range(Npts):
300        alpha = ia*2*np.pi/(Npts-1)
301        for ib in range(Npts/2):
302            beta = ib*np.pi/(2.*(Npts/2-1))
303            sphereup[:,ib,ia] = position_sphere(radii, alpha, beta)
304        for ib in range(Npts/2):
305            beta = -ib*np.pi/(2.*(Npts/2-1))
306            spheredown[:,ib,ia] = position_sphere(radii, alpha, beta)
307
308    return sphereup, spheredown
309
310def ellipse_polar(c, a, b, Nang=100):
311    """ Function to determine an ellipse from its center and polar coordinates
312        FROM: https://en.wikipedia.org/wiki/Ellipse
313      c= coordinates of the center
314      a= distance major axis
315      b= distance minor axis
316      Nang= number of angles to use
317    """
318    fname = 'ellipse_polar'
319
320    if np.mod(Nang,2) == 0: Nang=Nang+1
321 
322    dtheta = 2*np.pi/(Nang-1)
323
324    ellipse = np.zeros((Nang,2), dtype=np.float)
325    for ia in range(Nang):
326        theta = dtheta*ia
327        rad = a*b/np.sqrt( (b*np.cos(theta))**2 + (a*np.sin(theta))**2 )
328        x = rad*np.cos(theta)
329        y = rad*np.sin(theta)
330        ellipse[ia,:] = [y+c[0],x+c[1]]
331
332    return ellipse
333
334def hyperbola_polar(a, b, Nang=100):
335    """ Fcuntion to determine an hyperbola in polar coordinates
336        FROM: https://en.wikipedia.org/wiki/Hyperbola#Polar_coordinates
337          x^2/a^2 - y^2/b^2 = 1
338      a= x-parameter
339      y= y-parameter
340      Nang= number of angles to use
341      DOES NOT WORK!!!!
342    """
343    fname = 'hyperbola_polar'
344
345    dtheta = 2.*np.pi/(Nang-1)
346
347    # Positive branch
348    hyperbola_p = np.zeros((Nang,2), dtype=np.float)
349    for ia in range(Nang):
350        theta = dtheta*ia
351        x = a*np.cosh(theta)
352        y = b*np.sinh(theta)
353        hyperbola_p[ia,:] = [y,x]
354
355    # Negative branch
356    hyperbola_n = np.zeros((Nang,2), dtype=np.float)
357    for ia in range(Nang):
358        theta = dtheta*ia
359        x = -a*np.cosh(theta)
360        y = b*np.sinh(theta)
361        hyperbola_n[ia,:] = [y,x]
362
363    return hyperbola_p, hyperbola_n
364
365def circ_sec(ptA, ptB, radii, Nang=100):
366    """ Function union of point A and B by a section of a circle
367      ptA= coordinates od the point A [yA, xA]
368      ptB= coordinates od the point B [yB, xB]
369      radii= radi of the circle to use to unite the points
370      Nang= amount of angles to use
371    """
372    fname = 'circ_sec'
373
374    distAB = dist_points(ptA,ptB)
375
376    if distAB > radii:
377        print errormsg
378        print '  ' + fname + ': radii=', radii, " too small for the distance " +     \
379          "between points !!"
380        print '    distance between points:', distAB
381        quit(-1)
382
383    # Coordinate increments
384    dAB = np.abs(ptA-ptB)
385
386    # angle of the circular section joining points
387    alpha = 2.*np.arcsin((distAB/2.)/radii)
388
389    # center along coincident bisection of the union
390    xcc = -radii
391    ycc = 0.
392
393    # Getting the arc of the circle at the x-axis
394    dalpha = alpha/(Nang-1)
395    circ_sec = np.zeros((Nang,2), dtype=np.float)
396    for ia in range(Nang):
397        alpha = dalpha*ia
398        x = radii*np.cos(alpha)
399        y = radii*np.sin(alpha)
400
401        circ_sec[ia,:] = [y+ycc,x+xcc]
402   
403    # Angle of the points
404    theta = np.arctan2(ptB[0]-ptA[0],ptB[1]-ptA[1])
405
406    # rotating angle of the circ
407    rotangle = theta + 3.*np.pi/2. - alpha/2.
408
409    #print 'alpha:', alpha*180./np.pi, 'theta:', theta*180./np.pi, 'rotangle:', rotangle*180./np.pi
410 
411
412    # rotating the arc along the x-axis
413    rotcirc_sec = rotate_polygon_2D(circ_sec, rotangle)
414
415    # Moving arc to the ptA
416    circ_sec = rotcirc_sec + ptA
417
418    return circ_sec
419
420def p_square(face, N=5):
421    """ Function to get a polygon square
422      face: length of the face of the square
423      N: number of points of the polygon
424    """
425    fname = 'p_square'
426
427    square = np.zeros((N,2), dtype=np.float)
428
429    f2 = face/2.
430    N4 = N/4
431    df = face/(N4)
432    # SW-NW
433    for ip in range(N4):
434        square[ip,:] = [-f2+ip*df,-f2]
435    # NW-NE
436    for ip in range(N4):
437        square[ip+N4,:] = [f2,-f2+ip*df]
438    # NE-SE
439    for ip in range(N4):
440        square[ip+2*N4,:] = [f2-ip*df,f2]
441    N42 = N-3*N4-1
442    df = face/(N42)
443    # SE-SW
444    for ip in range(N42):
445        square[ip+3*N4,:] = [-f2,f2-ip*df]
446    square[N-1,:] = [-f2,-f2]
447
448    return square
449
450def p_circle(radii, N=50):
451    """ Function to get a polygon of a circle
452      radii: length of the radii of the circle
453      N: number of points of the polygon
454    """
455    fname = 'p_circle'
456
457    circle = np.zeros((N,2), dtype=np.float)
458
459    dangle = 2.*np.pi/(N-1)
460
461    for ia in range(N):
462        circle[ia,:] = [radii*np.sin(ia*dangle), radii*np.cos(ia*dangle)]
463
464    circle[N-1,:] = [0., radii]
465
466    return circle
467
468def p_triangle(p1, p2, p3, N=4):
469    """ Function to provide the polygon of a triangle from its 3 vertices
470      p1: vertex 1 [y,x]
471      p2: vertex 2 [y,x]
472      p3: vertex 3 [y,x]
473      N: number of vertices of the triangle
474    """
475    fname = 'p_triangle'
476
477    triangle = np.zeros((N,2), dtype=np.float)
478
479    N3 = N / 3
480    # 1-2
481    dx = (p2[1]-p1[1])/N3
482    dy = (p2[0]-p1[0])/N3
483    for ip in range(N3):
484        triangle[ip,:] = [p1[0]+ip*dy,p1[1]+ip*dx]
485    # 2-3
486    dx = (p3[1]-p2[1])/N3
487    dy = (p3[0]-p2[0])/N3
488    for ip in range(N3):
489        triangle[ip+N3,:] = [p2[0]+ip*dy,p2[1]+ip*dx]
490    # 3-1
491    N32 = N - 2*N/3
492    dx = (p1[1]-p3[1])/N32
493    dy = (p1[0]-p3[0])/N32
494    for ip in range(N32):
495        triangle[ip+2*N3,:] = [p3[0]+ip*dy,p3[1]+ip*dx]
496
497    triangle[N-1,:] = p1
498
499    return triangle
500
501def p_spiral(loops, eradii, N=1000):
502    """ Function to provide a polygon of an Archimedean spiral
503        FROM: https://en.wikipedia.org/wiki/Spiral
504      loops: number of loops of the spiral
505      eradii: length of the radii of the final spiral
506      N: number of points of the polygon
507    """
508    fname = 'p_spiral'
509
510    spiral = np.zeros((N,2), dtype=np.float)
511
512    dangle = 2.*np.pi*loops/(N-1)
513    dr = eradii*1./(N-1)
514
515    for ia in range(N):
516        radii = dr*ia
517        spiral[ia,:] = [radii*np.sin(ia*dangle), radii*np.cos(ia*dangle)]
518
519    return spiral
520
521def p_reg_polygon(Nv, lf, N=50):
522    """ Function to provide a regular polygon of Nv vertices
523      Nv: number of vertices
524      lf: length of the face
525      N: number of points
526    """
527    fname = 'p_reg_polygon'
528
529    reg_polygon = np.zeros((N,2), dtype=np.float)
530
531    # Number of points per vertex
532    Np = N/Nv
533    # Angle incremental between vertices
534    da = 2.*np.pi/Nv
535    # Radii of the circle according to lf
536    radii = lf*Nv/(2*np.pi)
537
538    iip = 0
539    for iv in range(Nv-1):
540        # Characteristics between vertices iv and iv+1
541        av1 = da*iv
542        v1 = [radii*np.sin(av1), radii*np.cos(av1)]
543        av2 = da*(iv+1)
544        v2 = [radii*np.sin(av2), radii*np.cos(av2)]
545        dx = (v2[1]-v1[1])/Np
546        dy = (v2[0]-v1[0])/Np
547        for ip in range(Np):
548            reg_polygon[ip+iv*Np,:] = [v1[0]+dy*ip,v1[1]+dx*ip]
549
550    # Characteristics between vertices Nv and 1
551
552    # Number of points per vertex
553    Np2 = N - Np*(Nv-1)
554
555    av1 = da*Nv
556    v1 = [radii*np.sin(av1), radii*np.cos(av1)]
557    av2 = 0.
558    v2 = [radii*np.sin(av2), radii*np.cos(av2)]
559    dx = (v2[1]-v1[1])/Np2
560    dy = (v2[0]-v1[0])/Np2
561    for ip in range(Np2):
562        reg_polygon[ip+(Nv-1)*Np,:] = [v1[0]+dy*ip,v1[1]+dx*ip]
563
564    return reg_polygon
565
566def p_reg_star(Nv, lf, freq, vs=0, N=50):
567    """ Function to provide a regular star of Nv vertices
568      Nv: number of vertices
569      lf: length of the face of the regular polygon
570      freq: frequency of union of vertices ('0', for just centered to zero arms)
571      vs: vertex from which start (0 being first [0,lf])
572      N: number of points
573    """
574    fname = 'p_reg_star'
575
576    reg_star = np.zeros((N,2), dtype=np.float)
577
578    # Number of arms of the star
579    if freq != 0 and np.mod(Nv,freq) == 0: 
580        Na = Nv/freq + 1
581    else:
582        Na = Nv
583
584    # Number of points per arm
585    Np = N/Na
586    # Angle incremental between vertices
587    da = 2.*np.pi/Nv
588    # Radii of the circle according to lf
589    radii = lf*Nv/(2*np.pi)
590
591    iip = 0
592    av1 = vs*da
593    for iv in range(Na-1):
594        # Characteristics between vertices iv and iv+1
595        v1 = [radii*np.sin(av1), radii*np.cos(av1)]
596        if freq != 0:
597            av2 = av1 + da*freq
598            v2 = [radii*np.sin(av2), radii*np.cos(av2)]
599        else:
600            v2 = [0., 0.]
601            av2 = av1 + da
602        dx = (v2[1]-v1[1])/(Np-1)
603        dy = (v2[0]-v1[0])/(Np-1)
604        for ip in range(Np):
605            reg_star[ip+iv*Np,:] = [v1[0]+dy*ip,v1[1]+dx*ip]
606        if av2 > 2.*np.pi: av1 = av2 - 2.*np.pi
607        else: av1 = av2 + 0.
608
609    iv = Na-1
610    # Characteristics between vertices Na and 1
611    Np2 = N-Np*iv
612    v1 = [radii*np.sin(av1), radii*np.cos(av1)]
613    if freq != 0:
614        av2 = vs*da
615        v2 = [radii*np.sin(av2), radii*np.cos(av2)]
616    else:
617        v2 = [0., 0.]
618    dx = (v2[1]-v1[1])/(Np2-1)
619    dy = (v2[0]-v1[0])/(Np2-1)
620    for ip in range(Np2):
621        reg_star[ip+iv*Np,:] = [v1[0]+dy*ip,v1[1]+dx*ip]
622
623    return reg_star
624
625def p_sinusiode(length=10., amp=5., lamb=3., ival=0., func='sin', N=100):
626    """ Function to get coordinates of a sinusoidal curve
627      length: length of the line (default 10.)
628      amp: amplitude of the peaks (default 5.)
629      lamb: wave longitude (defalult 3.)
630      ival: initial angle (default 0. in degree)
631      func: function to use: (default sinus)
632        'sin': sinus
633        'cos': cosinus
634      N: number of points (default 100)
635    """
636    fname = 'p_sinusiode'
637    availfunc = ['sin', 'cos']
638
639    dx = length/(N-1)
640    ia = ival*np.pi/180.
641    da = 2*np.pi*dx/lamb
642
643    sinusoide = np.zeros((N,2), dtype=np.float)
644    if func == 'sin':
645        for ix in range(N):
646            sinusoide[ix,:] = [amp*np.sin(ia+da*ix),dx*ix]
647    elif func == 'cos':
648        for ix in range(N):
649            sinusoide[ix,:] = [amp*np.cos(ia+da*ix),dx*ix]
650    else:
651        print errormsg
652        print '  ' + fname + ": function '" + func + "' not ready !!"
653        print '    available ones:', availfunc
654        quit(-1)
655
656    sinusoidesecs = ['sinusoide']
657    sinusoidedic = {'sinusoide': [sinusoide, '-', '#000000', 1.]}
658
659    return sinusoide, sinusoidesecs, sinusoidedic
660
661def p_doubleArrow(length=5., angle=45., width=1., alength=0.10, N=50):
662    """ Function to provide an arrow with double lines
663      length: length of the arrow (5. default)
664      angle: angle of the head of the arrow (45., default)
665      width: separation between the two lines (2., default)
666      alength: length of the head (as percentage in excess of width, 0.1 default)
667      N: number of points (50, default)
668    """
669    import numpy.ma as ma
670    function = 'p_doubleArrow'
671
672    doubleArrow = np.zeros((50,2), dtype=np.float)
673    N4 = int((N-3)/4)
674
675    doublearrowdic = {}
676
677    # Arms
678    dx = length/(N4-1)
679    for ix in range(N4-1):
680        doubleArrow[ix,:] = [dx*ix,-width/2.]
681    doublearrowdic['leftarm'] = [doubleArrow[0:N4-1,:], '-', '#000000', 2.]
682    doubleArrow[N4-1,:] = [gen.fillValueF,gen.fillValueF]
683    for ix in range(N4-1):
684        doubleArrow[N4+ix,:] = [dx*ix,width/2.]
685    doublearrowdic['rightarm'] = [doubleArrow[N4:2*N4-1,:], '-', '#000000', 2.]
686    doubleArrow[2*N4-1,:] = [gen.fillValueF,gen.fillValueF]
687
688    # Head
689    N42 = int((N-2 - 2*N4)/2)
690    dx = width*(1.+alength)/(N42-1)
691    for ix in range(N42):
692        doubleArrow[2*N4+ix,:] = [length-dx*ix,-dx*ix]
693    doublearrowdic['lefthead'] = [doubleArrow[2*N4:2*N4+N42,:], '-', '#000000', 2.]
694    doubleArrow[2*N4+N42,:] = [gen.fillValueF,gen.fillValueF]
695
696    N43 = N-2 - 2*N4 - N42 + 1
697    dx = width*(1.+alength)/(N43-1)
698    for ix in range(N43):
699        doubleArrow[2*N4+N42+1+ix,:] = [length-dx*ix,dx*ix]
700    doublearrowdic['rightthead'] = [doubleArrow[2*N4+N42:51,:], '-', '#000000', 2.]
701
702    doubleArrow = ma.masked_equal(doubleArrow, gen.fillValueF)
703    doublearrowsecs = ['leftarm', 'rightarm', 'lefthead', 'righthead']
704
705    return doubleArrow, doublearrowsecs, doublearrowdic
706
707# Combined objects
708##
709
710# FROM: http://www.photographers1.com/Sailing/NauticalTerms&Nomenclature.html
711def zboat(length=10., beam=1., lbeam=0.4, sternbp=0.5):
712    """ Function to define an schematic boat from the z-plane
713      length: length of the boat (without stern, default 10)
714      beam: beam of the boat (default 1)
715      lbeam: length at beam (as percentage of length, default 0.4)
716      sternbp: beam at stern (as percentage of beam, default 0.5)
717    """
718    fname = 'zboat'
719
720    bow = np.array([length, 0.])
721    maxportside = np.array([length*lbeam, -beam])
722    maxstarboardside = np.array([length*lbeam, beam])
723    portside = np.array([0., -beam*sternbp])
724    starboardside = np.array([0., beam*sternbp])
725
726    # forward section
727    fportside = circ_sec(bow,maxportside, length*2)
728    fstarboardside = circ_sec(maxstarboardside, bow, length*2)
729    # aft section
730    aportside = circ_sec(maxportside, portside, length*2)
731    astarboardside = circ_sec(starboardside, maxstarboardside, length*2)
732    # stern
733    stern = circ_sec(portside, starboardside, length*2)
734
735    dpts = stern.shape[0]
736    boat = np.zeros((dpts*5,2), dtype=np.float)
737
738    boat[0:dpts,:] = fportside
739    boat[dpts:2*dpts,:] = aportside
740    boat[2*dpts:3*dpts,:] = stern
741    boat[3*dpts:4*dpts,:] = astarboardside
742    boat[4*dpts:5*dpts,:] = fstarboardside
743
744    fname = 'boat_L' + str(int(length*100.)) + '_B' + str(int(beam*100.)) + '_lb' +  \
745      str(int(lbeam*100.)) + '_sb' + str(int(sternbp*100.)) + '.dat'
746    if not os.path.isfile(fname):
747        print infmsg
748        print '  ' + fname + ": writting boat coordinates file '" + fname + "' !!"
749        of = open(fname, 'w')
750        of.write('# boat file with Length: ' + str(length) +' max_beam: '+str(beam)+ \
751          'length_at_max_beam:' + str(lbeam) + '% beam at stern: ' + str(sternbp)+   \
752          ' %\n')
753        for ip in range(dpts*5):
754            of.write(str(boat[ip,0]) + ' ' + str(boat[ip,1]) + '\n')
755       
756        of.close()
757        print fname + ": Successfull written '" + fname + "' !!"
758 
759    return boat
760
761def zsailing_boat(length=10., beam=1., lbeam=0.4, sternbp=0.5, lmast=0.6, wmast=0.1, \
762  hsd=5., msd=5., lheads=0.38, lmains=0.55):
763    """ Function to define an schematic sailing boat from the z-plane with sails
764      length: length of the boat (without stern, default 10)
765      beam: beam of the boat (default 1)
766      lbeam: length at beam (as percentage of length, default 0.4)
767      sternbp: beam at stern (as percentage of beam, default 0.5)
768      lmast: position of the mast (as percentage of length, default 0.6)
769      wmast: width of the mast (default 0.1)
770      hsd: head sail direction respect to center line (default 5., -999.99 for upwind)
771      msd: main sail direction respect to center line (default 5., -999.99 for upwind)
772      lheads: length of head sail (as percentage of legnth, defaul 0.38)
773      lmains: length of main sail (as percentage of legnth, defaul 0.55)
774    """
775    import numpy.ma as ma
776    fname = 'zsailing_boat'
777
778    bow = np.array([length, 0.])
779    maxportside = np.array([length*lbeam, -beam])
780    maxstarboardside = np.array([length*lbeam, beam])
781    portside = np.array([0., -beam*sternbp])
782    starboardside = np.array([0., beam*sternbp])
783
784    # forward section
785    fportside = circ_sec(bow,maxportside, length*2)
786    fstarboardside = circ_sec(maxstarboardside, bow, length*2)
787    dpts = fportside.shape[0]
788
789    # aft section
790    aportside = circ_sec(maxportside, portside, length*2)
791    astarboardside = circ_sec(starboardside, maxstarboardside, length*2)
792    # stern
793    stern = circ_sec(portside, starboardside, length*2)
794    # mast
795    mast = p_circle(wmast,N=dpts)
796    mast = mast + [length*lmast, 0.]
797    # head sails
798    lsail = lheads*length
799    if hsd != -999.99:
800        sailsa = np.pi/2. - np.pi*hsd/180.
801        endsail = np.array([lsail*np.sin(sailsa), lsail*np.cos(sailsa)])
802        endsail[0] = length - endsail[0]
803        if bow[1] < endsail[1]:
804            hsail = circ_sec(endsail, bow, lsail*2.15)
805        else:
806            hsail = circ_sec(bow, endsail, lsail*2.15)
807    else:
808        hsail0 = p_sinusiode(length=lsail, amp=0.2, lamb=0.75, N=dpts)
809        hsail = np.zeros((dpts,2), dtype=np.float)
810        hsail[:,0] = hsail0[:,1]
811        hsail[:,1] = hsail0[:,0]
812        hsail = bow - hsail
813
814    # main sails
815    lsail = lmains*length
816    if msd != -999.99:
817        sailsa = np.pi/2. - np.pi*msd/180.
818        begsail = np.array([length*lmast, 0.])
819        endsail = np.array([lsail*np.sin(sailsa), lsail*np.cos(sailsa)])
820        endsail[0] = length*lmast - endsail[0]
821        if endsail[1] < begsail[1]:
822            msail = circ_sec(begsail, endsail, lsail*2.15)
823        else:
824            msail = circ_sec(endsail, begsail, lsail*2.15)
825    else:
826        msail0 = p_sinusiode(length=lsail, amp=0.25, lamb=1., N=dpts)
827        msail = np.zeros((dpts,2), dtype=np.float)
828        msail[:,0] = msail0[:,1]
829        msail[:,1] = msail0[:,0]
830        msail = [length*lmast,0] - msail
831
832    sailingboat = np.zeros((dpts*8+4,2), dtype=np.float)
833
834    sailingboat[0:dpts,:] = fportside
835    sailingboat[dpts:2*dpts,:] = aportside
836    sailingboat[2*dpts:3*dpts,:] = stern
837    sailingboat[3*dpts:4*dpts,:] = astarboardside
838    sailingboat[4*dpts:5*dpts,:] = fstarboardside
839    sailingboat[5*dpts,:] = [gen.fillValueF, gen.fillValueF]
840    sailingboat[5*dpts+1:6*dpts+1,:] = mast
841    sailingboat[6*dpts+1,:] = [gen.fillValueF, gen.fillValueF]
842    sailingboat[6*dpts+2:7*dpts+2,:] = hsail
843    sailingboat[7*dpts+2,:] = [gen.fillValueF, gen.fillValueF]
844    sailingboat[7*dpts+3:8*dpts+3,:] = msail
845    sailingboat[8*dpts+3,:] = [gen.fillValueF, gen.fillValueF]
846
847    sailingboat = ma.masked_equal(sailingboat, gen.fillValueF)
848
849    # Center line extending [fcl] percentage from length on aft and stern
850    fcl = 0.15
851    centerline = np.zeros((dpts,2), dtype=np.float)
852    dl = length*(1.+fcl*2.)/(dpts-1)
853    centerline[:,0] = np.arange(-length*fcl, length*(1. + fcl)+dl, dl)
854
855    # correct order of sections
856    sailingboatsecs = ['fportside', 'aportside', 'stern', 'astarboardside',          \
857      'fstarboardside', 'mast', 'hsail', 'msail', 'centerline']
858    # dictionary with sections [polygon_vertices, line_type, line_color, line_width]
859    dicsailingboat = {'fportside': [fportside, '-', '#8A5900', 2.],                  \
860      'aportside': [aportside, '-', '#8A5900', 2.],                                  \
861      'stern': [stern, '-', '#8A5900', 2.],                                          \
862      'astarboardside': [astarboardside, '-', '#8A5900', 2.],                        \
863      'fstarboardside': [fstarboardside, '-', '#8A5900', 2.],                        \
864      'mast': [mast, '-', '#8A5900', 2.], 'hsail': [hsail, '-', '#AAAAAA', 1.],      \
865      'msail': [msail, '-', '#AAAAAA', 1.],                                          \
866      'centerline': [centerline, '-.', '#AA6464', 1.5]}
867   
868    fname = 'sailboat_L' + str(int(length*100.)) + '_B' + str(int(beam*100.)) +      \
869      '_lb' + str(int(lbeam*100.)) + '_sb' + str(int(sternbp*100.)) +                \
870      '_lm' + str(int(lmast*100.)) + '_wm' + str(int(wmast)) +                       \
871      '_hsd' + str(int(hsd)) + '_hs' + str(int(lheads*100.)) +                       \
872      '_ms' + str(int(lheads*100.)) + '_msd' + str(int(msd)) +'.dat'
873    if not os.path.isfile(fname):
874        print infmsg
875        print '  ' + fname + ": writting boat coordinates file '" + fname + "' !!"
876        of = open(fname, 'w')
877        of.write('# boat file with Length: ' + str(length) +' max_beam: '+str(beam)+ \
878          'length_at_max_beam:' + str(lbeam) + '% beam at stern: ' + str(sternbp)+   \
879          ' % mast position: '+ str(lmast) + ' % mast width: ' + str(wmast) + ' ' +  \
880          ' head sail direction:' + str(hsd) + ' head sail length: ' + str(lheads) + \
881          ' %' + ' main sail length' + str(lmains) + ' main sail direction:' +       \
882          str(msd) +'\n')
883        for ip in range(dpts*5):
884            of.write(str(sailingboat[ip,0]) + ' ' + str(sailingboat[ip,1]) + '\n')
885       
886        of.close()
887        print fname + ": Successfull written '" + fname + "' !!"
888 
889    return sailingboat, sailingboatsecs, dicsailingboat
890
891def write_join_poly(polys, flname='join_polygons.dat'):
892    """ Function to write an ASCII file with the combination of polygons
893      polys: dictionary with the names of the different polygons
894      flname: name of the ASCII file
895    """
896    fname = 'write_join_poly'
897
898    of = open(flname, 'w')
899
900    for polyn in polys.keys():
901        vertices = polys[polyn]
902        Npts = vertices.shape[0]
903        for ip in range(Npts):
904            of.write(polyn+' '+str(vertices[ip,1]) + ' ' + str(vertices[ip,0]) + '\n')
905
906    of.close()
907
908    return
909
910def read_join_poly(flname='join_polygons.dat'):
911    """ Function to read an ASCII file with the combination of polygons
912      flname: name of the ASCII file
913    """
914    fname = 'read_join_poly'
915
916    of = open(flname, 'r')
917
918    polys = {}
919    polyn = ''
920    poly = []
921    for line in of:
922        if len(line) > 1: 
923            linevals = line.replace('\n','').split(' ')
924            if polyn != linevals[0]:
925                if len(poly) > 1:
926                    polys[polyn] = np.array(poly)
927                polyn = linevals[0]
928                poly = []
929                poly.append([np.float(linevals[2]), np.float(linevals[1])])
930            else:
931                poly.append([np.float(linevals[2]), np.float(linevals[1])])
932
933    of.close()
934    polys[polyn] = np.array(poly)
935
936    return polys
937
938####### ####### ##### #### ### ## #
939# Plotting
940
941def plot_sphere(iazm=-60., iele=30., dist=10., Npts=100, radii=10,                   \
942  drwsfc=[True,True], colsfc=['#AAAAAA','#646464'],                                  \
943  drwxline = True, linex=[':','b',2.], drwyline = True, liney=[':','r',2.],          \
944  drwzline = True, linez=['-.','g',2.], drwxcline=[True,True],                       \
945  linexc=[['-','#646400',1.],['--','#646400',1.]],                                   \
946  drwequator=[True,True], lineeq=[['-','#AA00AA',1.],['--','#AA00AA',1.]],           \
947  drwgreeenwhich=[True,True], linegw=[['-','k',1.],['--','k',1.]]):
948    """ Function to plot an sphere and determine which standard lines will be also
949        drawn
950      iazm: azimut of the camera form the sphere
951      iele: elevation of the camera form the sphere
952      dist: distance of the camera form the sphere
953      Npts: Resolution for the sphere
954      radii: radius of the sphere
955      drwsfc: whether 'up' and 'down' portions of the sphere should be drawn
956      colsfc: colors of the surface of the sphere portions ['up', 'down']
957      drwxline: whether x-axis line should be drawn
958      linex: properties of the x-axis line ['type', 'color', 'wdith']
959      drwyline: whether y-axis line should be drawn
960      liney: properties of the y-axis line ['type', 'color', 'wdith']
961      drwzline: whether z-axis line should be drawn
962      linez: properties of the z-axis line ['type', 'color', 'wdith']
963      drwequator: whether 'front' and 'back' portions of the Equator should be drawn
964      lineeq: properties of the lines 'front' and 'back' of the Equator
965      drwgreeenwhich: whether 'front', 'back' portions of Greenqhich should be drawn
966      linegw: properties of the lines 'front' and 'back' Greenwhich
967      drwxcline: whether 'front', 'back' 90 line (lon=90., lon=270.) should be drawn
968      linexc: properties of the lines 'front' and 'back' for the 90 line
969    """
970    fname = 'plot_sphere'
971
972    iazmrad = iazm*np.pi/180.
973    ielerad = iele*np.pi/180.
974
975    # 3D surface Sphere
976    sfcsphereu, sfcsphered = surface_sphere(radii,Npts)
977   
978    # greenwhich
979    if iazmrad > np.pi/2. and iazmrad < 3.*np.pi/2.:
980        ia=np.pi-ielerad
981    else:
982        ia=0.-ielerad
983    ea=ia+np.pi
984    da = (ea-ia)/(Npts-1)
985    beta = np.arange(ia,ea+da,da)[0:Npts]
986    alpha = np.zeros((Npts), dtype=np.float)
987    greenwhichc = spheric_line(radii,alpha,beta)
988    ia=ea+0.
989    ea=ia+np.pi
990    da = (ea-ia)/(Npts-1)
991    beta = np.arange(ia,ea+da,da)[0:Npts]
992    greenwhichd = spheric_line(radii,alpha,beta)
993
994    # Equator
995    ia=np.pi-iazmrad/2.
996    ea=ia+np.pi
997    da = (ea-ia)/(Npts-1)
998    alpha = np.arange(ia,ea+da,da)[0:Npts]
999    beta = np.zeros((Npts), dtype=np.float)
1000    equatorc = spheric_line(radii,alpha,beta)
1001    ia=ea+0.
1002    ea=ia+np.pi
1003    da = (ea-ia)/(Npts-1)
1004    alpha = np.arange(ia,ea+da,da)[0:Npts]
1005    equatord = spheric_line(radii,alpha,beta)
1006
1007    # 90 line
1008    if iazmrad > np.pi and iazmrad < 2.*np.pi:
1009        ia=3.*np.pi/2. + ielerad
1010    else:
1011        ia=np.pi/2. - ielerad
1012    if ielerad < 0.:
1013        ia = ia + np.pi
1014    ea=ia+np.pi
1015    da = (ea-ia)/(Npts-1)
1016    beta = np.arange(ia,ea+da,da)[0:Npts]
1017    alpha = np.ones((Npts), dtype=np.float)*np.pi/2.
1018    xclinec = spheric_line(radii,alpha,beta)
1019    ia=ea+0.
1020    ea=ia+np.pi
1021    da = (ea-ia)/(Npts-1)
1022    beta = np.arange(ia,ea+da,da)[0:Npts]
1023    xclined = spheric_line(radii,alpha,beta)
1024
1025    # x line
1026    xline = np.zeros((2,3), dtype=np.float)
1027    xline[0,:] = position_sphere(radii, 0., 0.)
1028    xline[1,:] = position_sphere(radii, np.pi, 0.)
1029
1030    # y line
1031    yline = np.zeros((2,3), dtype=np.float)
1032    yline[0,:] = position_sphere(radii, np.pi/2., 0.)
1033    yline[1,:] = position_sphere(radii, 3*np.pi/2., 0.)
1034
1035    # z line
1036    zline = np.zeros((2,3), dtype=np.float)
1037    zline[0,:] = position_sphere(radii, 0., np.pi/2.)
1038    zline[1,:] = position_sphere(radii, 0., -np.pi/2.)
1039
1040    fig = plt.figure()
1041    ax = fig.gca(projection='3d')
1042
1043    # Sphere surface
1044    if drwsfc[0]:
1045        ax.plot_surface(sfcsphereu[0,:,:], sfcsphereu[1,:,:], sfcsphereu[2,:,:],     \
1046          color=colsfc[0])
1047    if drwsfc[1]:
1048        ax.plot_surface(sfcsphered[0,:,:], sfcsphered[1,:,:], sfcsphered[2,:,:],     \
1049          color=colsfc[1])
1050
1051    # greenwhich
1052    linev = linegw[0]
1053    if drwgreeenwhich[0]:
1054        ax.plot(greenwhichc[:,0], greenwhichc[:,1], greenwhichc[:,2], linev[0],      \
1055          color=linev[1], linewidth=linev[2],  label='Greenwich')
1056    linev = linegw[1]
1057    if drwgreeenwhich[1]:
1058        ax.plot(greenwhichd[:,0], greenwhichd[:,1], greenwhichd[:,2], linev[0],      \
1059          color=linev[1], linewidth=linev[2])
1060
1061    # Equator
1062    linev = lineeq[0]
1063    if drwequator[0]:
1064        ax.plot(equatorc[:,0], equatorc[:,1], equatorc[:,2], linev[0],               \
1065          color=linev[1], linewidth=linev[2], label='Equator')
1066    linev = lineeq[1]
1067    if drwequator[1]:
1068        ax.plot(equatord[:,0], equatord[:,1], equatord[:,2], linev[0],               \
1069          color=linev[1], linewidth=linev[2])
1070
1071    # 90line
1072    linev = linexc[0]
1073    if drwxcline[0]:
1074        ax.plot(xclinec[:,0], xclinec[:,1], xclinec[:,2], linev[0], color=linev[1],  \
1075          linewidth=linev[2], label='90-line')
1076    linev = linexc[1]
1077    if drwxcline[1]:
1078        ax.plot(xclined[:,0], xclined[:,1], xclined[:,2], linev[0], color=linev[1],  \
1079          linewidth=linev[2])
1080
1081    # x line
1082    linev = linex
1083    if drwxline:
1084        ax.plot([xline[0,0],xline[1,0]], [xline[0,1],xline[1,1]],                    \
1085          [xline[0,2],xline[1,2]], linev[0], color=linev[1], linewidth=linev[2],  label='xline')
1086
1087    # y line
1088    linev = liney
1089    if drwyline:
1090        ax.plot([yline[0,0],yline[1,0]], [yline[0,1],yline[1,1]],                    \
1091          [yline[0,2],yline[1,2]], linev[0], color=linev[1], linewidth=linev[2],  label='yline')
1092
1093    # z line
1094    linev = linez
1095    if drwzline:
1096        ax.plot([zline[0,0],zline[1,0]], [zline[0,1],zline[1,1]],                    \
1097          [zline[0,2],zline[1,2]], linev[0], color=linev[1], linewidth=linev[2],  label='zline')
1098
1099    plt.legend()
1100
1101    return fig, ax
Note: See TracBrowser for help on using the repository browser.