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

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

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