Changeset 24254
- Timestamp:
- 10/18/19 00:13:27 (5 years ago)
- Location:
- issm/trunk-jpl/src/m
- Files:
-
- 19 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/m/coordsystems/ll2xy.m
r22339 r24254 1 function [x,y,scale_factor] = ll2xy(lat,lon,sgn,central_meridian,standard_parallel) 1 function [x,y,scale_factor] = ll2xy(lat,lon,sgn,central_meridian,standard_parallel) 2 2 %LL2XY - converts lat long to polar stereographic 3 3 % 4 % Converts from geodetic latitude and longitude to Polar 4 % Converts from geodetic latitude and longitude to Polar 5 5 % Stereographic (X,Y) coordinates for the polar regions. 6 6 % Optional scale factor provides the scaling factor needed to correct projection error … … 13 13 % [x,y] = ll2xy(lat,lon,sgn,central_meridian,standard_parallel) 14 14 % 15 % - sgn = Sign of latitude +1 : north latitude (default is mer=45 lat=70)16 % -1 : south latitude (default is mer=0 lat=71)15 % - sgn = Sign of latitude 1 : north latitude (default is mer=45 lat=70) 16 % -1 : south latitude (default is mer=0 lat=71) 17 17 18 18 %Get central_meridian and standard_parallel depending on hemisphere … … 55 55 T = tan(pi/4-latitude/2) ./ ((1-ex*sin(latitude))./(1+ex*sin(latitude))).^(ex/2); 56 56 57 if (90 - slat) < 1.e-5 57 if (90 - slat) < 1.e-5 58 58 rho = 2.*re*T/sqrt((1.+ex)^(1.+ex)*(1.-ex)^(1.-ex)); 59 59 else -
issm/trunk-jpl/src/m/coordsystems/ll2xy.py
r24213 r24254 14 14 x, y = ll2xy(lat, lon, sgn, central_meridian, standard_parallel) 15 15 16 - sgn = Sign of latitude +1 : north latitude (default is mer = 45 lat = 70)17 -1 : south latitude (default is mer = 0 lat = 71)16 - sgn = Sign of latitude 1 : north latitude (default is mer = 45 lat = 70) 17 -1 : south latitude (default is mer = 0 lat = 71) 18 18 ''' 19 19 20 assert sgn == 1 or sgn == - 1, 'error: sgn should be either + 1 or -1'20 assert sgn == 1 or sgn == -1, 'error: sgn should be either 1 or -1' 21 21 22 22 #Get central_meridian and standard_parallel depending on hemisphere … … 53 53 rho = re * mc * T / tc 54 54 55 y = - 55 y = -rho * sgn * np.cos(sgn * longitude) 56 56 x = rho * sgn * np.sin(sgn * longitude) 57 57 -
issm/trunk-jpl/src/m/coordsystems/xy2ll.m
r22493 r24254 14 14 % [lat,lon] = xy2ll(x,y,sgn,central_meridian,standard_parallel); 15 15 % 16 % - sgn = Sign of latitude +1 : north latitude (default is mer=45 lat=70)17 % -1 : south latitude (default is mer=0 lat=71)16 % - sgn = Sign of latitude 1 : north latitude (default is mer=45 lat=70) 17 % -1 : south latitude (default is mer=0 lat=71) 18 18 19 19 %Get central_meridian and standard_parallel depending on hemisphere … … 29 29 disp('Warning: expecting coordinates in polar stereographic (Std Latitude: 71ºS Meridian: 0º)'); 30 30 else 31 error('Sign should be either +1 or -1');31 error('Sign should be either 1 or -1'); 32 32 end 33 33 else … … 78 78 lon = lon * 180. / pi; 79 79 lat = lat * 180. / pi; 80 lon = lon - delta; 80 lon = lon - delta; 81 81 if nargout==3, 82 82 m=((1+sin(abs(slat)*pi/180))*ones(length(lat),1)./(1+sin(abs(lat)*pi/180))); -
issm/trunk-jpl/src/m/coordsystems/xy2ll.py
r24213 r24254 16 16 [lat, lon] = xy2ll(x, y, sgn, central_meridian, standard_parallel) 17 17 18 - sgn = Sign of latitude +1 : north latitude (default is mer = 45 lat = 70)19 -1 : south latitude (default is mer = 0 lat = 71)18 - sgn = Sign of latitude 1 : north latitude (default is mer = 45 lat = 70) 19 -1 : south latitude (default is mer = 0 lat = 71) 20 20 ''' 21 21 … … 29 29 slat = 70. 30 30 print(' xy2ll: creating coordinates in north polar stereographic (Std Latitude: 70degN Meridian: 45deg)') 31 elif sgn == - 31 elif sgn == -1: 32 32 delta = 0. 33 33 slat = 71. 34 34 print(' xy2ll: creating coordinates in south polar stereographic (Std Latitude: 71degS Meridian: 0deg)') 35 35 else: 36 raise ValueError('sgn should be either + 1 or -1')36 raise ValueError('sgn should be either 1 or -1') 37 37 else: 38 38 raise Exception('bad usage: type "help(xy2ll)" for details') … … 67 67 68 68 lat = sgn * lat 69 lon = np.arctan2(sgn * x, - 69 lon = np.arctan2(sgn * x, -sgn * y) 70 70 lon = sgn * lon 71 71 -
issm/trunk-jpl/src/m/dev/devpath.py
r24241 r24254 11 11 raise NameError('"ISSM_DIR" environment variable is empty! You should define ISSM_DIR in your .cshrc or .bashrc!') 12 12 13 #Go through src / m and append any directory that contains a *.py file to PATH13 #Go through src/m and append any directory that contains a *.py file to PATH 14 14 for root, dirs, files in os.walk(ISSM_DIR + '/src/m'): 15 15 if '.svn' in dirs: 16 16 dirs.remove('.svn') 17 17 for file in files: 18 if file.find(".py") != - 19 if file.find(".pyc") == - 18 if file.find(".py") != -1: 19 if file.find(".pyc") == -1: 20 20 if root not in sys.path: 21 21 sys.path.append(root) … … 35 35 sys.path.append(jpl_path) 36 36 else: 37 warnings.warn('cluster settings should be in, {} / usr /{}'.format(JPL_SVN, USERNAME))37 warnings.warn('cluster settings should be in, {}/usr/{}'.format(JPL_SVN, USERNAME)) 38 38 39 39 from runme import runme #first because plotmodel may fail -
issm/trunk-jpl/src/m/exp/expcoarsen.py
r24213 r24254 55 55 yi = np.linspace(y[j], y[j + 1], division) 56 56 57 x = np.hstack((x[0:j + 1], xi[1: -1], x[j + 1:]))58 y = np.hstack((y[0:j + 1], yi[1: -1], y[j + 1:]))57 x = np.hstack((x[0:j + 1], xi[1:-1], x[j + 1:])) 58 y = np.hstack((y[0:j + 1], yi[1:-1], y[j + 1:])) 59 59 60 60 #update current point -
issm/trunk-jpl/src/m/exp/expread.py
r24213 r24254 50 50 51 51 if len(A[1]) > 5: 52 contour['name'] = A[1][5: -1]52 contour['name'] = A[1][5:-1] 53 53 else: 54 54 contour['name'] = '' -
issm/trunk-jpl/src/m/exp/flowlines.py
r24213 r24254 91 91 flowindex = 0 92 92 elif flowdirection == 'downstream': 93 flowindex = - 93 flowindex = -1 94 94 95 95 while not all(done): … … 112 112 #velocity of the current triangle and norm it 113 113 if flowdirection == 'upstream': 114 ut = - 115 vt = - 114 ut = -u[tria] 115 vt = -v[tria] 116 116 if flowdirection == 'downstream': 117 117 ut = u[tria] -
issm/trunk-jpl/src/m/geometry/NowickiProfile.py
r24213 r24254 39 39 #downstream of the GL 40 40 for i in range(int(np.size(x) / 2), int(np.size(x))): 41 h[i] = (x[i] / (4. * (delta + 1) * q) + hg**(- 2))**(-0.5) # ice thickness for ice shelf from (3.1)41 h[i] = (x[i] / (4. * (delta + 1) * q) + hg**(-2))**(-0.5) # ice thickness for ice shelf from (3.1) 42 42 b[i] = sea - h[i] * (1. / (1 + delta)) 43 43 -
issm/trunk-jpl/src/m/geometry/SegIntersect.py
r24213 r24254 48 48 49 49 #if colinear 50 if test1 * test2 == 0 and test3 * test4 == 0 and np.linalg.det(np.hstack((n1.reshape((- 1, )), n2.reshape(-1, )))) == 0:50 if test1 * test2 == 0 and test3 * test4 == 0 and np.linalg.det(np.hstack((n1.reshape((-1, )), n2.reshape(-1, )))) == 0: 51 51 52 52 #projection on the axis O1O2 … … 57 57 O1D = np.dot(O2O1, O1D) 58 58 59 #test if one point is included in the other segment (- 59 #test if one point is included in the other segment (-> bval = 1) 60 60 if (O1C - O1A) * (O1D - O1A) < 0: 61 61 bval = 1 … … 71 71 return bval 72 72 73 #test if the 2 segments have the same middle (- 73 #test if the 2 segments have the same middle (-> bval = 1) 74 74 if O2O1 == 0: 75 75 bval = 1 -
issm/trunk-jpl/src/m/geometry/slope.py
r24213 r24254 34 34 35 35 summation = np.array([[1], [1], [1]]) 36 sx = np.dot(surf[index - 1, 0] * alpha, summation).reshape(- 37 sy = np.dot(surf[index - 1, 0] * beta, summation).reshape(- 36 sx = np.dot(surf[index - 1, 0] * alpha, summation).reshape(-1, ) 37 sy = np.dot(surf[index - 1, 0] * beta, summation).reshape(-1, ) 38 38 39 39 s = np.sqrt(sx**2 + sy**2) -
issm/trunk-jpl/src/m/mech/backstressfrominversion.py
r24213 r24254 62 62 T = 0.5 * md.materials.rho_ice * md.constants.g * (1 - md.materials.rho_ice / md.materials.rho_water) * md.geometry.thickness 63 63 n = averaging(md, md.materials.rheology_n, 0) 64 Bi = md.results.StressbalanceSolution.MaterialsRheologyBbar.reshape(- 64 Bi = md.results.StressbalanceSolution.MaterialsRheologyBbar.reshape(-1, ) 65 65 66 66 a0, b0, theta0, ex0 = thomasparams(md, eq='Thomas', smoothing=smoothing, coordsys=coordsys) -
issm/trunk-jpl/src/m/mech/damagefrominversion.py
r24213 r24254 26 26 raise Exception('Warning: the model has some non - SSA elements. These will be treated like SSA elements') 27 27 if np.ndim(md.results.StressbalanceSolution.MaterialsRheologyBbar) == 2: 28 Bi = md.results.StressbalanceSolution.MaterialsRheologyBbar.reshape(- 28 Bi = md.results.StressbalanceSolution.MaterialsRheologyBbar.reshape(-1, ) 29 29 else: 30 30 Bi = md.results.StressbalanceSolution.MaterialsRheologyBbar 31 31 if np.ndim(md.materials.rheology_B) == 2: 32 BT = md.materials.rheology_B.reshape(- 32 BT = md.materials.rheology_B.reshape(-1, ) 33 33 else: 34 34 BT = md.materials.rheology_B -
issm/trunk-jpl/src/m/mech/mechanicalproperties.py
r24213 r24254 60 60 vxlist = vx[index - 1] / md.constants.yts 61 61 vylist = vy[index - 1] / md.constants.yts 62 ux = np.dot((vxlist * alpha), summation).reshape(- 63 uy = np.dot((vxlist * beta), summation).reshape(- 64 vx = np.dot((vylist * alpha), summation).reshape(- 65 vy = np.dot((vylist * beta), summation).reshape(- 62 ux = np.dot((vxlist * alpha), summation).reshape(-1, ) 63 uy = np.dot((vxlist * beta), summation).reshape(-1, ) 64 vx = np.dot((vylist * alpha), summation).reshape(-1, ) 65 vy = np.dot((vylist * beta), summation).reshape(-1, ) 66 66 uyvx = (vx + uy) / 2. 67 67 #clear vxlist vylist … … 69 69 #compute viscosity 70 70 nu = np.zeros((numberofelements, )) 71 B_bar = np.dot(md.materials.rheology_B[index - 1], summation / 3.).reshape(- 72 power = ((md.materials.rheology_n - 1.) / (2. * md.materials.rheology_n)).reshape(- 73 second_inv = (ux**2. + vy**2. + ((uy + vx)**2.) / 4. + ux * vy).reshape(- 71 B_bar = np.dot(md.materials.rheology_B[index - 1], summation / 3.).reshape(-1, ) 72 power = ((md.materials.rheology_n - 1.) / (2. * md.materials.rheology_n)).reshape(-1, ) 73 second_inv = (ux**2. + vy**2. + ((uy + vx)**2.) / 4. + ux * vy).reshape(-1, ) 74 74 75 75 #some corrections … … 86 86 elif 'matdamageice' in md.materials.__module__ and damage is not None: 87 87 print('computing damage-dependent properties!') 88 Zinv = np.dot(1 - damage[index - 1], summation / 3.).reshape(- 88 Zinv = np.dot(1 - damage[index - 1], summation / 3.).reshape(-1, ) 89 89 location = np.nonzero(second_inv) 90 90 nu[location] = Zinv[location] * B_bar[location] / np.power(second_inv[location], power[location]) … … 109 109 #eigenvalues and vectors for stress 110 110 value, directions = np.linalg.eig(stress) 111 idx = value.argsort()[:: -1] # sort in descending algebraic (not absolute) order111 idx = value.argsort()[::-1] # sort in descending algebraic (not absolute) order 112 112 value = value[idx] 113 113 directions = directions[:, idx] … … 117 117 #eigenvalues and vectors for strain 118 118 value, directions = np.linalg.eig(strain) 119 idx = value.argsort()[:: -1] # sort in descending order119 idx = value.argsort()[::-1] # sort in descending order 120 120 value = value[idx] 121 121 directions = directions[:, idx] -
issm/trunk-jpl/src/m/mech/robintemperature.py
r24213 r24254 16 16 17 17 Parameters (SI units): 18 - heatflux Geothermal heat flux (W m^ -2)19 - accumrate Surface accumulation rate (m s^ -1 ice equivalent)18 - heatflux Geothermal heat flux (W m^-2) 19 - accumrate Surface accumulation rate (m s^-1 ice equivalent) 20 20 - thickness Ice thickness (m) 21 21 - surftemp Surface temperature (K) … … 30 30 31 31 # some constants (from Holland and Jenkins, 1999) 32 alphaT = 1.14e-6 # thermal diffusivity (m^2 s^ -1)33 c = 2009. # specific heat capacity (J kg^ - 1 K^ -1)34 rho = 917. # ice density (kg m^ -3)32 alphaT = 1.14e-6 # thermal diffusivity (m^2 s^-1) 33 c = 2009. # specific heat capacity (J kg^-1 K^-1) 34 rho = 917. # ice density (kg m^-3) 35 35 36 36 #create vertical coordinate variable 37 37 zstar = np.sqrt(2. * alphaT * thickness / accumrate) 38 38 39 tprofile = surftemp + np.sqrt(2. * thickness * np.pi / accumrate / alphaT) * (- 39 tprofile = surftemp + np.sqrt(2. * thickness * np.pi / accumrate / alphaT) * (-heatflux) / 2. / rho / c * (erf(z / zstar) - erf(thickness / zstar)) 40 40 41 41 return tprofile -
issm/trunk-jpl/src/m/mech/thomasparams.py
r24213 r24254 72 72 73 73 # average element strain rates onto vertices 74 e1 = averaging(md, md.results.strainrate.principalvalue1, smoothing) / md.constants.yts # convert to s^ -174 e1 = averaging(md, md.results.strainrate.principalvalue1, smoothing) / md.constants.yts # convert to s^-1 75 75 e2 = averaging(md, md.results.strainrate.principalvalue2, smoothing) / md.constants.yts 76 76 exx = averaging(md, md.results.strainrate.xx, smoothing) / md.constants.yts … … 81 81 pos = np.nonzero(e1 == 0) 82 82 if np.any(pos == 1): 83 print('WARNING: first principal strain rate equal to zero. Value set to 1e-13 s^ -1')83 print('WARNING: first principal strain rate equal to zero. Value set to 1e-13 s^-1') 84 84 e1[pos] = 1.e-13 85 85 pos = np.nonzero(e2 == 0) 86 86 if np.any(pos == 1): 87 print('WARNING: second principal strain rate equal to zero. Value set to 1e-13 s^ -1')87 print('WARNING: second principal strain rate equal to zero. Value set to 1e-13 s^-1') 88 88 e2[pos] = 1.e-13 89 89 … … 119 119 ex = 0.5 * (exx + eyy) + 0.5 * (exx - eyy) * np.cos(2. * velangle) + exy * np.sin(2. * velangle) 120 120 ey = exx + eyy - ex # trace of strain rate tensor is invariant 121 exy = - 121 exy = -0.5 * (exx - eyy) * np.sin(2. * velangle) + exy * np.cos(2. * velangle) 122 122 a = ey / ex 123 123 b = exy / ex … … 130 130 pos = np.nonzero(np.abs((np.abs(a) - 2.)) < 1.e-3) 131 131 if len(pos) > 0: 132 print(('Warning: ', len(pos), ' vertices have alpha within 1e-3 of - 132 print(('Warning: ', len(pos), ' vertices have alpha within 1e-3 of -2')) 133 133 a[pos] = -2 + 1e-3 134 134 -
issm/trunk-jpl/src/m/partition/adjacency.py
r24213 r24254 19 19 20 20 md.qmu.adjacency = m.sparse(indi, indj, values, md.mesh.numberofvertices, md.mesh.numberofvertices) 21 md.qmu.adjacency = np.logical_or(md.qmu.adjacency, md.qmu.adjacency.T).astype(float) #change to reshape(- 21 md.qmu.adjacency = np.logical_or(md.qmu.adjacency, md.qmu.adjacency.T).astype(float) #change to reshape(-1, 1) if needed 22 22 23 23 #now, build vwgt: … … 27 27 md.mesh.vertexconnectivity = NodeConnectivity(md.mesh.elements, md.mesh.numberofvertices) 28 28 29 connectivity = md.mesh.vertexconnectivity[0][:, 0: -1]29 connectivity = md.mesh.vertexconnectivity[0][:, 0:-1] 30 30 pos = np.where(connectivity) 31 31 connectivity[pos] = areas[connectivity[pos] - 1] / 3. -
issm/trunk-jpl/src/m/partition/partitioner.py
r24213 r24254 76 76 weights = [] 77 77 78 method = method.reshape(- 78 method = method.reshape(-1, 1) # transpose to 1x10 instead of 10 79 79 # partition into nparts 80 80 if isinstance(md.mesh, mesh2d): … … 120 120 part = project3d(md, 'vector', np.squeeze(part), 'type', 'node') 121 121 122 part = part.reshape(- 122 part = part.reshape(-1, 1) 123 123 124 124 if vectortype == 'element': -
issm/trunk-jpl/src/m/qmu/helpers.py
r24213 r24254 24 24 Example uses: 25 25 26 x = Lstruct(1, 2, 3, 4) - 26 x = Lstruct(1, 2, 3, 4) -> [1, 2, 3, 4] 27 27 x.a = 'hello' 28 len(x) - 28 len(x) -> 4 29 29 x.append(5) 30 len(x) - 31 x[2] - 32 x.a - 33 print x - 34 x.__dict__ - > {'a':'hello'}30 len(x) -> 5 31 x[2] -> 3 32 x.a -> 'hello' 33 print x -> [1, 2, 3, 4, 5] 34 x.__dict__ -> {'a':'hello'} 35 35 x.b = [6, 7, 8, 9] 36 x.b[-1] - 37 len(x.b) - 36 x.b[-1] -> 9 37 len(x.b) -> 4 38 38 39 39 Other valid constructors: 40 x = Lstruct(1, 2, 3, a = 'hello') - > x.a - > 'hello', x -> [1, 2, 3]40 x = Lstruct(1, 2, 3, a = 'hello') -> x.a -> 'hello', x -> [1, 2, 3] 41 41 x = Lstruct(1, 2, 3)(a = 'hello') 42 42 x = Lstruct([1, 2, 3], x = 'hello') 43 43 x = Lstruct((1, 2, 3), a = 'hello') 44 44 45 Credit: https: / / github.com / Vectorized / Python -Attribute-List45 Credit: https://github.com/Vectorized/Python-Attribute-List 46 46 ''' 47 47 … … 63 63 class OrderedStruct(object): 64 64 ''' 65 A form of dictionary -like structure that maintains the65 A form of dictionary-like structure that maintains the 66 66 ordering in which its fields / attributes and their 67 67 corresponding values were added. … … 70 70 can be used as an "ordered struct / class" giving 71 71 it much more flexibility in practice. It is 72 also easier to work with fixed valued keys in -code.72 also easier to work with fixed valued keys in-code. 73 73 74 74 Eg: … … 79 79 x.y = 5 80 80 OR 81 x['y'] = 5 # supports OrderedDict -style usage82 83 Supports: len(x), str(x), for -loop iteration.81 x['y'] = 5 # supports OrderedDict-style usage 82 83 Supports: len(x), str(x), for-loop iteration. 84 84 Has methods: x.keys(), x.values(), x.items(), x.iterkeys() 85 85 … … 94 94 # in the same order as the inputs 95 95 96 x.keys() - 97 x.values() - 98 x.items() - 99 x.__dict__ - 100 vars(x) - 101 102 x.y - 103 x['y'] - 104 x.z - 105 x['z'] - 96 x.keys() -> ['y', 'z'] 97 x.values() -> [5, 6] 98 x.items() -> [('y', 6), ('z', 6)] 99 x.__dict__ -> [('y', 6), ('z', 6)] 100 vars(x) -> [('y', 6), ('z', 6)] 101 102 x.y -> 5 103 x['y'] -> 5 104 x.z -> 6 105 x['z'] -> 6 106 106 107 107 for i in x: # same as x.items() 108 108 print i 109 - 109 -> 110 110 ('x', 5) 111 111 ('y', 6) … … 197 197 # same thing but call deepcopy recursively 198 198 # technically not how it should be done, 199 # (see https: / / docs.python.org / 2 / library /copy.html #copy.deepcopy )199 # (see https://docs.python.org/2/library/copy.html #copy.deepcopy ) 200 200 # but will generally work in this case 201 201 newInstance = type(self)() … … 227 227 def isempty(x): 228 228 ''' 229 returns true if object is + - 230 array /matrix composed only of such components (includes mixtures of "empty" types)'''229 returns true if object is + -infinity, NaN, None, '', has length 0, or is an 230 array/matrix composed only of such components (includes mixtures of "empty" types)''' 231 231 232 232 if type(x) in [list, np.ndarray, tuple]: … … 247 247 if x is None: 248 248 return True 249 if type(x) == str and x.lower() in ['', 'nan', 'none', 'inf', 'infinity', ' - inf', ' -infinity']:249 if type(x) == str and x.lower() in ['', 'nan', 'none', 'inf', 'infinity', '-inf', '-infinity']: 250 250 return True 251 251 … … 280 280 def fileparts(x): 281 281 ''' 282 given: "path / path / ... /file_name.ext"282 given: "path/path/.../file_name.ext" 283 283 returns: [path, file_name, ext] (list of strings)''' 284 284 try: … … 300 300 301 301 fullfile(path, path, ... , file_name + ext) 302 returns: "path / path / ... /file_name.ext"303 304 with all arguments as strings with no " /"s302 returns: "path/path/.../file_name.ext" 303 304 with all arguments as strings with no "/"s 305 305 306 306 regarding extensions and the '.': … … 331 331 def empty_nd_list(shape, filler=0., as_numpy_ndarray=False): 332 332 ''' 333 returns a python list of the size /shape given (shape must be int or tuple)333 returns a python list of the size/shape given (shape must be int or tuple) 334 334 the list will be filled with the optional second argument 335 335 … … 338 338 as_numpy_ndarray will return the result as a numpy.ndarray and is False by default 339 339 340 Note: the filler must be either None / np.nan / float('NaN'), float /double, or int341 other numpy and float values such as + /- np.inf will also work340 Note: the filler must be either None/np.nan/float('NaN'), float/double, or int 341 other numpy and float values such as +/- np.inf will also work 342 342 343 343 use:
Note:
See TracChangeset
for help on using the changeset viewer.