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

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

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