# Python tools to manage netCDF files. # L. Fita, CIMA. March 2019 # More information at: http://www.xn--llusfb-5va.cat/python/PyNCplot # # pyNCplot and its component geometry_tools.py comes with ABSOLUTELY NO WARRANTY. # This work is licendes under a Creative Commons # Attribution-ShareAlike 4.0 International License (http://creativecommons.org/licenses/by-sa/4.0) # ## Script for geometry calculations and operations as well as definition of different ### standard objects and shapes import numpy as np import matplotlib as mpl from mpl_toolkits.mplot3d import Axes3D import matplotlib.pyplot as plt import os import generic_tools as gen import numpy.ma as ma import module_ForSci as fsci errormsg = 'ERROR -- error -- ERROR -- error' infmsg = 'INFORMATION -- information -- INFORMATION -- information' ####### Contents: # add_secpolygon_list: Function to add a range of points of a polygon into a list # angle_vectors2D: Angle between two vectors with sign # cut_between_[x/y]polygon: Function to cut a polygon between 2 given value of the [x/y]-axis # cut_[x/y]polygon: Function to cut a polygon from a given value of the [x/y]-axis # deg_deci: Function to pass from degrees [deg, minute, sec] to decimal angles [rad] # displace_objdic_2D: Function to displace 2D plain the vertices of all polygons of an object # dist_points: Function to provide the distance between two points # join_circ_sec: Function to join aa series of points by circular segments # join_circ_sec_rand: Function to join aa series of points by circular segments with # random perturbations # max_coords_poly: Function to provide the extremes of the coordinates of a polygon # mirror_polygon: Function to reflex a polygon for a given axis # mod_vec: Function to compute the module of a vector # position_sphere: Function to tranform fom a point in lon, lat deg coordinates to # cartesian coordinates over an sphere # read_join_poly: Function to read an ASCII file with the combination of polygons # rm_consecpt_polygon: Function to remove consecutive same point of a polygon # rotate_2D: Function to rotate a vector by a certain angle in the plain # rotate_objdic_2D: Function to rotate 2D plain the vertices of all polygons of an object # rotate_polygon_2D: Function to rotate 2D plain the vertices of a polygon # rotate_line2D: Function to rotate a line given by 2 pairs of x,y coordinates by a # certain angle in the plain # rotate_lines2D: Function to rotate multiple lines given by mulitple pars of x,y # coordinates by a certain angle in the plain # spheric_line: Function to transform a series of locations in lon, lat coordinates # to x,y,z over an 3D spaceFunction to provide coordinates of a line on a 3D space # val_consec_between: Function to provide if a given value is between two consecutive ones # write_join_poly: Function to write an ASCII file with the combination of polygons ## Shapes/objects # circ_sec: Function union of point A and B by a section of a circle # ellipse_polar: Function to determine an ellipse from its center and polar coordinates # p_angle_triangle: Function to draw a triangle by an initial point and two # consecutive angles and the first length of face. The third angle and 2 and 3rd # face will be computed accordingly the provided values # p_doubleArrow: Function to provide an arrow with double lines # p_circle: Function to get a polygon of a circle # p_cross_width: Function to draw a cross with arms with a given width and an angle # p_prism: Function to get a polygon prism # p_reg_polygon: Function to provide a regular polygon of Nv vertices # p_reg_star: Function to provide a regular star of Nv vertices # p_sinusoide: Function to get coordinates of a sinusoidal curve # p_square: Function to get a polygon square # p_spiral: Function to provide a polygon of an Archimedean spiral # p_triangle: Function to provide the polygon of a triangle from its 3 vertices # surface_sphere: Function to provide an sphere as matrix of x,y,z coordinates ## Plotting # draw_secs: Function to draw an object according to its dictionary # paint_filled: Function to draw an object filling given sections # plot_sphere: Function to plot an sphere and determine which standard lines will be # also drawn def deg_deci(angle): """ Function to pass from degrees [deg, minute, sec] to decimal angles [rad] angle: list of [deg, minute, sec] to pass >>> deg_deci([41., 58., 34.]) 0.732621346072 """ fname = 'deg_deci' deg = np.abs(angle[0]) + np.abs(angle[1])/60. + np.abs(angle[2])/3600. if angle[0] < 0.: deg = -deg*np.pi/180. else: deg = deg*np.pi/180. return deg def position_sphere(radii, alpha, beta): """ Function to tranform fom a point in lon, lat deg coordinates to cartesian coordinates over an sphere radii: radii of the sphere alpha: longitude of the point beta: latitude of the point >>> position_sphere(10., 30., 45.) (0.81031678432964027, -5.1903473778327376, 8.5090352453411846 """ fname = 'position_sphere' xpt = radii*np.cos(beta)*np.cos(alpha) ypt = radii*np.cos(beta)*np.sin(alpha) zpt = radii*np.sin(beta) return xpt, ypt, zpt def spheric_line(radii,lon,lat): """ Function to transform a series of locations in lon, lat coordinates to x,y,z over an 3D space radii: radius of the sphere lon: array of angles along longitudes lat: array of angles along latitudes """ fname = 'spheric_line' Lint = lon.shape[0] coords = np.zeros((Lint,3), dtype=np.float) for iv in range(Lint): coords[iv,:] = position_sphere(radii, lon[iv], lat[iv]) return coords def rotate_2D(vector, angle): """ Function to rotate a vector by a certain angle in the plain vector= vector to rotate [y, x] angle= angle to rotate [rad] >>> rotate_2D(np.array([1.,0.]), np.pi/4.) [ 0.70710678 -0.70710678] """ fname = 'rotate_2D' rotmat = np.zeros((2,2), dtype=np.float) rotmat[0,0] = np.cos(angle) rotmat[0,1] = -np.sin(angle) rotmat[1,0] = np.sin(angle) rotmat[1,1] = np.cos(angle) rotvector = np.zeros((2), dtype=np.float) vecv = np.zeros((2), dtype=np.float) # Unifying vector modvec = vector[0]**2+vector[1]**2 if modvec != 0 and vector[0] != gen.fillValue: vecv[0] = vector[1]/modvec vecv[1] = vector[0]/modvec rotvec = np.matmul(rotmat, vecv) rotvec = np.where(np.abs(rotvec) < 1.e-7, 0., rotvec) rotvector[0] = rotvec[1]*modvec rotvector[1] = rotvec[0]*modvec else: rotvector = vector + 0. return rotvector def rotate_polygon_2D(vectors, angle): """ Function to rotate 2D plain the vertices of a polygon line= matrix of vectors to rotate angle= angle to rotate [rad] >>> square = np.zeros((4,2), dtype=np.float) >>> square[0,:] = [-0.5,-0.5] >>> square[1,:] = [0.5,-0.5] >>> square[2,:] = [0.5,0.5] >>> square[3,:] = [-0.5,0.5] >>> rotate_polygon_2D(square, np.pi/4.) [[-0.70710678 0. ] [ 0. -0.70710678] [ 0.70710678 0. ] [ 0. 0.70710678]] """ fname = 'rotate_polygon_2D' rotvecs = np.zeros(vectors.shape, dtype=np.float) mavec = False if type(vectors) == type(gen.mamat): mavec = True vectors = ma.filled(vectors,gen.fillValueF) Nvecs = vectors.shape[0] for iv in range(Nvecs): rotvecs[iv,:] = rotate_2D(vectors[iv,:], angle) if mavec: rotvecs = ma.masked_equal(rotvecs, gen.fillValueF) return rotvecs def displace_objdic_2D(objdic, distance): """ Function to displace 2D plain the vertices of all polygons of an object objdic= dictionary with all the polygons of the object distance= distance to displace [ydist, xdist] """ fname = 'displace_objdic_2D' disobjdic = dict(objdic) for secn in objdic.keys(): objv = objdic[secn] vectors = objv[0] lt = objv[1] lc = objv[2] lw = objv[3] disvecs = np.zeros(vectors.shape, dtype=np.float) disvecs = vectors + distance disobjdic[secn] = [disvecs, lt, lc, lw] return disobjdic def rotate_objdic_2D(objdic, angle): """ Function to rotate 2D plain the vertices of all polygons of an object objdic= dictionary with all the polygons of the object angle= angle to rotate [rad] """ fname = 'rotate_objdic_2D' rotobjdic = dict(objdic) for secn in objdic.keys(): objv = objdic[secn] vectors = objv[0] lt = objv[1] lc = objv[2] lw = objv[3] rotvecs = np.zeros(vectors.shape, dtype=np.float) Nvecs = vectors.shape[0] for iv in range(Nvecs): rotvecs[iv,:] = rotate_2D(vectors[iv,:], angle) rotobjdic[secn] = [rotvecs, lt, lc, lw] return rotobjdic def rotate_line2D(line, angle): """ Function to rotate a line given by 2 pairs of x,y coordinates by a certain angle in the plain line= line to rotate as couple of points [[y0,x0], [y1,x1]] angle= angle to rotate [rad] >>> rotate_line2D(np.array([[0.,0.], [1.,0.]]), np.pi/4.) [[ 0. 0. ] [0.70710678 -0.70710678]] """ fname = 'rotate_2D' rotline = np.zeros((2,2), dtype=np.float) rotline[0,:] = rotate_2D(line[0,:], angle) rotline[1,:] = rotate_2D(line[1,:], angle) return rotline def rotate_lines2D(lines, angle): """ Function to rotate multiple lines given by mulitple pars of x,y coordinates by a certain angle in the plain line= matrix of N couples of points [N, [y0,x0], [y1,x1]] angle= angle to rotate [rad] >>> square = np.zeros((4,2,2), dtype=np.float) >>> square[0,0,:] = [-0.5,-0.5] >>> square[0,1,:] = [0.5,-0.5] >>> square[1,0,:] = [0.5,-0.5] >>> square[1,1,:] = [0.5,0.5] >>> square[2,0,:] = [0.5,0.5] >>> square[2,1,:] = [-0.5,0.5] >>> square[3,0,:] = [-0.5,0.5] >>> square[3,1,:] = [-0.5,-0.5] >>> rotate_lines2D(square, np.pi/4.) [[[-0.70710678 0. ] [ 0. -0.70710678]] [[ 0. -0.70710678] [ 0.70710678 0. ]] [[ 0.70710678 0. ] [ 0. 0.70710678]] [[ 0. 0.70710678] [-0.70710678 0. ]]] """ fname = 'rotate_lines2D' rotlines = np.zeros(lines.shape, dtype=np.float) Nlines = lines.shape[0] for il in range(Nlines): line = np.zeros((2,2), dtype=np.float) line[0,:] = lines[il,0,:] line[1,:] = lines[il,1,:] rotlines[il,:,:] = rotate_line2D(line, angle) return rotlines def dist_points(ptA, ptB): """ Function to provide the distance between two points ptA: coordinates of the point A [yA, xA] ptB: coordinates of the point B [yB, xB] >>> dist_points([1.,1.], [-1.,-1.]) 2.82842712475 """ fname = 'dist_points' dist = np.sqrt( (ptA[0]-ptB[0])**2 + (ptA[1]-ptB[1])**2) return dist def mod_vec(vec): """ Function to compute the module of a vector vec: vector [y, x] >>> mod_vec([1., 1.]) 1.41421356237 """ fname = 'mod_vec' v = np.array(vec, dtype=np.float) vv = v*v mod = np.sqrt(np.sum(vv[:])) return mod def angle_vectors2D(veca, vecb): """ Angle between two vectors with sign FROM: https://stackoverflow.com/questions/5188561/signed-angle-between-two-3d-vectors-with-same-origin-within-the-same-plane veca: angle A [ya, xa] vecb: angle B [yb, xb] NOTE: angle from A to B >>> angle_vectors2D([1.,0.], [0.,1.]) 1.57079632679 >>> angle_vectors2D([0.,1.], [1.,0.]) -1.57079632679 """ fname = 'angle_vectors2D' v1 = np.array(veca, dtype=np.float) v2 = np.array(vecb, dtype=np.float) moda = mod_vec(v1) modb = mod_vec(v2) modab = mod_vec(v1*v2) vc = np.cross(v1,v2) theta = np.arcsin(vc/(moda*modb)) # Without sign #alpha = np.arccos(modab/(moda*modb)) return theta def max_coords_poly(polygon): """ Function to provide the extremes of the coordinates of a polygon polygon: coordinates [Nvertexs, 2] of a polygon >>> square = np.zeros((4,2), dtype=np.float) >>> square[0,:] = [-0.5,-0.5] >>> square[1,:] = [0.5,-0.5] >>> square[2,:] = [0.5,0.5] >>> square[3,:] = [-0.5,0.5] >>> max_coords_poly(square) [-0.5, 0.5], [-0.5, 0.5], [0.5, 0.5], 0.5 """ fname = 'max_coords_poly' # x-coordinate min/max nx = np.min(polygon[:,1]) xx = np.max(polygon[:,1]) # y-coordinate min/max ny = np.min(polygon[:,0]) xy = np.max(polygon[:,0]) # x/y-coordinate maximum of absolute values axx = np.max(np.abs(polygon[:,1])) ayx = np.max(np.abs(polygon[:,0])) # absolute maximum xyx = np.max([axx, ayx]) return [nx, xx], [ny, xy], [ayx, axx], xyx def mirror_polygon(polygon,axis): """ Function to reflex a polygon for a given axis polygon: polygon to mirror axis: axis at which mirror is located ('x' or 'y') """ fname = 'mirror_polygon' reflex = np.zeros(polygon.shape, dtype=np.float) N = polygon.shape[0] if axis == 'x': for iv in range(N): reflex[iv,:] = [-polygon[iv,0], polygon[iv,1]] elif axis == 'y': for iv in range(N): reflex[iv,:] = [polygon[iv,0], -polygon[iv,1]] return reflex def join_circ_sec(points, radfrac=3., arc='short', side='left', N=200): """ Function to join aa series of points by circular segments points: main points of the island (clockwise ordered, to be joined by circular segments of radii as the radfrac factor of the distance between consecutive points) radfrac: multiplicative factor of the distance between consecutive points to draw the circular segment (3., default) arc: type of arc ('short', default) pos: position of arc ('left', default) N: number of points (200, default) """ fname = 'join_circ_sec' jcirc_sec = np.ones((N,2), dtype=np.float) # main points lpoints = list(points) Npts = len(lpoints) Np = int(N/(Npts+1)) for ip in range(Npts-1): p1 = lpoints[ip] p2 = lpoints[ip+1] dps = dist_points(p1, p2) jcirc_sec[Np*ip:Np*(ip+1),:] = circ_sec(p1,p2,dps*radfrac, arc, side, Np) Np2 = N - (Npts-1)*Np p1 = lpoints[Npts-1] p2 = lpoints[0] dps = dist_points(p1, p2) jcirc_sec[(Npts-1)*Np:N,:] = circ_sec(p1, p2, dps*3., arc, side, Np2) return jcirc_sec def join_circ_sec_rand(points, radfrac=3., Lrand=0.1, arc='short', pos='left', N=200): """ Function to join aa series of points by circular segments with random perturbations points: main points of the island (clockwise ordered, to be joined by circular segments of radii as the radfrac factor of the distance between consecutive points) radfrac: multiplicative factor of the distance between consecutive points to draw the circular segment (3., default) Lrand: maximum length of the random perturbation to be added perpendicularly to the direction of the union line between points (0.1, default) arc: type of arc ('short', default) pos: position of arc ('left', default) N: number of points (200, default) """ import random fname = 'join_circ_sec_rand' jcirc_sec = np.ones((N,2), dtype=np.float) # main points lpoints = list(points) Npts = len(lpoints) Np = int(N/(Npts+1)) for ip in range(Npts-1): p1 = lpoints[ip] p2 = lpoints[ip+1] dps = dist_points(p1, p2) angle = np.arctan2(p2[0]-p1[0], p2[1]-p1[1]) + np.pi/2. jcirc_sec[Np*ip:Np*(ip+1),:] = circ_sec(p1, p2, dps*radfrac, arc, pos, Np) drand = Lrand*np.array([np.sin(angle), np.cos(angle)]) for iip in range(Np*ip,Np*(ip+1)): jcirc_sec[iip,:] = jcirc_sec[iip,:] + drand*random.uniform(-1.,1.) Np2 = N - (Npts-1)*Np p1 = lpoints[Npts-1] p2 = lpoints[0] dps = dist_points(p1, p2) angle = np.arctan2(p2[0]-p1[0], p2[1]-p1[1]) + np.pi/2. jcirc_sec[(Npts-1)*Np:N,:] = circ_sec(p1, p2, dps*3., arc, pos, Np2) drand = Lrand*np.array([np.sin(angle), np.cos(angle)]) for iip in range(Np*(Npts-1),N): jcirc_sec[iip,:] = jcirc_sec[iip,:] + drand*random.uniform(-1.,1.) return jcirc_sec def write_join_poly(polys, flname='join_polygons.dat'): """ Function to write an ASCII file with the combination of polygons polys: dictionary with the names of the different polygons flname: name of the ASCII file """ fname = 'write_join_poly' of = open(flname, 'w') for polyn in polys.keys(): vertices = polys[polyn] Npts = vertices.shape[0] for ip in range(Npts): of.write(polyn+' '+str(vertices[ip,1]) + ' ' + str(vertices[ip,0]) + '\n') of.close() return def read_join_poly(flname='join_polygons.dat'): """ Function to read an ASCII file with the combination of polygons flname: name of the ASCII file """ fname = 'read_join_poly' of = open(flname, 'r') polys = {} polyn = '' poly = [] for line in of: if len(line) > 1: linevals = line.replace('\n','').split(' ') if polyn != linevals[0]: if len(poly) > 1: polys[polyn] = np.array(poly) polyn = linevals[0] poly = [] poly.append([np.float(linevals[2]), np.float(linevals[1])]) else: poly.append([np.float(linevals[2]), np.float(linevals[1])]) of.close() polys[polyn] = np.array(poly) return polys def val_consec_between(valA, valB, val): """ Function to provide if a given value is between two consecutive ones valA: first value valB: second value val: value to determine if it is between >>> val_consec_between(0.5,1.5,0.8) True >>> val_consec_between(0.5,1.5.,-0.8) False >>> val_consec_between(0.5,1.5,0.5) True >>> val_consec_between(-1.58, -1.4, -1.5) True >>> val_consec_between(-1.48747753212, -1.57383530044, -1.5) False """ fname = 'val_consec_between' btw = False diffA = valA - val diffB = valB - val absdA = np.abs(diffA) absdB = np.abs(diffB) #if (diffA/absdA)* (diffB/absdB) < 0.: btw = True # if valA < 0. and valB < 0. and val < 0.: # if (valA >= val and valB < val) or (valA > val and valB <= val): btw =True # else: # if (valA <= val and valB > val) or (valA < val and valB >= val): btw =True if (valA <= val and valB > val) or (valA < val and valB >= val): btw =True return btw def add_secpolygon_list(listv, iip, eep, polygon): """ Function to add a range of points of a polygon into a list listv: list into which add values of the polygon iip: initial value of the range eep: ending value of the range polygon: array with the points of the polygon """ fname = 'add_secpolygon_list' if eep > iip: for ip in range(iip,eep): listv.append(polygon[ip,:]) else: for ip in range(iip,eep,-1): listv.append(polygon[ip,:]) return def rm_consecpt_polygon(polygon): """ Function to remove consecutive same point of a polygon poly: polygon >>> poly = np.ones((5,2), dtype=np.float) >>> poly[2,:] = [2., 1.] rm_consecpt_polygon(poly) [[ 1. 1.] [ 2. 1.] [ 1. 1.]] """ fname = 'rm_consecpt_polygon' newpolygon = [] prevpt = polygon[0,:] newpolygon.append(prevpt) for ip in range(1,polygon.shape[0]): if polygon[ip,0] != prevpt[0] or polygon[ip,1] != prevpt[1]: prevpt = polygon[ip,:] newpolygon.append(prevpt) newpolygon = np.array(newpolygon) return newpolygon def cut_ypolygon(polygon, yval, keep='below', Nadd=20): """ Function to cut a polygon from a given value of the y-axis polygon: polygon to cut yval: value to use to cut the polygon keep: part to keep from the height ('below', default) 'below': below the height 'above': above the height Nadd: additional points to add to draw the line (20, default) """ fname = 'cut_ypolygon' N = polygon.shape[0] availkeeps = ['below', 'above'] if not gen.searchInlist(availkeeps, keep): print errormsg print ' ' + fname + ": wring keep '" + keep + "' value !!" print ' available ones:', availkeeps quit(-1) ipt = None ept = None # There might be more than 1 cut... Ncuts = 0 icut = [] ecut = [] ipt = [] ept = [] if type(polygon) == type(gen.mamat) and type(polygon.mask) != \ type(gen.mamat.mask[1]): # Assuming clockwise polygons for ip in range(N-1): if not polygon.mask[ip,0]: eep = ip + 1 if eep == N: eep = 0 if val_consec_between(polygon[ip,0], polygon[eep,0], yval): icut.append(ip) dx = polygon[eep,1] - polygon[ip,1] dy = polygon[eep,0] - polygon[ip,0] dd = yval - polygon[ip,0] ipt.append([yval, polygon[ip,1]+dx*dd/dy]) if val_consec_between(polygon[eep,0], polygon[ip,0], yval): ecut.append(ip) dx = polygon[eep,1] - polygon[ip,1] dy = polygon[eep,0] - polygon[ip,0] dd = yval - polygon[ip,0] ept.append([yval, polygon[ip,1]+dx*dd/dy]) Ncuts = Ncuts + 1 else: # Assuming clockwise polygons for ip in range(N-1): eep = ip + 1 if eep == N: eep = 0 if val_consec_between(polygon[ip,0], polygon[eep,0], yval): icut.append(ip) dx = polygon[eep,1] - polygon[ip,1] dy = polygon[eep,0] - polygon[ip,0] dd = yval - polygon[ip,0] ipt.append([yval, polygon[ip,1]+dx*dd/dy]) if val_consec_between(polygon[eep,0], polygon[ip,0], yval): ecut.append(ip) dx = polygon[eep,1] - polygon[ip,1] dy = polygon[eep,0] - polygon[ip,0] dd = yval - polygon[ip,0] ept.append([yval, polygon[ip,1]+dx*dd/dy]) Ncuts = Ncuts + 1 # Looking for repeated newicut = icut + [] newecut = ecut + [] newipt = ipt + [] newept = ept + [] newNcuts = Ncuts for ic in range(newNcuts-1): for ic2 in range(ic+1,newNcuts): if newipt[ic] == newipt[ic2]: Ncuts = Ncuts-1 icut.pop(ic2) ecut.pop(ic2) ipt.pop(ic2) ept.pop(ic2) newNcuts = Ncuts + 0 if ipt is None or ept is None or Ncuts == 0: print errormsg print ' ' + fname + ': no cutting for polygon at y=', yval, '!!' else: print ' ' + fname + ': found ', Ncuts, ' Ncuts' if Ncuts > 1 and keep == 'below': # Re-shifting cuts by closest distance. xis = [] xes = [] for ic in range(Ncuts): xp = ipt[ic] xis.append(xp[1]) xp = ept[ic] xes.append(xp[1]) xs = xis + xes xs.sort() newicut = icut + [] newecut = ecut + [] newipt = ipt + [] newept = ept + [] icut = [] ecut = [] ipt = [] ept = [] for xv in xs: ic = xis.count(xv) if ic != 0: icc = xis.index(xv) if len(icut) > len(ecut): ecut.append(newicut[icc]) ept.append(newipt[icc]) else: icut.append(newicut[icc]) ipt.append(newipt[icc]) else: icc = xes.index(xv) if len(icut) > len(ecut): ecut.append(newecut[icc]) ept.append(newept[icc]) else: icut.append(newecut[icc]) ipt.append(newept[icc]) # # Re-shifting cuts. 1st icut --> last ecut; 1st ecut as 1st icut; # # 2nd icut --> last-1 ecut, .... # newicut = icut + [] # newecut = ecut + [] # newipt = ipt + [] # newept = ept + [] # for ic in range(Ncuts-1): # ecut[ic] = newecut[Ncuts-ic-1] # ept[ic] = newept[Ncuts-ic-1] # icut[ic+1] = newecut[ic] # ipt[ic+1] = newept[ic] # ecut[Ncuts-1] = newicut[Ncuts-1] # ept[Ncuts-1] = newipt[Ncuts-1] ## print ' yval=', yval, 'cut, ip; ipt ep; ept ________' ## for ic in range(Ncuts): ## print ' ', ic, icut[ic], ';', ipt[ic], ecut[ic], ';', ept[ic] # Length of joining lines Nadds = [] if Ncuts > 1: Naddc = (Nadd-Ncuts)/(Ncuts) if Naddc < 3: print errormsg print ' ' + fname + ': too few points for jioning lines !!' print ' increase Nadd at least to:', Ncuts*3+Ncuts quit(-1) for ic in range(Ncuts-1): Nadds.append(Naddc) Nadds.append(Nadd-Naddc*(Ncuts-1)) else: Nadds.append(Nadd) # Total points cut polygon Ntotpts = 0 Ncpts = [] for ic in range(Ncuts): if keep == 'below': if ic == 0: dpts = icut[ic] + Nadds[ic] + (N - ecut[ic]) else: dpts = ecut[ic] - icut[ic] + Nadds[ic] - 1 # Adding end of the polygon in 'left' keeps if ic == Ncuts - 1: dpts = dpts + N-ecut[ic] else: dpts = ecut[ic] - icut[ic] + Nadds[ic] - 1 Ntotpts = Ntotpts + dpts Ncpts.append(ecut[ic] - icut[ic]) cutpolygon = np.ones((Ntotpts+Ncuts,2), dtype=np.float)*gen.fillValue iipc = 0 for ic in range(Ncuts): dcpt = Ncpts[ic] if keep == 'below': if ic == 0: cutpolygon[0:icut[ic],:] = polygon[0:icut[ic],:] iipc = icut[ic] else: cutpolygon[iipc:iipc+dcpt-1,:] = polygon[icut[ic]+1:ecut[ic],:] iipc = iipc + dcpt -1 else: cutpolygon[iipc,:] = ipt[ic] cutpolygon[iipc:iipc+dcpt-1,:]=polygon[icut[ic]+1:ecut[ic],:] iipc = iipc+dcpt-1 # cutting line cutline = np.zeros((Nadds[ic],2), dtype=np.float) dx = (ept[ic][1] - ipt[ic][1])/(Nadds[ic]-1) dy = (ept[ic][0] - ipt[ic][0])/(Nadds[ic]-1) cutline[0,:] = ipt[ic] for ip in range(1,Nadds[ic]-1): cutline[ip,:] = ipt[ic] + np.array([dy*ip,dx*ip]) cutline[Nadds[ic]-1,:] = ept[ic] if keep == 'below': if ic == 0: cutpolygon[iipc:iipc+Nadds[ic],:] = cutline else: cutpolygon[iipc:iipc+Nadds[ic],:] = cutline[::-1,:] iipc = iipc+Nadds[ic] if ic == 0: cutpolygon[iipc:iipc+N-ecut[ic]-1,:] = polygon[ecut[ic]+1:N,:] iipc = iipc + N-ecut[ic]-1 cutpolygon[iipc,:] = polygon[0,:] iipc = iipc + 1 else: cutpolygon[iipc:iipc+Nadds[ic],:] = cutline[::-1,:] iipc = iipc+Nadds[ic] iipc = iipc + 1 rmpolygon = [] Npts = cutpolygon.shape[0] if keep == 'below': for ip in range(Npts): if cutpolygon[ip,0] > yval: rmpolygon.append([gen.fillValueF, gen.fillValueF]) else: rmpolygon.append(cutpolygon[ip,:]) else: for ip in range(Npts): if cutpolygon[ip,0] < yval: rmpolygon.append([gen.fillValueF, gen.fillValueF]) else: rmpolygon.append(cutpolygon[ip,:]) Npts = len(rmpolygon) cutpolygon = np.array(rmpolygon) cutpolygon = rm_consecpt_polygon(cutpolygon) cutpolygon = ma.masked_equal(cutpolygon, gen.fillValueF) return Npts, cutpolygon def cut_xpolygon(polygon, xval, keep='left', Nadd=20): """ Function to cut a polygon from a given value of the x-axis polygon: polygon to cut yval: value to use to cut the polygon keep: part to keep from the value ('left', default) 'left': left of the value 'right': right of the value Nadd: additional points to add to draw the line (20, default) """ fname = 'cut_xpolygon' N = polygon.shape[0] availkeeps = ['left', 'right'] if not gen.searchInlist(availkeeps, keep): print errormsg print ' ' + fname + ": wring keep '" + keep + "' value !!" print ' available ones:', availkeeps quit(-1) ipt = None ept = None # There might be more than 1 cut ... icut = [] ecut = [] ipt = [] ept = [] Ncuts = 0 if type(polygon) == type(gen.mamat) and type(polygon.mask) != \ type(gen.mamat.mask[1]): # Assuming clockwise polygons for ip in range(N-1): if not polygon.mask[ip,1]: eep = ip + 1 if eep == N: eep = 0 if val_consec_between(polygon[ip,1], polygon[eep,1], xval): icut.append(ip) dx = polygon[eep,1] - polygon[ip,1] dy = polygon[eep,0] - polygon[ip,0] dd = xval - polygon[ip,1] ipt.append([polygon[ip,0]+dy*dd/dx, xval]) if val_consec_between(polygon[eep,1], polygon[ip,1], xval): ecut.append(ip) dx = polygon[eep,1] - polygon[ip,1] dy = polygon[eep,0] - polygon[ip,0] dd = xval - polygon[ip,1] ept.append([polygon[ip,0]+dy*dd/dx, xval]) Ncuts = Ncuts + 1 else: # Assuming clockwise polygons for ip in range(N-1): eep = ip + 1 if eep == N: eep = 0 if val_consec_between(polygon[ip,1], polygon[eep,1], xval): icut.append(ip) dx = polygon[eep,1] - polygon[ip,1] dy = polygon[eep,0] - polygon[ip,0] dd = xval - polygon[ip,1] ipt.append([polygon[ip,0]+dy*dd/dx, xval]) if val_consec_between(polygon[eep,1], polygon[ip,1], xval): ecut.append(ip) dx = polygon[eep,1] - polygon[ip,1] dy = polygon[eep,0] - polygon[ip,0] dd = xval - polygon[ip,1] ept.append([polygon[ip,0]+dy*dd/dx, xval]) Ncuts = Ncuts + 1 # Looking for repeated newicut = icut + [] newecut = ecut + [] newipt = ipt + [] newept = ept + [] newNcuts = Ncuts for ic in range(newNcuts-1): for ic2 in range(ic+1,newNcuts): if newipt[ic] == newipt[ic2]: Ncuts = Ncuts-1 icut.pop(ic2) ecut.pop(ic2) ipt.pop(ic2) ept.pop(ic2) newNcuts = Ncuts + 0 if ipt is None or ept is None or Ncuts == 0: print errormsg print ' ' + fname + ': no cutting for polygon at x=', xval, '!!' else: ##print ' ' + fname + ': found ', Ncuts, ' Ncuts' if Ncuts >= 1 and keep == 'left': # Re-shifting cuts by closest heigth. yis = [] yes = [] for ic in range(Ncuts): yp = ipt[ic] yis.append(yp[0]) yp = ept[ic] yes.append(yp[0]) ys = yis + yes ys.sort() newicut = icut + [] newecut = ecut + [] newipt = ipt + [] newept = ept + [] icut = [] ecut = [] ipt = [] ept = [] for yv in ys: ic = yis.count(yv) if ic != 0: icc = yis.index(yv) if len(icut) > len(ecut): ecut.append(newicut[icc]) ept.append(newipt[icc]) else: icut.append(newicut[icc]) ipt.append(newipt[icc]) else: icc = yes.index(yv) if len(icut) > len(ecut): ecut.append(newecut[icc]) ept.append(newept[icc]) else: icut.append(newecut[icc]) ipt.append(newept[icc]) #print ' xval=', xval, 'cut, ip; ipt ep; ept ________' #for ic in range(Ncuts): # print ' ', ic, icut[ic], ';', ipt[ic], ecut[ic], ';', ept[ic] # Length of joining lines Nadds = [] if Ncuts > 1: Naddc = (Nadd-Ncuts)/(Ncuts) if Naddc < 3: print errormsg print ' ' + fname + ': too few points for jioning lines !!' print ' increase Nadd at least to:', Ncuts*3+Ncuts quit(-1) for ic in range(Ncuts-1): Nadds.append(Naddc) Nadds.append(Nadd-Naddc*(Ncuts-1)) else: Nadds.append(Nadd) # Total points cut polygon Ntotpts = 0 Ncpts = [] for ic in range(Ncuts): if keep == 'left': if ic == 0: dpts = icut[ic] + Nadds[ic] + (N - ecut[ic]) else: dpts = ecut[ic] - icut[ic] + Nadds[ic] - 1 # Adding end of the polygon in 'left' keeps if ic == Ncuts - 1: dpts = dpts + N-ecut[ic] else: dpts = ecut[ic] - icut[ic] + Nadds[ic] - 1 Ntotpts = Ntotpts + dpts Ncpts.append(ecut[ic] - icut[ic]) cutpolygon = [] iipc = 0 for ic in range(Ncuts): dcpt = Ncpts[ic] cutpolygon.append(ipt[ic]) if keep == 'left': if ic == 0: add_secpolygon_list(cutpolygon,icut[ic]+1,N,polygon) add_secpolygon_list(cutpolygon,0,ecut[ic],polygon) iipc = icut[ic] else: add_secpolygon_list(cutpolygon,icut[ic]+1,ecut[ic],polygon) else: add_secpolygon_list(cutpolygon,icut[ic]+1,ecut[ic],polygon) iipc = iipc+dcpt-1 # cutting line cutline = np.zeros((Nadds[ic],2), dtype=np.float) dx = (ept[ic][1] - ipt[ic][1])/(Nadds[ic]-1) dy = (ept[ic][0] - ipt[ic][0])/(Nadds[ic]-1) cutline[0,:] = ipt[ic] for ip in range(1,Nadds[ic]-1): cutline[ip,:] = ipt[ic] + np.array([dy*ip,dx*ip]) cutline[Nadds[ic]-1,:] = ept[ic] if keep == 'left': for ip in range(Nadds[ic]-1,-1,-1): cutpolygon.append(cutline[ip,:]) iipc = iipc+Nadds[ic] if ic == 0: add_secpolygon_list(cutpolygon,ecut[ic],N,polygon) cutpolygon.append(polygon[0,:]) iipc = iipc + 1 else: for ip in range(Nadds[ic]-1,-1,-1): cutpolygon.append(cutline[ip,:]) iipc = iipc+Nadds[ic] cutpolygon.append([gen.fillValueF, gen.fillValueF]) iipc = iipc + 1 cutpolygon = np.array(cutpolygon) rmpolygon = [] Npts = cutpolygon.shape[0] if keep == 'left': for ip in range(Npts): if cutpolygon[ip,1] > xval: rmpolygon.append([gen.fillValueF, gen.fillValueF]) else: rmpolygon.append(cutpolygon[ip,:]) else: for ip in range(Npts): if cutpolygon[ip,1] < xval: rmpolygon.append([gen.fillValueF, gen.fillValueF]) else: rmpolygon.append(cutpolygon[ip,:]) rmpolygon = np.array(rmpolygon) cutpolygon = rm_consecpt_polygon(rmpolygon) Npts = cutpolygon.shape[0] cutpolygon = ma.masked_equal(cutpolygon, gen.fillValueF) return Npts, cutpolygon def cut_between_ypolygon(polygon, yval1, yval2, Nadd=20): """ Function to cut a polygon between 2 given value of the y-axis polygon: polygon to cut yval1: first value to use to cut the polygon yval2: first value to use to cut the polygon Nadd: additional points to add to draw the line (20, default) """ fname = 'cut_betwen_ypolygon' N = polygon.shape[0] if yval1 > yval2: print errormsg print ' ' + fname + ': wrong between cut values !!' print ' it is expected yval1 < yval2' print ' values provided yval1: (', yval1, ')> yval2 (', yval2, ')' quit(-1) yvals = [yval1, yval2] ipt = None ept = None cuts = {} if type(polygon) == type(gen.mamat) and type(polygon.mask) != \ type(gen.mamat.mask[1]): for ic in range(2): yval = yvals[ic] # There might be more than 1 cut ... icut = [] ecut = [] ipt = [] ept = [] Ncuts = 0 # Assuming clockwise polygons for ip in range(N-1): if not polygon.mask[ip,0]: eep = ip + 1 if eep == N: eep = 0 if val_consec_between(polygon[ip,0], polygon[eep,0], yval): icut.append(ip) dx = polygon[eep,1] - polygon[ip,1] dy = polygon[eep,0] - polygon[ip,0] dd = yval - polygon[ip,0] ipt.append([yval, polygon[ip,1]+dx*dd/dy]) if val_consec_between(polygon[eep,0], polygon[ip,0], yval): ecut.append(ip) dx = polygon[eep,1] - polygon[ip,1] dy = polygon[eep,0] - polygon[ip,0] dd = yval - polygon[ip,0] ept.append([yval, polygon[ip,1]+dx*dd/dy]) Ncuts = Ncuts + 1 # Looking for repeated newicut = icut + [] newecut = ecut + [] newipt = ipt + [] newept = ept + [] newNcuts = Ncuts for icp in range(newNcuts-1): for ic2 in range(icp+1,newNcuts): if newipt[icp] == newipt[ic2]: Ncuts = Ncuts-1 icut.pop(ic2) ecut.pop(ic2) ipt.pop(ic2) ept.pop(ic2) newNcuts = Ncuts + 0 cuts[ic] = [icut, ecut, ipt, ept, Ncuts] else: for ic in range(2): yval = yvals[ic] # There might be more than 1 cut ... icut = [] ecut = [] ipt = [] ept = [] Ncuts = 0 # Assuming clockwise polygons for ip in range(N-1): eep = ip + 1 if eep == N: eep = 0 if val_consec_between(polygon[ip,0], polygon[eep,0], yval): icut.append(ip) dx = polygon[eep,1] - polygon[ip,1] dy = polygon[eep,0] - polygon[ip,0] dd = yval - polygon[ip,0] ipt.append([yval, polygon[ip,1]+dx*dd/dy]) if val_consec_between(polygon[eep,0], polygon[ip,0], yval): ecut.append(ip) dx = polygon[eep,1] - polygon[ip,1] dy = polygon[eep,0] - polygon[ip,0] dd = yval - polygon[ip,0] ept.append([yval, polygon[ip,1]+dx*dd/dy]) Ncuts = Ncuts + 1 # Looking for repeated newicut = icut + [] newecut = ecut + [] newipt = ipt + [] newept = ept + [] newNcuts = Ncuts for icp in range(newNcuts-1): for ic2 in range(icp+1,newNcuts): if newipt[icp] == newipt[ic2]: Ncuts = Ncuts-1 icut.pop(ic2) ecut.pop(ic2) ipt.pop(ic2) ept.pop(ic2) newNcuts = Ncuts + 0 cuts[ic] = [icut, ecut, ipt, ept, Ncuts] Naddlines = {} for icc in range(2): cutv = cuts[icc] Ncuts = cutv[4] # Length of joining lines Nadds = [] if Ncuts > 1: Naddc = (Nadd-Ncuts)/(Ncuts) if Naddc < 3: print errormsg print ' ' + fname + ': too few points for jioning lines !!' print ' increase Nadd at least to:', Ncuts*3+Ncuts quit(-1) for ic in range(Ncuts-1): Nadds.append(Naddc) Nadds.append(Nadd-Naddc*(Ncuts-1)) else: Nadds.append(Nadd) # Total points cut polygon Ntotpts = 0 Ncpts = [] for ic in range(Ncuts): dpts = ecut[ic] - icut[ic] + Nadds[ic] - 1 Ntotpts = Ntotpts + dpts Ncpts.append(ecut[ic] - icut[ic]) Naddlines[icc] = [Nadds, Ntotpts, Ncpts] cutv1 = cuts[0] addv1 = Naddlines[0] Nadds1 = addv1[0] Ncuts1 = cutv1[4] cutv2 = cuts[1] addv2 = Naddlines[1] Nadds2 = addv2[0] Ncuts2 = cutv2[4] if Ncuts1 != Ncuts2: print errormsg print ' ' + fname + ": different number of cuts !!" print ' yval1:', yval1, 'Ncuts=', Ncuts1 print ' yval2:', yval2, 'Ncuts=', Ncuts2 print ' I am not prepare to deal with it' quit(-1) #else: # print ' ' + fname + ' _______' # print ' yval1:', yval1, 'Ncuts=', Ncuts1 # print ' yval2:', yval2, 'Ncuts=', Ncuts2 icut1 = cutv1[0] ecut1 = cutv1[1] ipt1 = cutv1[2] ept1 = cutv1[3] icut2 = cutv2[0] ecut2 = cutv2[1] ipt2 = cutv2[2] ept2 = cutv2[3] # Looking for pairs of cuts. Grouping for smallest x distance between initial # points of each cut cutpolygons = [] for iic1 in range(Ncuts1): iip = 0 cutpolygon = [] ic1 = icut1[iic1] ec1 = ecut1[iic1] ip1 = ipt1[iic1] ep1 = ept1[iic1] ipx2s = [] for ip in range(Ncuts2): ip2 = ipt2[ip] ipx2s.append(ip2[1]) dxps = ipx2s - ip1[1] dxps = np.where(dxps < 0., gen.fillValueF, dxps) ndxps = np.min(dxps) iic12 = gen.index_vec(dxps,ndxps) ic2 = icut2[iic12] ec2 = ecut2[iic12] ip2 = ipt2[iic12] ep2 = ept2[iic12] #print 'Lluis iic1', iic1, 'ic1', ic1, 'ec1', ec1, 'ipt1', ip1, 'ept1', ep1, 'Nadds1', Nadds1 #print ' iic12', iic12, 'ic2', ic2, 'ec2', ec2, 'ipt2', ip2, 'ept2', ep2, 'Nadds2', Nadds2 cutpolygon.append(ip1) for ip in range(ic1+1,ic2-1): cutpolygon.append(polygon[ip,:]) iip = ic2-ic1 # cutting line 1 Nadd2 = Nadds1[iic1] cutlines = np.zeros((Nadd2,2), dtype=np.float) dx = (ep2[1] - ip2[1])/(Nadd2-2) dy = (ep2[0] - ip2[0])/(Nadd2-2) cutlines[0,:] = ip2 for ip in range(1,Nadd2-1): cutlines[ip,:] = ip2 + np.array([dy*ip,dx*ip]) cutlines[Nadd2-1,:] = ep2 for ip in range(Nadd2): cutpolygon.append(cutlines[ip,:]) iip = iip + Nadd2 for ip in range(ec2,ec1): cutpolygon.append(polygon[ip,:]) iip = iip + ec1-ec2 # cutting line 2 Nadd2 = Nadds2[iic12] cutlines = np.zeros((Nadd2,2), dtype=np.float) dx = (ep1[1] - ip1[1])/(Nadd2-2) dy = (ep1[0] - ip1[0])/(Nadd2-2) cutlines[0,:] = ip1 for ip in range(1,Nadd2-1): cutlines[ip,:] = ip1 + np.array([dy*ip,dx*ip]) cutlines[Nadd2-1,:] = ep1 for ip in range(Nadd2-1,0,-1): cutpolygon.append(cutlines[ip,:]) cutpolygon.append(ip1) cutpolygon.append([gen.fillValueF,gen.fillValueF]) if len(cutpolygons) == 0: cutpolygons = cutpolygon else: cutpolygons = cutpolygons + cutpolygon cutpolygons = np.array(cutpolygons) cutpolygons = rm_consecpt_polygon(cutpolygons) cutpolygons = ma.masked_equal(cutpolygons, gen.fillValueF) Npts = cutpolygons.shape[0] return Npts, cutpolygons def cut_between_xpolygon(polygon, xval1, xval2, Nadd=20): """ Function to cut a polygon between 2 given value of the x-axis polygon: polygon to cut xval1: first value to use to cut the polygon xval2: first value to use to cut the polygon Nadd: additional points to add to draw the line (20, default) """ fname = 'cut_betwen_xpolygon' N = polygon.shape[0] if xval1 > xval2: print errormsg print ' ' + fname + ': wrong between cut values !!' print ' it is expected xval1 < xval2' print ' values provided xval1: (', xval1, ')> xval2 (', xval2, ')' quit(-1) xvals = [xval1, xval2] ipt = None ept = None cuts = {} if type(polygon) == type(gen.mamat) and type(polygon.mask) != \ type(gen.mamat.mask[1]): for ic in range(2): xval = xvals[ic] # There might be more than 1 cut ... icut = [] ecut = [] ipt = [] ept = [] Ncuts = 0 # Assuming clockwise polygons for ip in range(N-1): if not polygon.mask[ip,0]: eep = ip + 1 if eep == N: eep = 0 if val_consec_between(polygon[ip,1], polygon[eep,1], xval): icut.append(ip) dx = polygon[eep,1] - polygon[ip,1] dy = polygon[eep,0] - polygon[ip,0] dd = xval - polygon[ip,1] ipt.append([polygon[ip,0]+dy*dd/dx, xval]) if val_consec_between(polygon[eep,1], polygon[ip,1], xval): ecut.append(ip) dx = polygon[eep,1] - polygon[ip,1] dy = polygon[eep,0] - polygon[ip,0] dd = xval - polygon[ip,1] ept.append([polygon[ip,0]+dy*dd/dx, xval]) Ncuts = Ncuts + 1 # Looking for repeated newicut = icut + [] newecut = ecut + [] newipt = ipt + [] newept = ept + [] newNcuts = Ncuts for icp in range(newNcuts-1): for ic2 in range(icp+1,newNcuts): if newipt[icp] == newipt[ic2]: Ncuts = Ncuts-1 icut.pop(ic2) ecut.pop(ic2) ipt.pop(ic2) ept.pop(ic2) newNcuts = Ncuts + 0 cuts[ic] = [icut, ecut, ipt, ept, Ncuts] else: for ic in range(2): xval = xvals[ic] # There might be more than 1 cut ... icut = [] ecut = [] ipt = [] ept = [] Ncuts = 0 # Assuming clockwise polygons for ip in range(N-1): eep = ip + 1 if eep == N: eep = 0 if val_consec_between(polygon[ip,1], polygon[eep,1], xval): icut.append(ip) dx = polygon[eep,1] - polygon[ip,1] dy = polygon[eep,0] - polygon[ip,0] dd = xval - polygon[ip,1] ipt.append([polygon[ip,0]+dy*dd/dx, xval]) if val_consec_between(polygon[eep,1], polygon[ip,1], xval): ecut.append(ip) dx = polygon[eep,1] - polygon[ip,1] dy = polygon[eep,0] - polygon[ip,0] dd = xval - polygon[ip,1] ept.append([polygon[ip,0]+dy*dd/dx, xval]) Ncuts = Ncuts + 1 # Looking for repeated newicut = icut + [] newecut = ecut + [] newipt = ipt + [] newept = ept + [] newNcuts = Ncuts for icp in range(newNcuts-1): for ic2 in range(icp+1,newNcuts): if newipt[icp] == newipt[ic2]: Ncuts = Ncuts-1 icut.pop(ic2) ecut.pop(ic2) ipt.pop(ic2) ept.pop(ic2) newNcuts = Ncuts + 0 cuts[ic] = [icut, ecut, ipt, ept, Ncuts] for iic in range(1): cutvs = cuts[iic] icut = cutvs[0] ecut = cutvs[1] ipt = cutvs[2] ept = cutvs[3] Ncuts = cutvs[4] if Ncuts > 0: # Re-shifting cuts by closest heigth. yis = [] yes = [] for ic in range(Ncuts): yp = ipt[ic] yis.append(yp[0]) yp = ept[ic] yes.append(yp[0]) ys = yis + yes ys.sort() newicut = icut + [] newecut = ecut + [] newipt = ipt + [] newept = ept + [] icut = [] ecut = [] ipt = [] ept = [] for yv in ys: ic = yis.count(yv) if ic != 0: icc = yis.index(yv) if len(icut) > len(ecut): ecut.append(newicut[icc]) ept.append(newipt[icc]) else: icut.append(newicut[icc]) ipt.append(newipt[icc]) else: icc = yes.index(yv) if len(icut) > len(ecut): ecut.append(newecut[icc]) ept.append(newept[icc]) else: icut.append(newecut[icc]) ipt.append(newept[icc]) cuts[iic] = [icut, ecut, ipt, ept, Ncuts] Naddlines = {} for icc in range(2): cutv = cuts[icc] Ncuts = cutv[4] if Ncuts == 0: print errormsg print ' ' + fname + ": no cuts for xval=", xvals[icc], '!!' quit(-1) #print ' icc:', icc, 'ic ec ipt ept _______' #for ic in range(Ncuts): # print ic, ':', cutv[0][ic], cutv[1][ic], cutv[2][ic], cutv[3][ic] # Length of joining lines Nadds = [] if Ncuts > 1: Naddc = (Nadd-Ncuts)/(Ncuts) if Naddc < 3: print errormsg print ' ' + fname + ': too few points for jioning lines !!' print ' increase Nadd at least to:', Ncuts*3+Ncuts quit(-1) for ic in range(Ncuts-1): Nadds.append(Naddc) Nadds.append(Nadd-Naddc*(Ncuts-1)) else: Nadds.append(Nadd) Naddlines[icc] = Nadds # sides sides = {} for iic in range(2): cutvs = cuts[iic] icut = cutvs[0] ecut = cutvs[1] ipt = cutvs[2] ept = cutvs[3] Ncuts = cutvs[4] Nadds = Naddlines[iic] cutpolygon = [] # left side if iic == 0: for ic in range(Ncuts-1): cutpolygon.append(ipt[ic]) dx = (ept[ic][1] - ipt[ic][1])/(Nadds[ic]-1) dy = (ept[ic][0] - ipt[ic][0])/(Nadds[ic]-1) for ip in range(1,Nadds[ic]-1): cutpolygon.append([ipt[ic][0]+dy*ip, ipt[ic][1]+dx*ip]) cutpolygon.append(ept[ic]) for ip in range(ecut[ic]+1,icut[ic+1]): cutpolygon.append(polygon[ip,:]) ic = Ncuts-1 cutpolygon.append(ipt[ic]) dx = (ept[ic][1] - ipt[ic][1])/(Nadds[ic]-1) dy = (ept[ic][0] - ipt[ic][0])/(Nadds[ic]-1) for ip in range(1,Nadds[ic]-1): cutpolygon.append([ipt[ic][0]+dy*ip, ipt[ic][1]+dx*ip]) # right side else: for ic in range(Ncuts-1): cutpolygon.append(ipt[ic]) # line dx = (ept[ic][1] - ipt[ic][1])/(Nadds[ic]-1) dy = (ept[ic][0] - ipt[ic][0])/(Nadds[ic]-1) for ip in range(1,Nadds[ic]-1): cutpolygon.append([ipt[ic][0]+dy*ip, ipt[ic][1]+dx*ip]) cutpolygon.append(ept[ic]) for ip in range(ecut[ic],icut[ic+1]): cutpolygon.append(polygon[ip,:]) ic = Ncuts-1 cutpolygon.append(ipt[ic]) # line dx = (ept[ic][1] - ipt[ic][1])/(Nadds[ic]-1) dy = (ept[ic][0] - ipt[ic][0])/(Nadds[ic]-1) for ip in range(1,Nadds[ic]-1): cutpolygon.append([ipt[ic][0]+dy*ip, ipt[ic][1]+dx*ip]) cutpolygon.append(ept[ic]) sides[iic] = cutpolygon # joining sides by e1[Ncuts1-1] --> i2[0]; e2[Ncuts2-1] --> i1[0] cutv1 = cuts[0] Ncuts1 = cutv1[4] ec1 = cutv1[1][np.max([0,Ncuts1-1])] ic1 = cutv1[0][0] ept1 = cutv1[3][np.max([0,Ncuts1-1])] ipt1 = cutv1[2][0] cutv2 = cuts[1] Ncuts2 = cutv2[4] ec2 = cutv2[1][np.max([0,Ncuts2-1])] ic2 = cutv2[0][0] ept2 = cutv2[3][np.max([0,Ncuts2-1])] ipt2 = cutv2[2][0] finalcutpolygon = sides[0] for ip in range(ec1+1,ic2): finalcutpolygon.append(polygon[ip,:]) finalcutpolygon = finalcutpolygon + sides[1] for ip in range(ec2+1,ic1): finalcutpolygon.append(polygon[ip,:]) finalcutpolygon.append(ipt1) finalcutpolygon = np.array(finalcutpolygon) finalcutpolygon = rm_consecpt_polygon(finalcutpolygon) finalcutpolygon = ma.masked_equal(finalcutpolygon, gen.fillValueF) Npts = finalcutpolygon.shape[0] return Npts, finalcutpolygon def pile_polygons(polyns, polygons): """ Function to pile polygons one over the following one polyns: ordered list of polygons. First over all. last below all polygons: dictionary with the polygons >>> pns = ['sqra', 'sqrb'] >>> polya = np.array([[-0.5, -0.75], [0.5, -0.75], [0.5, 0.75], [-0.5, 0.75]]) >>> polyb = np.array([[-0.75, -0.5], [0.75, -0.5], [0.75, 0.5], [-0.75, 0.5]]) >>> plgs = {'sqra': polya, 'sqrb': polyb} >>> pile_polygons(pns, plgs) # sqrb : [[-0.75 -0.5] [-0.5 -0.5] [-- --] [0.5 -0.5] [0.75 -0.5] [0.75 0.5] [0.5 0.5] [-- --] [-0.5 0.5] [-0.75 0.5] [-0.75 -0.5]] # sqra : [[-0.5 -0.75] [ 0.5 -0.75] [ 0.5 0.75] [-0.5 0.75] [-0.5 -0.75]] """ fname = 'pile_polygons' pilepolygons = dict(polygons) Npolys = len(polyns) for ipolyp in range(Npolys-2,-1,-1): polyn = polyns[ipolyp] poly = pilepolygons[polyn] Npts = poly.shape[0] for ipolyi in range(ipolyp+1,Npolys,1): ipolyn = polyns[ipolyi] #print ' Lluis ' + polyn + ' above ' + ipolyn ipoly = pilepolygons[ipolyn] iNpts = ipoly.shape[0] newipoly = [] Nint, inti, intp, pts = fsci.module_scientific.crossingpoints_polys( \ nvertexa=iNpts, nvertexb=Npts, nvertexab=iNpts*Npts, polya=ipoly, \ polyb=poly) # We're in C-mode ! inti = inti-1 intp = intp-1 # Re-constructing below polygon looking for respective crossings linti = list(inti) for ip in range(iNpts): iip1 = ip+1 if ip == iNpts-1: iip1 = 0 #print ip, ipoly[ip,:], ':', ipoly[iip1,:] Nc = linti.count(ip) ldists = [] ddists = {} if Nc > 0: iic = gen.multi_index_vec(inti,ip) mindist = 1000000. # Sorting from distance respect the vertex ip for ic in range(Nc): ddists[iic[ic]] = dist_points(ipoly[ip,:], pts[iic[ic],:]) ldists.append(ddists[iic[ic]]) #print ' ', ic, ';', iic[ic], '=', pts[iic[ic],:], ddists[iic[ic]] ldists.sort() #print ' ldists', ldists newipoly.append(ipoly[ip,:]) for ic in range(Nc): iic = gen.dictionary_key(ddists,ldists[ic]) newipoly.append(pts[iic,:]) #print ' ', ic, '|', iic, ';', pts[iic,:] if ic < Nc-1: newipoly.append([gen.fillValueF, gen.fillValueF]) newipoly.append(ipoly[iip1,:]) else: newipoly.append(ipoly[ip,:]) newipoly = np.array(newipoly) pilepolygons[polyns[Npolys-1]] = rm_consecpt_polygon(newipoly) for polyn in polyns: poly = pilepolygons[polyn] poly = ma.masked_equal(poly, gen.fillValueF) pilepolygons[polyn] = poly return pilepolygons ####### ###### ##### #### ### ## # # Shapes/objects def surface_sphere(radii,Npts): """ Function to provide an sphere as matrix of x,y,z coordinates radii: radii of the sphere Npts: number of points to discretisize longitues (half for latitudes) """ fname = 'surface_sphere' sphereup = np.zeros((3,Npts/2,Npts), dtype=np.float) spheredown = np.zeros((3,Npts/2,Npts), dtype=np.float) for ia in range(Npts): alpha = ia*2*np.pi/(Npts-1) for ib in range(Npts/2): beta = ib*np.pi/(2.*(Npts/2-1)) sphereup[:,ib,ia] = position_sphere(radii, alpha, beta) for ib in range(Npts/2): beta = -ib*np.pi/(2.*(Npts/2-1)) spheredown[:,ib,ia] = position_sphere(radii, alpha, beta) return sphereup, spheredown def ellipse_polar(c, a, b, Nang=100): """ Function to determine an ellipse from its center and polar coordinates FROM: https://en.wikipedia.org/wiki/Ellipse c= coordinates of the center a= distance major axis b= distance minor axis Nang= number of angles to use """ fname = 'ellipse_polar' if np.mod(Nang,2) == 0: Nang=Nang+1 dtheta = 2*np.pi/(Nang-1) ellipse = np.zeros((Nang,2), dtype=np.float) for ia in range(Nang): theta = dtheta*ia rad = a*b/np.sqrt( (b*np.cos(theta))**2 + (a*np.sin(theta))**2 ) x = rad*np.cos(theta) y = rad*np.sin(theta) ellipse[ia,:] = [y+c[0],x+c[1]] return ellipse def hyperbola_polar(a, b, Nang=100): """ Fcuntion to determine an hyperbola in polar coordinates FROM: https://en.wikipedia.org/wiki/Hyperbola#Polar_coordinates x^2/a^2 - y^2/b^2 = 1 a= x-parameter y= y-parameter Nang= number of angles to use DOES NOT WORK!!!! """ fname = 'hyperbola_polar' dtheta = 2.*np.pi/(Nang-1) # Positive branch hyperbola_p = np.zeros((Nang,2), dtype=np.float) for ia in range(Nang): theta = dtheta*ia x = a*np.cosh(theta) y = b*np.sinh(theta) hyperbola_p[ia,:] = [y,x] # Negative branch hyperbola_n = np.zeros((Nang,2), dtype=np.float) for ia in range(Nang): theta = dtheta*ia x = -a*np.cosh(theta) y = b*np.sinh(theta) hyperbola_n[ia,:] = [y,x] return hyperbola_p, hyperbola_n def circ_sec(ptA, ptB, radii, arc='short', pos='left', Nang=100): """ Function union of point A and B by a section of a circle ptA= coordinates od the point A [yA, xA] ptB= coordinates od the point B [yB, xB] radii= radi of the circle to use to unite the points arc= which arc to be used ('short', default) 'short': shortest angle between points 'long': largest angle between points pos= orientation of the arc following clockwise union of points ('left', default) 'left': to the left of union 'right': to the right of union Nang= amount of angles to use """ fname = 'circ_sec' availarc = ['short', 'long'] availpos = ['left', 'right'] distAB = dist_points(ptA,ptB) if distAB > radii: print errormsg print ' ' + fname + ': radii=', radii, " too small for the distance " + \ "between points !!" print ' distance between points:', distAB quit(-1) # Coordinate increments dAB = np.abs(ptA-ptB) # angle of the circular section joining points alpha = 2.*np.arcsin((distAB/2.)/radii) # center along coincident bisection of the union xcc = -radii ycc = 0. # Getting the arc of the circle at the x-axis if arc == 'short': dalpha = alpha/(Nang-1) elif arc == 'long': dalpha = (2.*np.pi - alpha)/(Nang-1) else: print errormsg print ' ' + fname + ": arc '" + arc + "' not ready !!" print ' available ones:', availarc quit(-1) if pos == 'left': sign=-1. elif pos == 'right': sign=1. else: print errormsg print ' ' + fname + ": position '" + pos + "' not ready !!" print ' available ones:', availpos quit(-1) circ_sec = np.zeros((Nang,2), dtype=np.float) for ia in range(Nang): alpha = sign*dalpha*ia x = radii*np.cos(alpha) y = radii*np.sin(alpha) circ_sec[ia,:] = [y+ycc,x+xcc] # Angle of the points theta = np.arctan2(ptB[0]-ptA[0],ptB[1]-ptA[1]) # rotating angle of the circ if pos == 'left': rotangle = theta + np.pi/2. - alpha/2. elif pos == 'right': rotangle = theta + 3.*np.pi/2. - alpha/2. else: print errormsg print ' ' + fname + ": position '" + pos + "' not ready !!" print ' available ones:', availpos quit(-1) #print 'alpha:', alpha*180./np.pi, 'theta:', theta*180./np.pi, 'rotangle:', rotangle*180./np.pi # rotating the arc along the x-axis rotcirc_sec = rotate_polygon_2D(circ_sec, rotangle) # Moving arc to the ptA circ_sec = rotcirc_sec + ptA return circ_sec def p_square(face, N=5): """ Function to get a polygon square face: length of the face of the square N: number of points of the polygon """ fname = 'p_square' square = np.zeros((N,2), dtype=np.float) f2 = face/2. N4 = N/4 df = face/(N4) # SW-NW for ip in range(N4): square[ip,:] = [-f2+ip*df,-f2] # NW-NE for ip in range(N4): square[ip+N4,:] = [f2,-f2+ip*df] # NE-SE for ip in range(N4): square[ip+2*N4,:] = [f2-ip*df,f2] N42 = N-3*N4-1 df = face/(N42) # SE-SW for ip in range(N42): square[ip+3*N4,:] = [-f2,f2-ip*df] square[N-1,:] = [-f2,-f2] return square def p_prism(base, height, N=5): """ Function to get a polygon prism base: length of the base of the prism height: length of the height of the prism N: number of points of the polygon """ fname = 'p_prism' prism = np.zeros((N,2), dtype=np.float) b2 = base/2. h2 = height/2. N4 = N/4 dh = height/(N4) db = base/(N4) # SW-NW for ip in range(N4): prism[ip,:] = [-h2+ip*dh,-b2] # NW-NE for ip in range(N4): prism[ip+N4,:] = [h2,-b2+ip*db] # NE-SE for ip in range(N4): prism[ip+2*N4,:] = [h2-ip*dh,b2] N42 = N-3*N4-1 db = base/(N42) # SE-SW for ip in range(N42): prism[ip+3*N4,:] = [-h2,b2-ip*db] prism[N-1,:] = [-h2,-b2] return prism def p_circle(radii, N=50): """ Function to get a polygon of a circle radii: length of the radii of the circle N: number of points of the polygon """ fname = 'p_circle' circle = np.zeros((N,2), dtype=np.float) dangle = 2.*np.pi/(N-1) for ia in range(N): circle[ia,:] = [radii*np.sin(ia*dangle), radii*np.cos(ia*dangle)] circle[N-1,:] = [0., radii] return circle def p_triangle(p1, p2, p3, N=4): """ Function to provide the polygon of a triangle from its 3 vertices p1: vertex 1 [y,x] p2: vertex 2 [y,x] p3: vertex 3 [y,x] N: number of vertices of the triangle """ fname = 'p_triangle' triangle = np.zeros((N,2), dtype=np.float) N3 = N / 3 # 1-2 dx = (p2[1]-p1[1])/N3 dy = (p2[0]-p1[0])/N3 for ip in range(N3): triangle[ip,:] = [p1[0]+ip*dy,p1[1]+ip*dx] # 2-3 dx = (p3[1]-p2[1])/N3 dy = (p3[0]-p2[0])/N3 for ip in range(N3): triangle[ip+N3,:] = [p2[0]+ip*dy,p2[1]+ip*dx] # 3-1 N32 = N - 2*N/3 dx = (p1[1]-p3[1])/N32 dy = (p1[0]-p3[0])/N32 for ip in range(N32): triangle[ip+2*N3,:] = [p3[0]+ip*dy,p3[1]+ip*dx] triangle[N-1,:] = p1 return triangle def p_spiral(loops, eradii, N=1000): """ Function to provide a polygon of an Archimedean spiral FROM: https://en.wikipedia.org/wiki/Spiral loops: number of loops of the spiral eradii: length of the radii of the final spiral N: number of points of the polygon """ fname = 'p_spiral' spiral = np.zeros((N,2), dtype=np.float) dangle = 2.*np.pi*loops/(N-1) dr = eradii*1./(N-1) for ia in range(N): radii = dr*ia spiral[ia,:] = [radii*np.sin(ia*dangle), radii*np.cos(ia*dangle)] return spiral def p_reg_polygon(Nv, lf, N=50): """ Function to provide a regular polygon of Nv vertices Nv: number of vertices lf: length of the face N: number of points """ fname = 'p_reg_polygon' reg_polygon = np.zeros((N,2), dtype=np.float) # Number of points per vertex Np = N/Nv # Angle incremental between vertices da = 2.*np.pi/Nv # Radii of the circle according to lf radii = lf*Nv/(2*np.pi) iip = 0 for iv in range(Nv-1): # Characteristics between vertices iv and iv+1 av1 = da*iv v1 = [radii*np.sin(av1), radii*np.cos(av1)] av2 = da*(iv+1) v2 = [radii*np.sin(av2), radii*np.cos(av2)] dx = (v2[1]-v1[1])/Np dy = (v2[0]-v1[0])/Np for ip in range(Np): reg_polygon[ip+iv*Np,:] = [v1[0]+dy*ip,v1[1]+dx*ip] # Characteristics between vertices Nv and 1 # Number of points per vertex Np2 = N - Np*(Nv-1) av1 = da*Nv v1 = [radii*np.sin(av1), radii*np.cos(av1)] av2 = 0. v2 = [radii*np.sin(av2), radii*np.cos(av2)] dx = (v2[1]-v1[1])/Np2 dy = (v2[0]-v1[0])/Np2 for ip in range(Np2): reg_polygon[ip+(Nv-1)*Np,:] = [v1[0]+dy*ip,v1[1]+dx*ip] return reg_polygon def p_reg_star(Nv, lf, freq, vs=0, N=50): """ Function to provide a regular star of Nv vertices Nv: number of vertices lf: length of the face of the regular polygon freq: frequency of union of vertices ('0', for just centered to zero arms) vs: vertex from which start (0 being first [0,lf]) N: number of points """ fname = 'p_reg_star' reg_star = np.zeros((N,2), dtype=np.float) # Number of arms of the star if freq != 0 and np.mod(Nv,freq) == 0: Na = Nv/freq + 1 else: Na = Nv # Number of points per arm Np = N/Na # Angle incremental between vertices da = 2.*np.pi/Nv # Radii of the circle according to lf radii = lf*Nv/(2*np.pi) iip = 0 av1 = vs*da for iv in range(Na-1): # Characteristics between vertices iv and iv+1 v1 = [radii*np.sin(av1), radii*np.cos(av1)] if freq != 0: av2 = av1 + da*freq v2 = [radii*np.sin(av2), radii*np.cos(av2)] else: v2 = [0., 0.] av2 = av1 + da dx = (v2[1]-v1[1])/(Np-1) dy = (v2[0]-v1[0])/(Np-1) for ip in range(Np): reg_star[ip+iv*Np,:] = [v1[0]+dy*ip,v1[1]+dx*ip] if av2 > 2.*np.pi: av1 = av2 - 2.*np.pi else: av1 = av2 + 0. iv = Na-1 # Characteristics between vertices Na and 1 Np2 = N-Np*iv v1 = [radii*np.sin(av1), radii*np.cos(av1)] if freq != 0: av2 = vs*da v2 = [radii*np.sin(av2), radii*np.cos(av2)] else: v2 = [0., 0.] dx = (v2[1]-v1[1])/(Np2-1) dy = (v2[0]-v1[0])/(Np2-1) for ip in range(Np2): reg_star[ip+iv*Np,:] = [v1[0]+dy*ip,v1[1]+dx*ip] return reg_star def p_sinusoide(length=10., amp=5., lamb=3., ival=0., func='sin', N=100): """ Function to get coordinates of a sinusoidal curve length: length of the line (default 10.) amp: amplitude of the peaks (default 5.) lamb: wave longitude (defalult 3.) ival: initial angle (default 0. in degree) func: function to use: (default sinus) 'sin': sinus 'cos': cosinus N: number of points (default 100) """ fname = 'p_sinusoide' availfunc = ['sin', 'cos'] dx = length/(N-1) ia = ival*np.pi/180. da = 2*np.pi*dx/lamb sinusoide = np.zeros((N,2), dtype=np.float) if func == 'sin': for ix in range(N): sinusoide[ix,:] = [amp*np.sin(ia+da*ix),dx*ix] elif func == 'cos': for ix in range(N): sinusoide[ix,:] = [amp*np.cos(ia+da*ix),dx*ix] else: print errormsg print ' ' + fname + ": function '" + func + "' not ready !!" print ' available ones:', availfunc quit(-1) sinusoidesecs = ['sinusoide'] sinusoidedic = {'sinusoide': [sinusoide, '-', '#000000', 1.]} return sinusoide, sinusoidesecs, sinusoidedic def p_doubleArrow(length=5., angle=45., width=1., alength=0.10, N=50): """ Function to provide an arrow with double lines length: length of the arrow (5. default) angle: angle of the head of the arrow (45., default) width: separation between the two lines (2., default) alength: length of the head (as percentage in excess of width, 0.1 default) N: number of points (50, default) """ function = 'p_doubleArrow' doubleArrow = np.zeros((50,2), dtype=np.float) N4 = int((N-3)/4) doublearrowdic = {} ddy = width*np.tan(angle*np.pi/180.)/2. # Arms dx = (length-ddy)/(N4-1) for ix in range(N4): doubleArrow[ix,:] = [dx*ix,-width/2.] doublearrowdic['leftarm'] = [doubleArrow[0:N4,:], '-', '#000000', 2.] doubleArrow[N4,:] = [gen.fillValueF,gen.fillValueF] for ix in range(N4): doubleArrow[N4+1+ix,:] = [dx*ix,width/2.] doublearrowdic['rightarm'] = [doubleArrow[N4+1:2*N4+1,:], '-', '#000000', 2.] doubleArrow[2*N4+1,:] = [gen.fillValueF,gen.fillValueF] # Head N42 = int((N-2 - 2*N4)/2) dx = width*(1.+alength)*np.cos(angle*np.pi/180.)/(N42-1) dy = width*(1.+alength)*np.sin(angle*np.pi/180.)/(N42-1) for ix in range(N42): doubleArrow[2*N4+2+ix,:] = [length-dy*ix,-dx*ix] doublearrowdic['lefthead'] = [doubleArrow[2*N4:2*N4+N42,:], '-', '#000000', 2.] doubleArrow[2*N4+2+N42,:] = [gen.fillValueF,gen.fillValueF] N43 = N-3 - 2*N4 - N42 + 1 dx = width*(1.+alength)*np.cos(angle*np.pi/180.)/(N43-1) dy = width*(1.+alength)*np.sin(angle*np.pi/180.)/(N43-1) for ix in range(N43): doubleArrow[2*N4+N42+2+ix,:] = [length-dy*ix,dx*ix] doublearrowdic['rightthead'] = [doubleArrow[2*N4+N42+2:51,:], '-', '#000000', 2.] doubleArrow = ma.masked_equal(doubleArrow, gen.fillValueF) doublearrowsecs = ['leftarm', 'rightarm', 'lefthead', 'righthead'] return doubleArrow, doublearrowsecs, doublearrowdic def p_angle_triangle(pi=np.array([0.,0.]), angle1=60., length1=1., angle2=60., N=100): """ Function to draw a triangle by an initial point and two consecutive angles and the first length of face. The third angle and 2 and 3rd face will be computed accordingly the provided values: length1 / sin(angle1) = length2 / sin(angle2) = length3 / sin(angle3) angle1 + angle2 + angle3 = 180. pi: initial point ([0., 0.], default) angle1: first angle from pi clockwise (60., default) length1: length of face from pi by angle1 (1., default) angle2: second angle from second point (60., default) length2: length of face from p2 by angle2 (1., default) N: number of points (100, default) """ fname = 'p_angle_triangle' angle3 = 180. - angle1 - angle2 length2 = np.sin(angle2*np.pi/180.)*length1/np.sin(angle1*np.pi/180.) length3 = np.sin(angle3*np.pi/180.)*length1/np.sin(angle1*np.pi/180.) triangle = np.zeros((N,2), dtype=np.float) N3 = int(N/3) # first face ix = pi[1] iy = pi[0] dx = length1*np.cos(angle1*np.pi/180.)/(N3-1) dy = length1*np.sin(angle1*np.pi/180.)/(N3-1) for ip in range(N3): triangle[ip,:] = [iy+dy*ip, ix+dx*ip] # second face ia = -90. - (90.-angle1) ix = triangle[N3-1,1] iy = triangle[N3-1,0] dx = length2*np.cos((ia+angle2)*np.pi/180.)/(N3-1) dy = length2*np.sin((ia+angle2)*np.pi/180.)/(N3-1) for ip in range(N3): triangle[N3+ip,:] = [iy+dy*ip, ix+dx*ip] # third face N32 = N - 2*N3 ia = -180. - (90.-angle2) ix = triangle[2*N3-1,1] iy = triangle[2*N3-1,0] angle3 = np.arctan2(pi[0]-iy, pi[1]-ix) dx = (pi[1]-ix)/(N32-1) dy = (pi[0]-iy)/(N32-1) for ip in range(N32): triangle[2*N3+ip,:] = [iy+dy*ip, ix+dx*ip] return triangle def p_cross_width(larm=5., width=1., Narms=4, N=200): """ Function to draw a cross with arms with a given width and an angle larm: legnth of the arms (5., default) width: width of the arms (1., default) Narms: Number of arms (4, default) N: number of points to us (200, default) """ fname = 'p_cross_width' Narm = int((N-Narms)/Narms) larm2 = larm/2. width2 = width/2. cross = np.ones((N,2), dtype=np.float)*gen.fillValueF da = np.pi/Narms N1 = int(Narm*3./8.) N2 = int((Narm - 2*N1)/2.) N21 = Narm - 2*N1 - N2 if N2 < 3: print errormsg print ' ' + fname + ": too few points for ", Narms, " arms !!" print " increase number 'N' at least up to '", 25*Narms quit(-1) crosssecs = [] crossdic = {} Npot = int(np.log10(Narms))+1 iip = 0 for iarm in range(Narms-1): a = da*iarm iip0 = iip # bottom coordinate bx = larm*np.cos(a+np.pi) by = larm*np.sin(a+np.pi) # upper coordinate ux = larm*np.cos(a) uy = larm*np.sin(a) rela = a+np.pi*3./2. # SW-NW ix = bx + width2*np.cos(rela) iy = by + width2*np.sin(rela) ex = ux + width2*np.cos(rela) ey = uy + width2*np.sin(rela) dx = (ex-ix)/(N1-1) dy = (ey-iy)/(N1-1) for ip in range(N1): cross[iip+ip,:] = [iy+dy*ip,ix+dx*ip] iip = iip + N1 # NW-NE ix = ex + 0. iy = ey + 0. ex = ux - width2*np.cos(rela) ey = uy - width2*np.sin(rela) dx = (ex-ix)/(N2-1) dy = (ey-iy)/(N2-1) for ip in range(N2): cross[iip+ip,:] = [iy+dy*ip,ix+dx*ip] iip = iip + N2 # NW-SW ix = ex + 0. iy = ey + 0. ex = bx - width2*np.cos(rela) ey = by - width2*np.sin(rela) dx = (ex-ix)/(N1-1) dy = (ey-iy)/(N1-1) for ip in range(N1): cross[iip+ip,:] = [iy+dy*ip,ix+dx*ip] iip = iip + N1 # SW-SE ix = ex + 0. iy = ey + 0. ex = bx + width2*np.cos(rela) ey = by + width2*np.sin(rela) dx = (ex-ix)/(N21-1) dy = (ey-iy)/(N21-1) for ip in range(N21): cross[iip+ip,:] = [iy+dy*ip,ix+dx*ip] iip = iip + N21 + 1 iarmS = str(iarm).zfill(Npot) crosssecs.append(iarmS) crossdic[iarmS] = [cross[iip0:iip0+iip-1], '-', 'k', '1.'] iip0 = iip Narm = N - Narm*(Narms-1) - Narms N1 = int(Narm*3./8.) N2 = int((Narm - 2*N1)/2.) N21 = Narm - 2*N1 - N2 iarm = Narms-1 a = da*iarm # bottom coordinate bx = larm*np.cos(a+np.pi) by = larm*np.sin(a+np.pi) # upper coordinate ux = larm*np.cos(a) uy = larm*np.sin(a) rela = a+np.pi*3./2. # SW-NW ix = bx + width2*np.cos(rela) iy = by + width2*np.sin(rela) ex = ux + width2*np.cos(rela) ey = uy + width2*np.sin(rela) dx = (ex-ix)/(N1-1) dy = (ey-iy)/(N1-1) for ip in range(N1): cross[iip+ip,:] = [iy+dy*ip,ix+dx*ip] iip = iip + N1 # NW-NE ix = ex + 0. iy = ey + 0. ex = ux - width2*np.cos(rela) ey = uy - width2*np.sin(rela) dx = (ex-ix)/(N2-1) dy = (ey-iy)/(N2-1) for ip in range(N2): cross[iip+ip,:] = [iy+dy*ip,ix+dx*ip] iip = iip + N2 # NW-SW ix = ex + 0. iy = ey + 0. ex = bx - width2*np.cos(rela) ey = by - width2*np.sin(rela) dx = (ex-ix)/(N1-1) dy = (ey-iy)/(N1-1) for ip in range(N1): cross[iip+ip,:] = [iy+dy*ip,ix+dx*ip] iip = iip + N1 # SW-SE ix = ex + 0. iy = ey + 0. ex = bx + width2*np.cos(rela) ey = by + width2*np.sin(rela) dx = (ex-ix)/(N21-1) dy = (ey-iy)/(N21-1) for ip in range(N21): cross[iip+ip,:] = [iy+dy*ip,ix+dx*ip] iip = iip + N21 iarmS = str(iarm).zfill(Npot) crosssecs.append(iarmS) crossdic[iarmS] = [cross[iip0:iip0+iip-1], '-', 'k', '1.'] cross = ma.masked_equal(cross, gen.fillValueF) return cross, crosssecs, crossdic ####### ####### ##### #### ### ## # # Plotting def plot_sphere(iazm=-60., iele=30., dist=10., Npts=100, radii=10, \ drwsfc=[True,True], colsfc=['#AAAAAA','#646464'], \ drwxline = True, linex=[':','b',2.], drwyline = True, liney=[':','r',2.], \ drwzline = True, linez=['-.','g',2.], drwxcline=[True,True], \ linexc=[['-','#646400',1.],['--','#646400',1.]], \ drwequator=[True,True], lineeq=[['-','#AA00AA',1.],['--','#AA00AA',1.]], \ drwgreeenwhich=[True,True], linegw=[['-','k',1.],['--','k',1.]]): """ Function to plot an sphere and determine which standard lines will be also drawn iazm: azimut of the camera form the sphere iele: elevation of the camera form the sphere dist: distance of the camera form the sphere Npts: Resolution for the sphere radii: radius of the sphere drwsfc: whether 'up' and 'down' portions of the sphere should be drawn colsfc: colors of the surface of the sphere portions ['up', 'down'] drwxline: whether x-axis line should be drawn linex: properties of the x-axis line ['type', 'color', 'wdith'] drwyline: whether y-axis line should be drawn liney: properties of the y-axis line ['type', 'color', 'wdith'] drwzline: whether z-axis line should be drawn linez: properties of the z-axis line ['type', 'color', 'wdith'] drwequator: whether 'front' and 'back' portions of the Equator should be drawn lineeq: properties of the lines 'front' and 'back' of the Equator drwgreeenwhich: whether 'front', 'back' portions of Greenqhich should be drawn linegw: properties of the lines 'front' and 'back' Greenwhich drwxcline: whether 'front', 'back' 90 line (lon=90., lon=270.) should be drawn linexc: properties of the lines 'front' and 'back' for the 90 line """ fname = 'plot_sphere' iazmrad = iazm*np.pi/180. ielerad = iele*np.pi/180. # 3D surface Sphere sfcsphereu, sfcsphered = surface_sphere(radii,Npts) # greenwhich if iazmrad > np.pi/2. and iazmrad < 3.*np.pi/2.: ia=np.pi-ielerad else: ia=0.-ielerad ea=ia+np.pi da = (ea-ia)/(Npts-1) beta = np.arange(ia,ea+da,da)[0:Npts] alpha = np.zeros((Npts), dtype=np.float) greenwhichc = spheric_line(radii,alpha,beta) ia=ea+0. ea=ia+np.pi da = (ea-ia)/(Npts-1) beta = np.arange(ia,ea+da,da)[0:Npts] greenwhichd = spheric_line(radii,alpha,beta) # Equator ia=np.pi-iazmrad/2. ea=ia+np.pi da = (ea-ia)/(Npts-1) alpha = np.arange(ia,ea+da,da)[0:Npts] beta = np.zeros((Npts), dtype=np.float) equatorc = spheric_line(radii,alpha,beta) ia=ea+0. ea=ia+np.pi da = (ea-ia)/(Npts-1) alpha = np.arange(ia,ea+da,da)[0:Npts] equatord = spheric_line(radii,alpha,beta) # 90 line if iazmrad > np.pi and iazmrad < 2.*np.pi: ia=3.*np.pi/2. + ielerad else: ia=np.pi/2. - ielerad if ielerad < 0.: ia = ia + np.pi ea=ia+np.pi da = (ea-ia)/(Npts-1) beta = np.arange(ia,ea+da,da)[0:Npts] alpha = np.ones((Npts), dtype=np.float)*np.pi/2. xclinec = spheric_line(radii,alpha,beta) ia=ea+0. ea=ia+np.pi da = (ea-ia)/(Npts-1) beta = np.arange(ia,ea+da,da)[0:Npts] xclined = spheric_line(radii,alpha,beta) # x line xline = np.zeros((2,3), dtype=np.float) xline[0,:] = position_sphere(radii, 0., 0.) xline[1,:] = position_sphere(radii, np.pi, 0.) # y line yline = np.zeros((2,3), dtype=np.float) yline[0,:] = position_sphere(radii, np.pi/2., 0.) yline[1,:] = position_sphere(radii, 3*np.pi/2., 0.) # z line zline = np.zeros((2,3), dtype=np.float) zline[0,:] = position_sphere(radii, 0., np.pi/2.) zline[1,:] = position_sphere(radii, 0., -np.pi/2.) fig = plt.figure() ax = fig.gca(projection='3d') # Sphere surface if drwsfc[0]: ax.plot_surface(sfcsphereu[0,:,:], sfcsphereu[1,:,:], sfcsphereu[2,:,:], \ color=colsfc[0]) if drwsfc[1]: ax.plot_surface(sfcsphered[0,:,:], sfcsphered[1,:,:], sfcsphered[2,:,:], \ color=colsfc[1]) # greenwhich linev = linegw[0] if drwgreeenwhich[0]: ax.plot(greenwhichc[:,0], greenwhichc[:,1], greenwhichc[:,2], linev[0], \ color=linev[1], linewidth=linev[2], label='Greenwich') linev = linegw[1] if drwgreeenwhich[1]: ax.plot(greenwhichd[:,0], greenwhichd[:,1], greenwhichd[:,2], linev[0], \ color=linev[1], linewidth=linev[2]) # Equator linev = lineeq[0] if drwequator[0]: ax.plot(equatorc[:,0], equatorc[:,1], equatorc[:,2], linev[0], \ color=linev[1], linewidth=linev[2], label='Equator') linev = lineeq[1] if drwequator[1]: ax.plot(equatord[:,0], equatord[:,1], equatord[:,2], linev[0], \ color=linev[1], linewidth=linev[2]) # 90line linev = linexc[0] if drwxcline[0]: ax.plot(xclinec[:,0], xclinec[:,1], xclinec[:,2], linev[0], color=linev[1], \ linewidth=linev[2], label='90-line') linev = linexc[1] if drwxcline[1]: ax.plot(xclined[:,0], xclined[:,1], xclined[:,2], linev[0], color=linev[1], \ linewidth=linev[2]) # x line linev = linex if drwxline: ax.plot([xline[0,0],xline[1,0]], [xline[0,1],xline[1,1]], \ [xline[0,2],xline[1,2]], linev[0], color=linev[1], linewidth=linev[2], \ label='xline') # y line linev = liney if drwyline: ax.plot([yline[0,0],yline[1,0]], [yline[0,1],yline[1,1]], \ [yline[0,2],yline[1,2]], linev[0], color=linev[1], linewidth=linev[2], \ label='yline') # z line linev = linez if drwzline: ax.plot([zline[0,0],zline[1,0]], [zline[0,1],zline[1,1]], \ [zline[0,2],zline[1,2]], linev[0], color=linev[1], linewidth=linev[2], \ label='zline') plt.legend() return fig, ax def draw_secs(objdic): """ Function to draw an object according to its dictionary objdic: dictionary with the parts to draw [polygon, ltype, lcol, lw] """ fname = 'draw_secs' for secn in objdic.keys(): secv = objdic[secn] poly = secv[0] lt = secv[1] lc = secv[2] lw = secv[3] plt.plot(poly[:,1], poly[:,0], lt, color=lc, linewidth=lw) return def paint_filled(objdic, fillsecs): """ Function to draw an object filling given sections objdic: dictionary of the object filesecs: list of sections to be filled """ fname = 'paint_filled' Nsecs = len(fillsecs) for secn in fillsecs: secvals=objdic[secn] pvals = secvals[0] fillsecs = [] Nvals = pvals.shape[0] # re-sectionning to plot without masked values for ip in range(Nvals-1): if type(pvals[ip][0]) == type(gen.mamat[1]): fillsecs.append(ip) Nsecs = len(fillsecs) iisc = 0 for isc in range(Nsecs): plt.fill(pvals[iisc:fillsecs[isc],1], pvals[iisc:fillsecs[isc],0], \ color=secvals[2]) iisc = fillsecs[isc]+1 plt.fill(pvals[iisc:Nvals-1,1], pvals[iisc:Nvals-1,0], color=secvals[2]) return