Changeset 21255


Ignore:
Timestamp:
10/11/16 01:27:08 (8 years ago)
Author:
bdef
Message:

CHG: replacing npy by np for numpy abreviation

Location:
issm/trunk-jpl/src/py3
Files:
34 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/py3/boundaryconditions/PattynSMB.py

    r19895 r21255  
    11import os
    2 import numpy as npy
     2import numpy as np
    33def PattynSMB(md,Tf):
    44        """
    … …  
    2525
    2626        #Double check lat and long exist:
    27         if npy.any(npy.isnan(md.mesh.lat)):
     27        if np.any(np.isnan(md.mesh.lat)):
    2828                raise IOError('PattynSMB error message: md.mesh.lat field required')
    2929
    … …  
    3636        # Calculate summer temperature, Eqn (12)
    3737        # No melting at PIG in mean conditions - need about 6 degress Tf to start having a negative yearly SMB
    38         Tms = 16.81 - 0.00692*md.geometry.surface - 0.27937*npy.abs(md.mesh.lat) + Tf
     38        Tms = 16.81 - 0.00692*md.geometry.surface - 0.27937*np.abs(md.mesh.lat) + Tf
    3939        Tms= Tms[0]
    4040
    … …  
    4343
    4444        # Calculate Ablation, Eqn (10) (use for both Antarctica & Greenland), max melt is 10m/a
    45         ABLdot=0.*npy.ones(md.mesh.numberofvertices)
    46         pos=npy.nonzero(Tms>=0)
    47         ABLdot[pos]=npy.minimum(1.4*Tms[pos],10)
     45        ABLdot=0.*np.ones(md.mesh.numberofvertices)
     46        pos=np.nonzero(Tms>=0)
     47        ABLdot[pos]=np.minimum(1.4*Tms[pos],10)
    4848
    4949        smb=ACCdot-ABLdot

  • issm/trunk-jpl/src/py3/classes/outputdefinition.py

    r19895 r21255  
    44from checkfield import checkfield
    55from WriteData import WriteData
    6 import numpy as npy
     6import numpy as np
    77
    88class outputdefinition(object):
    … …  
    3636        def marshall(self,md,fid):    # {{{
    3737               
    38                 enums=npy.zeros(len(self.definitions),)
     38                enums=np.zeros(len(self.definitions),)
    3939               
    4040                for i in range(len(self.definitions)):
    … …  
    4444                        enums[i]=StringToEnum(classdefinition)[0]
    4545               
    46                 enums=npy.unique(enums);
     46                enums=np.unique(enums);
    4747               
    4848                WriteData(fid,'data',enums,'enum',OutputdefinitionListEnum(),'format','DoubleMat','mattype',1);

  • issm/trunk-jpl/src/py3/coordsystems/ll2xy.py

    r19895 r21255  
    1 import numpy as npy
     1import numpy as np
    22
    33def ll2xy(lat,lon,sgn=-1,central_meridian=0,standard_parallel=71):
    … …  
    3636        ex2 = .006693883
    3737        # Eccentricity of the Hughes ellipsoid
    38         ex = npy.sqrt(ex2)
     38        ex = np.sqrt(ex2)
    3939       
    40         latitude = npy.abs(lat) * npy.pi/180.
    41         longitude = (lon + delta) * npy.pi/180.
     40        latitude = np.abs(lat) * np.pi/180.
     41        longitude = (lon + delta) * np.pi/180.
    4242       
    4343        # compute X and Y in grid coordinates.
    44         T = npy.tan(npy.pi/4-latitude/2) / ((1-ex*npy.sin(latitude))/(1+ex*npy.sin(latitude)))**(ex/2)
     44        T = np.tan(np.pi/4-latitude/2) / ((1-ex*np.sin(latitude))/(1+ex*np.sin(latitude)))**(ex/2)
    4545       
    4646        if (90 - slat) <  1.e-5:
    47                 rho = 2.*re*T/npy.sqrt((1.+ex)**(1.+ex)*(1.-ex)**(1.-ex))
     47                rho = 2.*re*T/np.sqrt((1.+ex)**(1.+ex)*(1.-ex)**(1.-ex))
    4848        else:
    49                 sl  = slat*npy.pi/180.
    50                 tc  = npy.tan(npy.pi/4.-sl/2.)/((1.-ex*npy.sin(sl))/(1.+ex*npy.sin(sl)))**(ex/2.)
    51                 mc  = npy.cos(sl)/npy.sqrt(1.0-ex2*(npy.sin(sl)**2))
     49                sl  = slat*np.pi/180.
     50                tc  = np.tan(np.pi/4.-sl/2.)/((1.-ex*np.sin(sl))/(1.+ex*np.sin(sl)))**(ex/2.)
     51                mc  = np.cos(sl)/np.sqrt(1.0-ex2*(np.sin(sl)**2))
    5252                rho = re*mc*T/tc
    5353       
    54         y = -rho * sgn * npy.cos(sgn*longitude)
    55         x =  rho * sgn * npy.sin(sgn*longitude)
     54        y = -rho * sgn * np.cos(sgn*longitude)
     55        x =  rho * sgn * np.sin(sgn*longitude)
    5656
    57         cnt1=npy.nonzero(latitude>= npy.pi/2.)[0]
     57        cnt1=np.nonzero(latitude>= np.pi/2.)[0]
    5858       
    5959        if cnt1:

  • issm/trunk-jpl/src/py3/coordsystems/xy2ll.py

    r19895 r21255  
    1 import numpy as npy
     1import numpy as np
    22from math import pi
    33
    … …  
    3939        # if x,y passed as lists, convert to numpy arrays
    4040        if type(x) != "numpy.ndarray":
    41                 x=npy.array(x)
     41                x=np.array(x)
    4242        if type(y) != "numpy.ndarray":
    43                 y=npy.array(y)
     43                y=np.array(y)
    4444
    4545        ## Conversion constant from degrees to radians
    … …  
    5050        ex2 = .006693883
    5151        ## Eccentricity of the Hughes ellipsoid
    52         ex = npy.sqrt(ex2)
     52        ex = np.sqrt(ex2)
    5353       
    5454        sl = slat*pi/180.
    55         rho = npy.sqrt(x**2 + y**2)
    56         cm = npy.cos(sl) / npy.sqrt(1.0 - ex2 * (npy.sin(sl)**2))
    57         T = npy.tan((pi/4.0) - (sl/2.0)) / ((1.0 - ex*npy.sin(sl)) / (1.0 + ex*npy.sin(sl)))**(ex / 2.0)
     55        rho = np.sqrt(x**2 + y**2)
     56        cm = np.cos(sl) / np.sqrt(1.0 - ex2 * (np.sin(sl)**2))
     57        T = np.tan((pi/4.0) - (sl/2.0)) / ((1.0 - ex*np.sin(sl)) / (1.0 + ex*np.sin(sl)))**(ex / 2.0)
    5858       
    5959        if abs(slat-90.) < 1.e-5:
    60                 T = rho*npy.sqrt((1. + ex)**(1. + ex) * (1. - ex)**(1. - ex)) / 2. / re
     60                T = rho*np.sqrt((1. + ex)**(1. + ex) * (1. - ex)**(1. - ex)) / 2. / re
    6161        else:
    6262                T = rho * T / (re * cm)
    6363       
    64         chi = (pi / 2.0) - 2.0 * npy.arctan(T)
     64        chi = (pi / 2.0) - 2.0 * np.arctan(T)
    6565        lat = chi + ((ex2 / 2.0) + (5.0 * ex2**2.0 / 24.0) + (ex2**3.0 / 12.0)) * \
    66                 npy.sin(2 * chi) + ((7.0 * ex2**2.0 / 48.0) + (29.0 * ex2**3 / 240.0)) * \
    67                 npy.sin(4.0 * chi) + (7.0 * ex2**3.0 / 120.0) * npy.sin(6.0 * chi)
     66                np.sin(2 * chi) + ((7.0 * ex2**2.0 / 48.0) + (29.0 * ex2**3 / 240.0)) * \
     67                np.sin(4.0 * chi) + (7.0 * ex2**3.0 / 120.0) * np.sin(6.0 * chi)
    6868       
    6969        lat = sgn * lat
    70         lon = npy.arctan2(sgn * x,-sgn * y)
     70        lon = np.arctan2(sgn * x,-sgn * y)
    7171        lon = sgn * lon
    7272       
    73         res1 = npy.nonzero(rho <= 0.1)[0]
     73        res1 = np.nonzero(rho <= 0.1)[0]
    7474        if len(res1) > 0:
    7575                lat[res1] = 90. * sgn

  • issm/trunk-jpl/src/py3/exp/expcoarsen.py

    r19895 r21255  
    11import os.path
    2 import numpy as npy
     2import numpy as np
    33from collections import OrderedDict
    44from expread import expread
    … …  
    3434        for contour in  contours:
    3535               
    36                 numpoints=npy.size(contour['x'])
     36                numpoints=np.size(contour['x'])
    3737
    3838                j=0
    … …  
    4444
    4545                        #see whether we keep this point or not
    46                         distance=npy.sqrt((x[j]-x[j+1])**2+(y[j]-y[j+1])**2)
     46                        distance=np.sqrt((x[j]-x[j+1])**2+(y[j]-y[j+1])**2)
    4747                        if distance<resolution and j<numpoints-2:   #do not remove last point
    48                                 x=npy.delete(x,j+1,0)
    49                                 y=npy.delete(y,j+1,0)
     48                                x=np.delete(x,j+1,0)
     49                                y=np.delete(y,j+1,0)
    5050                                numpoints=numpoints-1
    5151                        else:
    52                                 division=int(npy.floor(distance/resolution)+1)
     52                                division=int(np.floor(distance/resolution)+1)
    5353                                if division>=2:
    54                                         xi=npy.linspace(x[j],x[j+1],division)
    55                                         yi=npy.linspace(y[j],y[j+1],division)
     54                                        xi=np.linspace(x[j],x[j+1],division)
     55                                        yi=np.linspace(y[j],y[j+1],division)
    5656                                       
    57                                         x=npy.hstack((x[0:j+1],xi[1:-1],x[j+1:]))
    58                                         y=npy.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:]))
    5959
    6060                                        #update current point
    … …  
    6565                                        j=j+1
    6666               
    67                 if npy.size(x)>1:
     67                if np.size(x)>1:
    6868                        #keep the (x,y) contour arond
    6969                        newcontour=OrderedDict()
    70                         newcontour['nods']=npy.size(x)
     70                        newcontour['nods']=np.size(x)
    7171                        newcontour['density']=contour['density']
    7272                        newcontour['x']=x

  • issm/trunk-jpl/src/py3/exp/expdisp.py

    r19895 r21255  
    11from expread import expread
    2 import numpy as npy
     2import numpy as np
    33
    44def expdisp(domainoutline,ax,linestyle='--k',linewidth=1,unitmultiplier=1.):

  • issm/trunk-jpl/src/py3/extrusion/DepthAverage.py

    r19895 r21255  
    1 import numpy as npy
     1import numpy as np
    22from project2d import project2d
    33
    … …  
    1919
    2020        # coerce to array in case float is passed
    21         if type(vector)!=npy.ndarray:
     21        if type(vector)!=np.ndarray:
    2222                print('coercing array')
    23                 vector=npy.array(value)
     23                vector=np.array(value)
    2424
    2525        vec2d=False
    … …  
    3030        #nods data
    3131        if vector.shape[0]==md.mesh.numberofvertices:
    32                 vector_average=npy.zeros(md.mesh.numberofvertices2d)
     32                vector_average=np.zeros(md.mesh.numberofvertices2d)
    3333                for i in range(1,md.mesh.numberoflayers):
    3434                        vector_average=vector_average+(project2d(md,vector,i)+project2d(md,vector,i+1))/2.*(project2d(md,md.mesh.z,i+1)-project2d(md,md.mesh.z,i))
    … …  
    3737        #element data
    3838        elif vector.shape[0]==md.mesh.numberofelements:
    39                 vector_average=npy.zeros(md.mesh.numberofelements2d)
     39                vector_average=np.zeros(md.mesh.numberofelements2d)
    4040                for i in range(1,md.mesh.numberoflayers):
    4141                        vector_average=vector_average+project2d(md,vector,i)*(project2d(md,md.mesh.z,i+1)-project2d(md,md.mesh.z,i))

  • issm/trunk-jpl/src/py3/extrusion/project2d.py

    r19895 r21255  
    1 import numpy as npy
     1import numpy as np
    22
    33def project2d(md3d,value,layer):
    … …  
    2525       
    2626        # coerce to array in case float is passed
    27         if type(value)!=npy.ndarray:
     27        if type(value)!=np.ndarray:
    2828                print('coercing array')
    29                 value=npy.array(value)
     29                value=np.array(value)
    3030
    3131        vec2d=False

  • issm/trunk-jpl/src/py3/geometry/slope.py

    r19895 r21255  
    1 import numpy as npy
     1import numpy as np
    22from GetNodalFunctionsCoeff import  GetNodalFunctionsCoeff
    33
    … …  
    3333        alpha,beta=GetNodalFunctionsCoeff(index,x,y)[0:2]
    3434
    35         summation=npy.array([[1],[1],[1]])
    36         sx=npy.dot(surf[index-1]*alpha,summation).reshape(-1,)
    37         sy=npy.dot(surf[index-1]*beta,summation).reshape(-1,)
     35        summation=np.array([[1],[1],[1]])
     36        sx=np.dot(surf[index-1]*alpha,summation).reshape(-1,)
     37        sy=np.dot(surf[index-1]*beta,summation).reshape(-1,)
    3838
    39         s=npy.sqrt(sx**2+sy**2)
     39        s=np.sqrt(sx**2+sy**2)
    4040
    4141        if md.mesh.dimension()==3:
    4242                sx=project3d(md,'vector',sx,'type','element')
    4343                sy=project3d(md,'vector',sy,'type','element')
    44                 s=npy.sqrt(sx**2+sy**2)
     44                s=np.sqrt(sx**2+sy**2)
    4545
    4646        return (sx,sy,s)

  • issm/trunk-jpl/src/py3/interp/SectionValues.py

    r19895 r21255  
    11import os
    22from expread import expread
    3 import numpy as npy
     3import numpy as np
    44from project2d import project2d
    55#from InterpFromMesh2d import InterpFromMesh2d
    … …  
    4141
    4242        #initialization
    43         X=npy.array([]) #X-coordinate
    44         Y=npy.array([]) #Y-coordinate
    45         S=npy.array([0.])  #curvilinear coordinate
     43        X=np.array([]) #X-coordinate
     44        Y=np.array([]) #Y-coordinate
     45        S=np.array([0.])  #curvilinear coordinate
    4646       
    4747        for i in range(nods-1):
    … …  
    5353                s_start=S[-1]
    5454       
    55                 length_segment=npy.sqrt((x_end-x_start)**2+(y_end-y_start)**2)
    56                 portion=npy.ceil(length_segment/res_h)
     55                length_segment=np.sqrt((x_end-x_start)**2+(y_end-y_start)**2)
     56                portion=np.ceil(length_segment/res_h)
    5757       
    58                 x_segment=npy.zeros(portion)
    59                 y_segment=npy.zeros(portion)
    60                 s_segment=npy.zeros(portion)
     58                x_segment=np.zeros(portion)
     59                y_segment=np.zeros(portion)
     60                s_segment=np.zeros(portion)
    6161
    6262                for j in range(int(portion)):
    … …  
    6666       
    6767                #plug into X and Y
    68                 X=npy.append(X,x_segment)
    69                 Y=npy.append(Y,y_segment)
    70                 S=npy.append(S,s_segment)
     68                X=np.append(X,x_segment)
     69                Y=np.append(Y,y_segment)
     70                S=np.append(S,s_segment)
    7171
    72         X=npy.append(X,x[nods-1])
    73         Y=npy.append(Y,y[nods-1])
     72        X=np.append(X,x[nods-1])
     73        Y=np.append(Y,y[nods-1])
    7474       
    7575        #Number of nodes:
    … …  
    7777       
    7878        #Compute Z
    79         Z=npy.zeros(numberofnodes)
     79        Z=np.zeros(numberofnodes)
    8080       
    8181        #New mesh and Data interpolation
    … …  
    8888       
    8989                #Compute index
    90                 index=npy.array([list(range(1,numberofnodes)),list(range(2,numberofnodes+1))]).T
     90                index=np.array([list(range(1,numberofnodes)),list(range(2,numberofnodes+1))]).T
    9191       
    9292        else:
    … …  
    102102       
    103103                #Some useful parameters
    104                 layers=int(npy.ceil(npy.mean(md.geometry.thickness)/res_v))
     104                layers=int(np.ceil(np.mean(md.geometry.thickness)/res_v))
    105105                nodesperlayer=int(numberofnodes)
    106106                nodestot=int(nodesperlayer*layers)
    … …  
    109109       
    110110                #initialization
    111                 X3=npy.zeros(nodesperlayer*layers)
    112                 Y3=npy.zeros(nodesperlayer*layers)
    113                 Z3=npy.zeros(nodesperlayer*layers)
    114                 S3=npy.zeros(nodesperlayer*layers)
    115                 index3=npy.zeros((elementstot,4))
     111                X3=np.zeros(nodesperlayer*layers)
     112                Y3=np.zeros(nodesperlayer*layers)
     113                Z3=np.zeros(nodesperlayer*layers)
     114                S3=np.zeros(nodesperlayer*layers)
     115                index3=np.zeros((elementstot,4))
    116116       
    117117                #Get new coordinates in 3d
    … …  
    123123       
    124124                        if i<layers-1:  #Build index3 with quads
    125                                 ids=npy.vstack((npy.arange(i,nodestot-layers,layers),npy.arange(i+1,nodestot-layers,layers),npy.arange(i+layers+1,nodestot,layers),npy.arange(i+layers,nodestot,layers))).T
     125                                ids=np.vstack((np.arange(i,nodestot-layers,layers),np.arange(i+1,nodestot-layers,layers),np.arange(i+layers+1,nodestot,layers),np.arange(i+layers,nodestot,layers))).T
    126126                                index3[(i-1)*elementsperlayer:i*elementsperlayer,:]=ids
    127127
    128128                #Interpolation of data on specified points
    129                 data_interp=InterpFromMeshToMesh3d(md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.z,data,X3,Y3,Z3,npy.nan)
     129                data_interp=InterpFromMeshToMesh3d(md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.z,data,X3,Y3,Z3,np.nan)
    130130       
    131131                #build outputs

  • issm/trunk-jpl/src/py3/interp/averaging.py

    r19895 r21255  
    1 import numpy as npy
     1import numpy as np
    22from GetAreas import GetAreas
    33from scipy.sparse import csc_matrix
    … …  
    3636        #initialization
    3737        if layer==0:
    38                 weights=npy.zeros(md.mesh.numberofvertices,)
     38                weights=np.zeros(md.mesh.numberofvertices,)
    3939                data=data.flatten(1)
    4040        else:
    41                 weights=npy.zeros(md.mesh.numberofvertices2d,)
     41                weights=np.zeros(md.mesh.numberofvertices2d,)
    4242                data=data[(layer-1)*md.mesh.numberofvertices2d+1:layer*md.mesh.numberofvertices2d,:]
    4343       
    … …  
    6666        index=index-1 # since python indexes starting from zero
    6767        line=index.flatten(1)
    68         areas=npy.vstack(areas).reshape(-1,)
    69         summation=1./rep*npy.ones(rep,)
     68        areas=np.vstack(areas).reshape(-1,)
     69        summation=1./rep*np.ones(rep,)
    7070        linesize=rep*numberofelements
    7171       
    7272        #update weights that holds the volume of all the element holding the node i
    73         weights=csc_matrix( (npy.tile(areas,(rep,1)).reshape(-1,),(line,npy.zeros(linesize,))), shape=(numberofnodes,1))
     73        weights=csc_matrix( (np.tile(areas,(rep,1)).reshape(-1,),(line,np.zeros(linesize,))), shape=(numberofnodes,1))
    7474       
    7575        #initialization
    7676        if len(data)==numberofelements:
    77                 average_node=csc_matrix( (npy.tile(areas*data,(rep,1)).reshape(-1,),(line,npy.zeros(linesize,))), shape=(numberofnodes,1))
     77                average_node=csc_matrix( (np.tile(areas*data,(rep,1)).reshape(-1,),(line,np.zeros(linesize,))), shape=(numberofnodes,1))
    7878                average_node=average_node/weights
    7979                average_node = csc_matrix(average_node)
    … …  
    8282
    8383        #loop over iteration
    84         for i in npy.arange(1,iterations+1):
    85                 average_el=npy.asarray(npy.dot(average_node.todense()[index].reshape(numberofelements,rep),npy.vstack(summation))).reshape(-1,)
    86                 average_node=csc_matrix( (npy.tile(areas*average_el.reshape(-1),(rep,1)).reshape(-1,),(line,npy.zeros(linesize,))), shape=(numberofnodes,1))
     84        for i in np.arange(1,iterations+1):
     85                average_el=np.asarray(np.dot(average_node.todense()[index].reshape(numberofelements,rep),np.vstack(summation))).reshape(-1,)
     86                average_node=csc_matrix( (np.tile(areas*average_el.reshape(-1),(rep,1)).reshape(-1,),(line,np.zeros(linesize,))), shape=(numberofnodes,1))
    8787                average_node=average_node/weights
    8888                average_node=csc_matrix(average_node)
    8989       
    9090        #return output as a full matrix (C code do not like sparse matrices)
    91         average=npy.asarray(average_node.todense()).reshape(-1,)
     91        average=np.asarray(average_node.todense()).reshape(-1,)
    9292
    9393        return average

  • issm/trunk-jpl/src/py3/interp/holefiller.py

    r19895 r21255  
    1 import numpy as npy
     1import numpy as np
    22from scipy.spatial import cKDTree
    33
    … …  
    3030        filled=data
    3131       
    32         XYGood=npy.dstack([x[goodids],y[goodids]])[0]
    33         XYBad=npy.dstack([x[badids],y[badids]])[0]
     32        XYGood=np.dstack([x[goodids],y[goodids]])[0]
     33        XYBad=np.dstack([x[badids],y[badids]])[0]
    3434        tree=cKDTree(XYGood)
    3535        nearest=tree.query(XYBad,k=knn)[1]
    … …  
    4242                        for j in range(knn):
    4343                                neardat.append(filled[goodids][nearest[i][j]])
    44                                 filled[badids[i]]=npy.mean(neardat)
     44                                filled[badids[i]]=np.mean(neardat)
    4545                               
    4646        return filled

  • issm/trunk-jpl/src/py3/interp/interp.py

    r19895 r21255  
    11# module for inperpolating/smoothing data
    2 import numpy as npy
     2import numpy as np
    33from scipy.interpolate import CloughTocher2DInterpolator, Rbf
    44from scipy.spatial import cKDTree
    … …  
    5151        xlim=[min(xi)-dx,max(xi)+dx]
    5252        ylim=[min(yi)-dy,max(yi)+dy]
    53         xflag=npy.logical_and(x>xlim[0],x<xlim[1])
    54         yflag=npy.logical_and(y>ylim[0],y<ylim[1])
    55         bothind=npy.nonzero(npy.logical_and(xflag,yflag))
     53        xflag=np.logical_and(x>xlim[0],x<xlim[1])
     54        yflag=np.logical_and(y>ylim[0],y<ylim[1])
     55        bothind=np.nonzero(np.logical_and(xflag,yflag))
    5656        subdata=data[bothind]
    5757        subx=x[bothind]
    5858        suby=y[bothind]
    59         points=npy.array([subx,suby]).T
     59        points=np.array([subx,suby]).T
    6060
    6161        # mask out any nan's in the data and corresponding coordinate points
    62         mask=npy.isnan(subdata)
    63         ind=npy.nonzero(mask)[0]
     62        mask=np.isnan(subdata)
     63        ind=np.nonzero(mask)[0]
    6464        if len(ind) and fill_nans:
    6565                print("         WARNING: filling nans using spline fit through good data points, which may or may not be appropriate. Check results carefully.")
    66         subdata=npy.delete(subdata,ind)
    67         points=npy.delete(points,ind,axis=0)
     66        subdata=np.delete(subdata,ind)
     67        points=np.delete(points,ind,axis=0)
    6868
    6969        if maxiter:
    … …  
    7676        if not fill_nans:
    7777                # identify nan's in xi,yi using nearest neighbors
    78                 xyinterp=npy.dstack([xi,yi])[0]
    79                 xg,yg=npy.meshgrid(subx,suby)
    80                 xydata=npy.dstack([subx,suby])[0]
     78                xyinterp=np.dstack([xi,yi])[0]
     79                xg,yg=np.meshgrid(subx,suby)
     80                xydata=np.dstack([subx,suby])[0]
    8181                tree=cKDTree(xydata)
    8282                nearest=tree.query(xyinterp)[1]
    83                 pos=npy.nonzero(npy.isnan(subdata[nearest]))
     83                pos=np.nonzero(np.isnan(subdata[nearest]))
    8484                interpdata[pos]=subdata[nearest][pos]
    8585
    8686        return interpdata
    8787#}}}
    88 def GridSplineToMesh2d(x,y,data,xi,yi,default_value=npy.nan,plotonly=False,fill_nans=False):#{{{
     88def GridSplineToMesh2d(x,y,data,xi,yi,default_value=np.nan,plotonly=False,fill_nans=False):#{{{
    8989        '''
    9090        python analog to InterpFromGridToMesh.  This routine uses
    … …  
    108108
    109109        Usage:
    110                 interpdata=GridToMesh(x,y,data,xi,yi,default_value=npy.nan,plotonly=False,fill_nans=False)
     110                interpdata=GridToMesh(x,y,data,xi,yi,default_value=np.nan,plotonly=False,fill_nans=False)
    111111
    112112        Examples:
    … …  
    114114        '''
    115115
    116         if npy.ndim(x)==2:
     116        if np.ndim(x)==2:
    117117                x=x.reshape(-1,)
    118         if npy.ndim(y)==2:
     118        if np.ndim(y)==2:
    119119                y=y.reshape(-1,)
    120120        if len(x) != data.shape[1]+1 and len(x) != data.shape[1]:
    … …  
    134134        if len(x)==data.shape[1] and len(y)==data.shape[0]:
    135135                print('         x,y taken to define the center of data grid cells')
    136                 xind=npy.nonzero(npy.logical_and(x>xlim[0],x<xlim[1]))[0]
    137                 yind=npy.nonzero(npy.logical_and(y>ylim[0],y<ylim[1]))[0]
    138                 xg,yg=npy.meshgrid(x[xind],y[yind])
     136                xind=np.nonzero(np.logical_and(x>xlim[0],x<xlim[1]))[0]
     137                yind=np.nonzero(np.logical_and(y>ylim[0],y<ylim[1]))[0]
     138                xg,yg=np.meshgrid(x[xind],y[yind])
    139139                subdata=data[yind[0]:yind[-1]+1,xind[0]:xind[-1]+1]
    140140        elif len(x)==data.shape[1]+1 and len(y)==data.shape[0]+1:
    141141                print('         x,y taken to define the corners of data grid cells')
    142                 xcenter=npy.fromiter(((x[i]+x[i+1])/2 for i in range(len(x)-1)),npy.float)
    143                 ycenter=npy.fromiter(((y[i]+y[i+1])/2 for i in range(len(y)-1)),npy.float)
    144                 xind=npy.nonzero(npy.logical_and(xcenter>xlim[0],xcenter<xlim[1]))[0]
    145                 yind=npy.nonzero(npy.logical_and(ycenter>ylim[0],ycenter<ylim[1]))[0]
    146                 xg,yg=npy.meshgrid(xcenter[xind],ycenter[yind])
     142                xcenter=np.fromiter(((x[i]+x[i+1])/2 for i in range(len(x)-1)),np.float)
     143                ycenter=np.fromiter(((y[i]+y[i+1])/2 for i in range(len(y)-1)),np.float)
     144                xind=np.nonzero(np.logical_and(xcenter>xlim[0],xcenter<xlim[1]))[0]
     145                yind=np.nonzero(np.logical_and(ycenter>ylim[0],ycenter<ylim[1]))[0]
     146                xg,yg=np.meshgrid(xcenter[xind],ycenter[yind])
    147147                subdata=data[yind[0]:yind[-1]+1,xind[0]:xind[-1]+1]
    148148        else:
    149149                raise ValueError('x and y have inconsistent sizes: both should have length ncols(data)/nrows(data) or ncols(data)+1/nrows(data)+1')
    150150
    151         points=npy.array([xg.ravel(),yg.ravel()]).T
     151        points=np.array([xg.ravel(),yg.ravel()]).T
    152152        flatsubdata=subdata.ravel()
    153153
    154154        if plotonly:
    155                 plt.imshow(npy.flipud(subdata),origin='upper')
     155                plt.imshow(np.flipud(subdata),origin='upper')
    156156                plt.show()
    157157                return
    158158
    159159        # mask out any nan's in the data and corresponding coordinate points
    160         mask=npy.isnan(flatsubdata)
    161         ind=npy.nonzero(mask)[0]
     160        mask=np.isnan(flatsubdata)
     161        ind=np.nonzero(mask)[0]
    162162        if len(ind) and fill_nans:
    163163                print("         WARNING: filling nans using spline fit through good data points, which may or may not be appropriate. Check results carefully.")
    164         goodsubdata=npy.delete(flatsubdata,ind)
    165         goodpoints=npy.delete(points,ind,axis=0)
     164        goodsubdata=np.delete(flatsubdata,ind)
     165        goodpoints=np.delete(points,ind,axis=0)
    166166
    167167        # create spline and index spline at mesh points
    … …  
    171171        if not fill_nans:
    172172                # identify nan's in xi,yi using nearest neighbors
    173                 xyinterp=npy.dstack([xi,yi])[0]
    174                 xydata=npy.dstack([xg.ravel(),yg.ravel()])[0]
     173                xyinterp=np.dstack([xi,yi])[0]
     174                xydata=np.dstack([xg.ravel(),yg.ravel()])[0]
    175175                tree=cKDTree(xydata)
    176176                nearest=tree.query(xyinterp)[1]
    177                 pos=npy.nonzero(npy.isnan(flatsubdata[nearest]))
     177                pos=np.nonzero(np.isnan(flatsubdata[nearest]))
    178178                interpdata[pos]=flatsubdata[nearest][pos]
    179179

  • issm/trunk-jpl/src/py3/inversions/parametercontroldrag.py

    r19895 r21255  
    2626
    2727        #weights
    28         weights=options.getfieldvalue('weights',npy.ones(md.mesh.numberofvertices))
    29         if npy.size(weights)!=md.mesh.numberofvertices:
     28        weights=options.getfieldvalue('weights',np.ones(md.mesh.numberofvertices))
     29        if np.size(weights)!=md.mesh.numberofvertices:
    3030                md.inversion.cost_functions_coefficients=ones(md.mesh.numberofvertices)
    3131        else:
    … …  
    3434        #nsteps
    3535        nsteps=options.getfieldvalue('nsteps',100);
    36         if (npy.size(nsteps)!=1) | (nsteps<=0) | (floor(nsteps)!=nsteps):
     36        if (np.size(nsteps)!=1) | (nsteps<=0) | (floor(nsteps)!=nsteps):
    3737                md.inversion.nsteps=100
    3838        else:
    … …  
    4141        #cm_min
    4242        cm_min=options.getfieldvalue('cm_min',ones(md.mesh.numberofvertices))
    43         if (npy.size(cm_min)==1):
     43        if (np.size(cm_min)==1):
    4444                md.inversion.min_parameters=cm_min*ones(md.mesh.numberofvertices)
    45         elif (npy.size(cm_min)==md.mesh.numberofvertices):
     45        elif (np.size(cm_min)==md.mesh.numberofvertices):
    4646                md.inversion.min_parameters=cm_min
    4747        else:
    … …  
    5050        #cm_max
    5151        cm_max=options.getfieldvalue('cm_max',250*ones(md.mesh.numberofvertices))
    52         if (npy.size(cm_max)==1):
     52        if (np.size(cm_max)==1):
    5353                md.inversion.max_parameters=cm_max*ones(md.mesh.numberofvertices)
    54         elif (npy.size(cm_max)==md.mesh.numberofvertices):
     54        elif (np.size(cm_max)==md.mesh.numberofvertices):
    5555                md.inversion.max_parameters=cm_max
    5656        else:
    … …  
    5959        #eps_cm
    6060        eps_cm=optoins.getfieldvalue('eps_cm',float('nan'))
    61         if (npy.size(eps_cm)~=1 | eps_cm<0 ):
     61        if (np.size(eps_cm)~=1 | eps_cm<0 ):
    6262                md.inversion.cost_function_threshold=float('nan')
    6363        else:
    … …  
    6666        #maxiter
    6767        maxiter=options.getfieldvalue('maxiter',10*ones(md.inversion.nsteps))
    68         if (npy.any(maxiter<0) | npy.any(floor(maxiter)~=maxiter)):
     68        if (np.any(maxiter<0) | np.any(floor(maxiter)~=maxiter)):
    6969                md.inversion.maxiter_per_step=10*ones(md.inversion.nsteps)
    7070        else:
    … …  
    7575        #cm_jump
    7676        cm_jump=options.getfieldvalue('cm_jump',0.8*ones(md.inversion.nsteps))
    77         if !npy.isreal(cm_jump):
     77        if !np.isreal(cm_jump):
    7878                md.inversion.step_threshold=0.8*ones(md.inversion.nsteps)
    7979        else:

  • issm/trunk-jpl/src/py3/materials/DepthAvgTempCond.py

    r19895 r21255  
    1 import numpy as npy
     1import numpy as np
    22from TMeltingPoint  import TMeltingPoint
    33
    … …  
    1616   alpha=G*H/k
    1717
    18    Tbar=npy.zeros(md.mesh.numberofvertices,)
     18   Tbar=np.zeros(md.mesh.numberofvertices,)
    1919
    2020   #find temperature average when we are below melting point:
    21    pos=npy.nonzero( Ts+alpha < Tpmp)
     21   pos=np.nonzero( Ts+alpha < Tpmp)
    2222   if pos:
    2323           Tbar[pos]=Ts[pos]+alpha[pos]/2
    2424
    25    pos=npy.nonzero( Ts+alpha>= Tpmp)
     25   pos=np.nonzero( Ts+alpha>= Tpmp)
    2626   if pos:
    2727           Tbar[pos]=Tpmp+(Tpmp**2-Ts[pos]**2)/2/alpha[pos]+ Tpmp*(Ts[pos]-Tpmp)/alpha[pos]
    2828   
    2929   #on ice shelf, easier:
    30    pos=npy.nonzero(md.mask.groundedice_levelset[0]<=0)
     30   pos=np.nonzero(md.mask.groundedice_levelset[0]<=0)
    3131   if pos:
    3232           Tbar[pos]=(Ts[pos]+Tpmp)/2

  • issm/trunk-jpl/src/py3/materials/TMeltingPoint.py

    r19895 r21255  
    1 import numpy as npy
     1import numpy as np
    22
    33def TMeltingPoint(reftemp,pressure):
    … …  
    1717
    1818        #ensure ref is same dimension as pressure
    19         ref=reftemp*npy.ones_like(pressure)
     19        ref=reftemp*np.ones_like(pressure)
    2020
    2121        return reftemp-beta*pressure

  • issm/trunk-jpl/src/py3/mech/analyticaldamage.py

    r19895 r21255  
    1 import numpy as npy
     1import numpy as np
    22from averaging import averaging
    33#from plotmodel import plotmodel
    … …  
    5454
    5555        if isinstance(sigmab,(int,float)):
    56                 sigmab=sigmab*npy.ones((md.mesh.numberofvertices,))
     56                sigmab=sigmab*np.ones((md.mesh.numberofvertices,))
    5757
    5858        # check inputs
    … …  
    6161        if not '2d' in md.mesh.__doc__:
    6262                raise Exception('only 2d (planview) model supported currently')
    63         if npy.any(md.flowequation.element_equation!=2):
     63        if np.any(md.flowequation.element_equation!=2):
    6464                print('Warning: the model has some non SSA elements. These will be treated like SSA elements')
    6565
    … …  
    7676        n=averaging(md,md.materials.rheology_n,0)
    7777       
    78         D=1.-(1.+a+a**2+b**2)**((n-1.)/(2.*n))/npy.abs(ex)**(1./n)*(T-sigmab)/B/(2.+a)/npy.sign(ex)
     78        D=1.-(1.+a+a**2+b**2)**((n-1.)/(2.*n))/np.abs(ex)**(1./n)*(T-sigmab)/B/(2.+a)/np.sign(ex)
    7979       
    8080        # D>1 where (2+a).*sign(ex)<0, compressive regions where high backstress needed
    81         pos=npy.nonzero(D>1)
     81        pos=np.nonzero(D>1)
    8282        D[pos]=0
    8383       
    84         backstress=npy.zeros((md.mesh.numberofvertices,))
     84        backstress=np.zeros((md.mesh.numberofvertices,))
    8585
    8686        # backstress to bring D down to one
    87         backstress[pos]=T[pos]-(1.-D[pos])*B[pos]*npy.sign(ex[pos])*(2.+a[pos])*npy.abs(ex[pos])**(1./n[pos])/(1.+a[pos]+a[pos]**2)**((n[pos]-1.)/2./n[pos])
     87        backstress[pos]=T[pos]-(1.-D[pos])*B[pos]*np.sign(ex[pos])*(2.+a[pos])*np.abs(ex[pos])**(1./n[pos])/(1.+a[pos]+a[pos]**2)**((n[pos]-1.)/2./n[pos])
    8888       
    89         pos=npy.nonzero(D<0)
     89        pos=np.nonzero(D<0)
    9090        #mask=ismember(1:md.mesh.numberofvertices,pos);
    9191        D[pos]=0
    9292       
    9393        # backstress to bring negative damage to zero
    94         backstress[pos]=T[pos]-(1.-D[pos])*B[pos]*npy.sign(ex[pos])*(2.+a[pos])*npy.abs(ex[pos])**(1./n[pos])/(1.+a[pos]+a[pos]**2)**((n[pos]-1.)/2./n[pos])
     94        backstress[pos]=T[pos]-(1.-D[pos])*B[pos]*np.sign(ex[pos])*(2.+a[pos])*np.abs(ex[pos])**(1./n[pos])/(1.+a[pos]+a[pos]**2)**((n[pos]-1.)/2./n[pos])
    9595       
    96         pos=npy.nonzero(backstress<0)
     96        pos=np.nonzero(backstress<0)
    9797        backstress[pos]=0
    9898       
    9999        # rigidity from Thomas relation for D=0 and backstress=0
    100         B=npy.sign(ex)/(2.+a)*(1.+a+a**2)**((n-1.)/2./n)*T/(npy.abs(ex)**(1./n))
    101         pos=npy.nonzero(B<0)
     100        B=np.sign(ex)/(2.+a)*(1.+a+a**2)**((n-1.)/2./n)*T/(np.abs(ex)**(1./n))
     101        pos=np.nonzero(B<0)
    102102        B[pos]=md.materials.rheology_B[pos]
    103103       

  • issm/trunk-jpl/src/py3/mech/backstressfrominversion.py

    r19895 r21255  
    1 import numpy as npy
     1import numpy as np
    22from averaging import averaging
    33from thomasparams import thomasparams
    … …  
    6565        if tempmask:
    6666                Bi=md.results.StressbalanceSolution.MaterialsRheologyBbar
    67                 pos=npy.nonzero(Bi>md.materials.rheology_B)
     67                pos=np.nonzero(Bi>md.materials.rheology_B)
    6868                Bi[pos]=md.materials.rheology_B[pos]
    6969       
    7070        # analytical backstress solution
    71         backstress=T-Bi*npy.sign(ex0)*(2+a0)*npy.abs(ex0)**(1./n)/((1+a0+a0**2+b0**2)**((n-1.)/2./n))
    72         backstress[npy.nonzero(backstress<0)]=0
     71        backstress=T-Bi*np.sign(ex0)*(2+a0)*np.abs(ex0)**(1./n)/((1+a0+a0**2+b0**2)**((n-1.)/2./n))
     72        backstress[np.nonzero(backstress<0)]=0
    7373
    7474        return backstress

  • issm/trunk-jpl/src/py3/mech/calcbackstress.py

    r19895 r21255  
    1 import numpy as npy
     1import numpy as np
    22from averaging import averaging
    33from thomasparams import thomasparams
    … …  
    6161       
    6262        # analytical backstress solution
    63         backstress=T-(1.-D)*B*npy.sign(ex0)*(2+a0)*npy.abs(ex0)**(1./n)/((1+a0+a0**2+b0**2)**((n-1.)/2./n))
    64         backstress[npy.nonzero(backstress<0)]=0
     63        backstress=T-(1.-D)*B*np.sign(ex0)*(2+a0)*np.abs(ex0)**(1./n)/((1+a0+a0**2+b0**2)**((n-1.)/2./n))
     64        backstress[np.nonzero(backstress<0)]=0
    6565
    6666        return backstress

  • issm/trunk-jpl/src/py3/mech/damagefrominversion.py

    r19895 r21255  
    1 import numpy as npy
     1import numpy as np
    22
    33def damagefrominversion(md):
    … …  
    2424        if any(md.flowequation.element_equation!=2):
    2525                raise Exception('Warning: the model has some non-SSA elements.  These will be treated like SSA elements')
    26         if npy.ndim(md.results.StressbalanceSolution.MaterialsRheologyBbar)==2:
     26        if np.ndim(md.results.StressbalanceSolution.MaterialsRheologyBbar)==2:
    2727                Bi=md.results.StressbalanceSolution.MaterialsRheologyBbar.reshape(-1,)
    2828        else:
    2929                Bi=md.results.StressbalanceSolution.MaterialsRheologyBbar
    30         if npy.ndim(md.materials.rheology_B)==2:
     30        if np.ndim(md.materials.rheology_B)==2:
    3131                BT=md.materials.rheology_B.reshape(-1,)
    3232        else:
    3333                BT=md.materials.rheology_B
    3434
    35         damage=npy.zeros_like(Bi)
     35        damage=np.zeros_like(Bi)
    3636
    3737        # Damage where Bi softer than B(T)
    38         pos=npy.nonzero(Bi<BT)[0]
     38        pos=np.nonzero(Bi<BT)[0]
    3939        damage[pos]=1.-Bi[pos]/BT[pos]
    4040       
    41         pos=npy.nonzero(damage<0)
     41        pos=np.nonzero(damage<0)
    4242        damage[pos]=0
    4343

  • issm/trunk-jpl/src/py3/mech/mechanicalproperties.py

    r19895 r21255  
    1 import numpy as npy
     1import numpy as np
    22from GetNodalFunctionsCoeff import GetNodalFunctionsCoeff
    33from results import results
    … …  
    2727        #       raise StandardError('only 2D model supported currently')
    2828
    29         if npy.any(md.flowequation.element_equation!=2):
     29        if np.any(md.flowequation.element_equation!=2):
    3030                print('Warning: the model has some non SSA elements. These will be treated like SSA elements')
    3131
    … …  
    3535            if len(damage)!=md.mesh.numberofvertices:
    3636                raise ValueError('if damage is supplied it should be of size ' + md.mesh.numberofvertices)
    37             if npy.ndim(damage)==2:
     37            if np.ndim(damage)==2:
    3838                damage=damage.reshape(-1,)
    3939        else: damage=None
    4040
    41         if npy.ndim(vx)==2:
     41        if np.ndim(vx)==2:
    4242                vx=vx.reshape(-1,)
    43         if npy.ndim(vy)==2:
     43        if np.ndim(vy)==2:
    4444                vy=vy.reshape(-1,)
    4545       
    … …  
    4848        numberofvertices=md.mesh.numberofvertices
    4949        index=md.mesh.elements
    50         summation=npy.array([[1],[1],[1]])
    51         directionsstress=npy.zeros((numberofelements,4))
    52         directionsstrain=npy.zeros((numberofelements,4))
    53         valuesstress=npy.zeros((numberofelements,2))
    54         valuesstrain=npy.zeros((numberofelements,2))
     50        summation=np.array([[1],[1],[1]])
     51        directionsstress=np.zeros((numberofelements,4))
     52        directionsstrain=np.zeros((numberofelements,4))
     53        valuesstress=np.zeros((numberofelements,2))
     54        valuesstrain=np.zeros((numberofelements,2))
    5555       
    5656        #compute nodal functions coefficients N(x,y)=alpha x + beta y +gamma
    … …  
    6060        vxlist=vx[index-1]/md.constants.yts
    6161        vylist=vy[index-1]/md.constants.yts
    62         ux=npy.dot((vxlist*alpha),summation).reshape(-1,)
    63         uy=npy.dot((vxlist*beta),summation).reshape(-1,)
    64         vx=npy.dot((vylist*alpha),summation).reshape(-1,)
    65         vy=npy.dot((vylist*beta),summation).reshape(-1,)
     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,)
    6666        uyvx=(vx+uy)/2.
    6767        #clear vxlist vylist
    6868       
    6969        #compute viscosity
    70         nu=npy.zeros((numberofelements,))
    71         B_bar=npy.dot(md.materials.rheology_B[index-1],summation/3.).reshape(-1,)
     70        nu=np.zeros((numberofelements,))
     71        B_bar=np.dot(md.materials.rheology_B[index-1],summation/3.).reshape(-1,)
    7272        power=((md.materials.rheology_n-1.)/(2.*md.materials.rheology_n)).reshape(-1,)
    7373        second_inv=(ux**2.+vy**2.+((uy+vx)**2.)/4.+ux*vy).reshape(-1,)
    7474       
    7575        #some corrections
    76         location=npy.nonzero(npy.logical_and(second_inv==0,power!=0))
     76        location=np.nonzero(np.logical_and(second_inv==0,power!=0))
    7777        nu[location]=10^18      #arbitrary maximum viscosity to apply where there is no effective shear
    7878       
    7979        if 'matice' in md.materials.__module__:
    80                 location=npy.nonzero(second_inv)
     80                location=np.nonzero(second_inv)
    8181                nu[location]=B_bar[location]/(second_inv[location]**power[location])
    82                 location=npy.nonzero(npy.logical_and(second_inv==0,power==0))
     82                location=np.nonzero(np.logical_and(second_inv==0,power==0))
    8383                nu[location]=B_bar[location]
    84                 location=npy.nonzero(npy.logical_and(second_inv==0,power!=0))
     84                location=np.nonzero(np.logical_and(second_inv==0,power!=0))
    8585                nu[location]=10^18
    8686        elif 'matdamageice' in md.materials.__module__ and damage is not None:
    8787                print('computing damage-dependent properties!')
    88                 Zinv=npy.dot(1-damage[index-1],summation/3.).reshape(-1,)
    89                 location=npy.nonzero(second_inv)
    90                 nu[location]=Zinv[location]*B_bar[location]/npy.power(second_inv[location],power[location])
    91                 location=npy.nonzero(npy.logical_and(second_inv==0,power==0))
     88                Zinv=np.dot(1-damage[index-1],summation/3.).reshape(-1,)
     89                location=np.nonzero(second_inv)
     90                nu[location]=Zinv[location]*B_bar[location]/np.power(second_inv[location],power[location])
     91                location=np.nonzero(np.logical_and(second_inv==0,power==0))
    9292                nu[location]=Zinv[location]*B_bar[location]
    9393                #clear Zinv
    … …  
    101101       
    102102        #compute principal properties of stress
    103         for i in npy.arange(numberofelements):
     103        for i in np.arange(numberofelements):
    104104       
    105105                #compute stress and strainrate matrices
    106                 stress=npy.array([ [tau_xx[i], tau_xy[i]], [tau_xy[i], tau_yy[i]] ])
    107                 strain=npy.array([ [ux[i], uyvx[i]], [uyvx[i], vy[i]] ])
     106                stress=np.array([ [tau_xx[i], tau_xy[i]], [tau_xy[i], tau_yy[i]] ])
     107                strain=np.array([ [ux[i], uyvx[i]], [uyvx[i], vy[i]] ])
    108108       
    109109                #eigenvalues and vectors for stress
    110                 value,directions=npy.linalg.eig(stress);
     110                value,directions=np.linalg.eig(stress);
    111111                idx=abs(value).argsort()[::-1] # sort in descending order
    112112                value=value[idx]
    … …  
    116116
    117117                #eigenvalues and vectors for strain
    118                 value,directions=npy.linalg.eig(strain);
     118                value,directions=np.linalg.eig(strain);
    119119                idx=abs(value).argsort()[::-1] # sort in descending order
    120120                value=value[idx]
    … …  
    133133        stress.principalvalue2=valuesstress[:,1]
    134134        stress.principalaxis2=directionsstress[:,2:4]
    135         stress.effectivevalue=1./npy.sqrt(2.)*npy.sqrt(stress.xx**2+stress.yy**2+2.*stress.xy**2)
     135        stress.effectivevalue=1./np.sqrt(2.)*np.sqrt(stress.xx**2+stress.yy**2+2.*stress.xy**2)
    136136        md.results.stress=stress
    137137       
    … …  
    144144        strainrate.principalvalue2=valuesstrain[:,1]*md.constants.yts
    145145        strainrate.principalaxis2=directionsstrain[:,2:4]
    146         strainrate.effectivevalue=1./npy.sqrt(2.)*npy.sqrt(strainrate.xx**2+strainrate.yy**2+2.*strainrate.xy**2)
     146        strainrate.effectivevalue=1./np.sqrt(2.)*np.sqrt(strainrate.xx**2+strainrate.yy**2+2.*strainrate.xy**2)
    147147        md.results.strainrate=strainrate
    148148       
    … …  
    155155        deviatoricstress.principalvalue2=valuesstress[:,1]
    156156        deviatoricstress.principalaxis2=directionsstress[:,2:4]
    157         deviatoricstress.effectivevalue=1./npy.sqrt(2.)*npy.sqrt(stress.xx**2+stress.yy**2+2.*stress.xy**2)
     157        deviatoricstress.effectivevalue=1./np.sqrt(2.)*np.sqrt(stress.xx**2+stress.yy**2+2.*stress.xy**2)
    158158        md.results.deviatoricstress=deviatoricstress
    159159

  • issm/trunk-jpl/src/py3/mech/robintemperature.py

    r19895 r21255  
    1 import numpy as npy
     1import numpy as np
    22from scipy.special import erf
    33
    … …  
    3434       
    3535        #create vertical coordinate variable
    36         zstar=npy.sqrt(2.*alphaT*thickness/accumrate)
     36        zstar=np.sqrt(2.*alphaT*thickness/accumrate)
    3737       
    38         tprofile=surftemp+npy.sqrt(2.*thickness*npy.pi/accumrate/alphaT)*(-heatflux)/2./rho/c*(erf(z/zstar)-erf(thickness/zstar))
     38        tprofile=surftemp+np.sqrt(2.*thickness*np.pi/accumrate/alphaT)*(-heatflux)/2./rho/c*(erf(z/zstar)-erf(thickness/zstar))
    3939
    4040        return tprofile

  • issm/trunk-jpl/src/py3/mech/steadystateiceshelftemp.py

    r19895 r21255  
    1 import numpy as npy
     1import numpy as np
    22
    33def steadystateiceshelftemp(md,surfacetemp,basaltemp):
    … …  
    4242        temperature=(Ts+Tb)/2  # where wi~=0
    4343       
    44         pos=npy.nonzero(abs(wi)>=1e-4) # to avoid division by zero
     44        pos=np.nonzero(abs(wi)>=1e-4) # to avoid division by zero
    4545
    46         npy.seterr(over='raise',divide='raise') # raise errors if floating point exceptions are encountered in following calculation
     46        np.seterr(over='raise',divide='raise') # raise errors if floating point exceptions are encountered in following calculation
    4747        #calculate depth-averaged temperature (in Celsius)
    4848        try:
    49                 temperature[pos]=-( (Tb[pos]-Ts[pos])*ki/wi[pos] + Hi[pos]*Tb[pos] - (Hi[pos]*Ts[pos] + (Tb[pos]-Ts[pos])*ki/wi[pos])*npy.exp(Hi[pos]*wi[pos]/ki) )/( Hi[pos]*(npy.exp(Hi[pos]*wi[pos]/ki)-1))
     49                temperature[pos]=-( (Tb[pos]-Ts[pos])*ki/wi[pos] + Hi[pos]*Tb[pos] - (Hi[pos]*Ts[pos] + (Tb[pos]-Ts[pos])*ki/wi[pos])*np.exp(Hi[pos]*wi[pos]/ki) )/( Hi[pos]*(np.exp(Hi[pos]*wi[pos]/ki)-1))
    5050        except FloatingPointError:
    5151                print('steadystateiceshelf warning: overflow encountered in multipy/divide/exp, trying another formulation.')
    52                 temperature[pos]=-( ((Tb[pos]-Ts[pos])*ki/wi[pos] + Hi[pos]*Tb[pos])/npy.exp(Hi[pos]*wi[pos]/ki) - Hi[pos]*Ts[pos] + (Tb[pos]-Ts[pos])*ki/wi[pos])/( Hi[pos]*(1-npy.exp(-Hi[pos]*wi[pos]/ki)))
     52                temperature[pos]=-( ((Tb[pos]-Ts[pos])*ki/wi[pos] + Hi[pos]*Tb[pos])/np.exp(Hi[pos]*wi[pos]/ki) - Hi[pos]*Ts[pos] + (Tb[pos]-Ts[pos])*ki/wi[pos])/( Hi[pos]*(1-np.exp(-Hi[pos]*wi[pos]/ki)))
    5353       
    5454        #temperature should not be less than surface temp
    55         pos=npy.nonzero(temperature<Ts)
     55        pos=np.nonzero(temperature<Ts)
    5656        temperature[pos]=Ts[pos]
    5757       
    5858        # NaN where melt rates are too high (infinity/infinity in exponential)
    59         pos=npy.nonzero(npy.isnan(temperature))
     59        pos=np.nonzero(np.isnan(temperature))
    6060        temperature[pos]=Ts[pos]
    6161       

  • issm/trunk-jpl/src/py3/mech/thomasparams.py

    r19895 r21255  
    1 import numpy as npy
     1import numpy as np
    22from averaging import averaging
    33
    … …  
    7575       
    7676        # checks: any of e1 or e2 equal to zero?
    77         pos=npy.nonzero(e1==0)
    78         if npy.any(pos==1):
     77        pos=np.nonzero(e1==0)
     78        if np.any(pos==1):
    7979                print('WARNING: first principal strain rate equal to zero.  Value set to 1e-13 s^-1')
    8080                e1[pos]=1.e-13
    81         pos=npy.nonzero(e2==0)
    82         if npy.any(pos==1):
     81        pos=np.nonzero(e2==0)
     82        if np.any(pos==1):
    8383                print('WARNING: second principal strain rate equal to zero.  Value set to 1e-13 s^-1')
    8484                e2[pos]=1.e-13
    … …  
    8989       
    9090        if coordsys=='principal':
    91                 b=npy.zeros((md.mesh.numberofvertices,))
     91                b=np.zeros((md.mesh.numberofvertices,))
    9292                ex=e1
    9393                a=e2/e1
    94                 pos=npy.nonzero(npy.logical_and(e1<0,e2>0)) # longitudinal compression and lateral tension
     94                pos=np.nonzero(np.logical_and(e1<0,e2>0)) # longitudinal compression and lateral tension
    9595                a[pos]=e1[pos]/e2[pos]
    9696                ex[pos]=e2[pos]
    97                 pos2=npy.nonzero(e1<0 & e2<0 & npy.abs(e1)<npy.abs(e2)) # lateral and longitudinal compression
     97                pos2=np.nonzero(e1<0 & e2<0 & np.abs(e1)<np.abs(e2)) # lateral and longitudinal compression
    9898                a[pos2]=e1[pos2]/e2[pos2]
    9999                ex[pos2]=e2[pos2]
    100                 pos3=npy.nonzero(e1>0 & e2>0 & npy.abs(e1)<npy.abs(e2)) # lateral and longitudinal tension
     100                pos3=np.nonzero(e1>0 & e2>0 & np.abs(e1)<np.abs(e2)) # lateral and longitudinal tension
    101101                a[pos3]=e1[pos3]/e2[pos3]
    102102                ex[pos3]=e2[pos3]
    103                 ind=npy.nonzero(e1<0 & e2<0)
     103                ind=np.nonzero(e1<0 & e2<0)
    104104                a[ind]=-a[ind] # where both strain rates are compressive, enforce negative alpha
    105                 sigxx=(npy.abs(ex)/((1.+a+a**2)**((n-1.)/2.)))**(1./n)*B
     105                sigxx=(np.abs(ex)/((1.+a+a**2)**((n-1.)/2.)))**(1./n)*B
    106106        elif coordsys=='xy':
    107107                ex=exx
    … …  
    110110        elif coordsys=='longitudinal':
    111111                # using longitudinal strain rates defined by observed velocity vector
    112                 velangle=npy.arctan(md.initialization.vy/md.initialization.vx)
    113                 pos=npy.nonzero(md.initialization.vx==0)
    114                 velangle[pos]=npy.pi/2
    115                 ex=0.5*(exx+eyy)+0.5*(exx-eyy)*npy.cos(2.*velangle)+exy*npy.sin(2.*velangle)
     112                velangle=np.arctan(md.initialization.vy/md.initialization.vx)
     113                pos=np.nonzero(md.initialization.vx==0)
     114                velangle[pos]=np.pi/2
     115                ex=0.5*(exx+eyy)+0.5*(exx-eyy)*np.cos(2.*velangle)+exy*np.sin(2.*velangle)
    116116                ey=exx+eyy-ex # trace of strain rate tensor is invariant
    117                 exy=-0.5*(exx-eyy)*npy.sin(2.*velangle)+exy*npy.cos(2.*velangle)
     117                exy=-0.5*(exx-eyy)*np.sin(2.*velangle)+exy*np.cos(2.*velangle)
    118118                a=ey/ex
    119119                b=exy/ex
    … …  
    124124        # a < -1 in areas of strong lateral compression or longitudinal compression and
    125125        # theta flips sign at a = -2
    126         pos=npy.nonzero(npy.abs((npy.abs(a)-2.))<1.e-3)
     126        pos=np.nonzero(np.abs((np.abs(a)-2.))<1.e-3)
    127127        if len(pos)>0:
    128128                print('Warning: ', len(pos), ' vertices have alpha within 1e-3 of -2')
    … …  
    131131        if eq=='Weertman1D':
    132132                theta=1./8
    133                 a=npy.zeros((md.mesh.numberofvertices,))
     133                a=np.zeros((md.mesh.numberofvertices,))
    134134        elif eq=='Weertman2D':
    135135                theta=1./9
    136                 a=npy.ones((md.mesh.numberofvertices,))
     136                a=np.ones((md.mesh.numberofvertices,))
    137137        elif eq=='Thomas':
    138                 theta=((1.+a+a**2+b**2)**((n-1.)/2.))/(npy.abs(2.+a)**n)
     138                theta=((1.+a+a**2+b**2)**((n-1.)/2.))/(np.abs(2.+a)**n)
    139139        else:
    140140                raise ValueError('argument passed to "eq" not valid')

  • issm/trunk-jpl/src/py3/plot/applyoptions.py

    r19895 r21255  
    1 import numpy as npy
     1import numpy as np
    22from cmaptools import truncate_colormap
    33from plot_contour import plot_contour
    … …  
    230230        cb.locator=MaxNLocator(nbins=5) # default 5 ticks
    231231        if options.exist('colorbartickspacing'):
    232                 locs=npy.arange(lims[0],lims[1]+1,options.getfieldvalue('colorbartickspacing'))
     232                locs=np.arange(lims[0],lims[1]+1,options.getfieldvalue('colorbartickspacing'))
    233233                cb.set_ticks(locs)
    234234                if options.exist('colorbarlines'):
    235                         locs=npy.arange(lims[0],lims[1]+1,options.getfieldvalue('colorbarlines'))
    236                         cb.add_lines(locs,['k' for i in range(len(locs))],npy.ones_like(locs))
     235                        locs=np.arange(lims[0],lims[1]+1,options.getfieldvalue('colorbarlines'))
     236                        cb.add_lines(locs,['k' for i in range(len(locs))],np.ones_like(locs))
    237237                        if options.exist('colorbarlineatvalue'):
    238238                                locs=options.getfieldvalue('colorbarlineatvalue')
    239239                                colors=options.getfieldvalue('colorbarlineatvaluecolor',['k' for i in range (len(locs))])
    240                                 widths=options.getfieldvalue('colorbarlineatvaluewidth',npy.ones_like(locs))
     240                                widths=options.getfieldvalue('colorbarlineatvaluewidth',np.ones_like(locs))
    241241                                cb.add_lines(locs,colors,widths)
    242242                                if options.exist('colorbartitle'):

  • issm/trunk-jpl/src/py3/plot/checkplotoptions.py

    r19895 r21255  
    1 import numpy as npy
     1import numpy as np
    22
    33def checkplotoptions(md,options):
    … …  
    7878                sizelist=[12]
    7979        if len(sizelist)==1:
    80                 sizelist=npy.tile(sizelist,numtext)
     80                sizelist=np.tile(sizelist,numtext)
    8181
    8282                # font color
    … …  
    8888                        colorlist=['k']
    8989                if len(colorlist)==1:
    90                         colorlist=npy.tile(colorlist,numtext)
     90                        colorlist=np.tile(colorlist,numtext)
    9191
    9292                # textweight
    … …  
    9898                        weightlist=['normal']
    9999                if len(weightlist)==1:
    100                         weightlist=npy.tile(weightlist,numtext)
     100                        weightlist=np.tile(weightlist,numtext)
    101101
    102102                # text rotation
    … …  
    108108                        rotationlist=[0]
    109109                if len(rotationlist)==1:
    110                         rotationlist=npy.tile(rotationlist,numtext)
     110                        rotationlist=np.tile(rotationlist,numtext)
    111111
    112112                options.changefieldvalue('text',textlist)
    … …  
    130130                if type(expdispvalues)==str:
    131131                        expdispvalues=[expdispvalues]
    132                 for i in npy.arange(len(expdispvalues)):
     132                for i in np.arange(len(expdispvalues)):
    133133                        expdispvaluesarray.append(expdispvalues[i])
    134134                        if len(expstylevalues)>i:

  • issm/trunk-jpl/src/py3/plot/colormaps/cmaptools.py

    r19895 r21255  
    1 import numpy as npy
     1import numpy as np
    22
    33try:
    … …  
    2121
    2222        new_cmap = mpl.colors.LinearSegmentedColormap.from_list('trunc({n},{a:.2f},{b:.2f})'.format(n=cmap.name,
    23                 a=minval, b=maxval), cmap(npy.linspace(minval, maxval, n)))
     23                a=minval, b=maxval), cmap(np.linspace(minval, maxval, n)))
    2424       
    2525        return new_cmap

  • issm/trunk-jpl/src/py3/plot/plot_contour.py

    r19895 r21255  
    2323                pass
    2424        elif datatype==3: # quiver (vector) data
    25                 data=npy.sqrt(datain**2)
     25                data=np.sqrt(datain**2)
    2626        else:
    2727                raise ValueError('datatype not supported in call to plot_contour')

  • issm/trunk-jpl/src/py3/plot/plot_overlay.py

    r19895 r21255  
    1 import numpy as npy
     1import numpy as np
    22from processmesh import processmesh
    33from processdata import processdata
    … …  
    2828        if data=='none' or data==None:
    2929                imageonly=1
    30                 data=npy.float('nan')*npy.ones((md.mesh.numberofvertices,))
     30                data=np.float('nan')*np.ones((md.mesh.numberofvertices,))
    3131                datatype=1
    3232        else:
    … …  
    7070
    7171        # normalize array
    72         arr=arr/npy.float(npy.max(arr.ravel()))
     72        arr=arr/np.float(np.max(arr.ravel()))
    7373        arr=1.-arr # somehow the values got flipped
    7474
    … …  
    9696        dy=trans[5]     
    9797       
    98         xarr=npy.arange(xmin,xmax,dx)
    99         yarr=npy.arange(ymin,ymax,-dy) # -dy since origin='upper' (not sure how robust this is)
    100         xg,yg=npy.meshgrid(xarr,yarr)
     98        xarr=np.arange(xmin,xmax,dx)
     99        yarr=np.arange(ymin,ymax,-dy) # -dy since origin='upper' (not sure how robust this is)
     100        xg,yg=np.meshgrid(xarr,yarr)
    101101        overlaylims=options.getfieldvalue('overlaylims',[min(arr.ravel()),max(arr.ravel())])
    102102        norm=mpl.colors.Normalize(vmin=overlaylims[0],vmax=overlaylims[1])
    … …  
    118118                #test
    119119                #m.ax=ax
    120                 meridians=npy.arange(-180.,181.,1.)
    121                 parallels=npy.arange(-80.,80.,1.)
     120                meridians=np.arange(-180.,181.,1.)
     121                parallels=np.arange(-80.,80.,1.)
    122122                m.drawparallels(parallels,labels=[0,0,1,1]) # labels=[left,right,top,bottom]
    123123                m.drawmeridians(meridians,labels=[1,1,0,0])
    124124                m.drawcoastlines()
    125                 pc=m.pcolormesh(xg, yg, npy.flipud(arr), cmap=mpl.cm.Greys, norm=norm, ax=ax)
     125                pc=m.pcolormesh(xg, yg, np.flipud(arr), cmap=mpl.cm.Greys, norm=norm, ax=ax)
    126126
    127127        else:
    128                 pc=ax.pcolormesh(xg, yg, npy.flipud(arr), cmap=mpl.cm.Greys, norm=norm)
     128                pc=ax.pcolormesh(xg, yg, np.flipud(arr), cmap=mpl.cm.Greys, norm=norm)
    129129       
    130130        #rasterization?

  • issm/trunk-jpl/src/py3/plot/plot_streamlines.py

    r19895 r21255  
    1 import numpy as npy
     1import numpy as np
    22import matplotlib.pyplot as plt
    33import matplotlib.tri as tri
    … …  
    4242
    4343    # format data for matplotlib streamplot function
    44     yg,xg=npy.mgrid[min(md.mesh.y):max(md.mesh.y):100j,min(md.mesh.x):max(md.mesh.x):100j]
     44    yg,xg=np.mgrid[min(md.mesh.y):max(md.mesh.y):100j,min(md.mesh.x):max(md.mesh.x):100j]
    4545    ug=griddata((x,y),u,(xg,yg),method='linear')
    4646    vg=griddata((x,y),v,(xg,yg),method='linear')
    … …  
    5757    if linewidth=='vel':
    5858        scale=options.getfieldvalue('streamlineswidthscale',3)
    59         vel=npy.sqrt(ug**2+vg**2)
    60         linewidth=scale*vel/npy.amax(vel)
     59        vel=np.sqrt(ug**2+vg**2)
     60        linewidth=scale*vel/np.amax(vel)
    6161        linewidth[linewidth<0.5]=0.5
    6262

  • issm/trunk-jpl/src/py3/plot/plot_unit.py

    r19895 r21255  
    44    import matplotlib as mpl
    55    import matplotlib.pyplot as plt
    6     import numpy as npy
     6    import numpy as np
    77except ImportError:
    88    print("could not import pylab, matplotlib has not been installed, no plotting capabilities enabled")
    … …  
    3939    #normalize colormap if clim/caxis specified
    4040    if options.exist('clim'):
    41         lims=options.getfieldvalue('clim',[npy.amin(data),npy.amax(data)])
     41        lims=options.getfieldvalue('clim',[np.amin(data),np.amax(data)])
    4242    elif options.exist('caxis'):
    43         lims=options.getfieldvalue('caxis',[npy.amin(data),npy.amax(data)])
     43        lims=options.getfieldvalue('caxis',[np.amin(data),np.amax(data)])
    4444    else:
    45         if npy.amin(data)==npy.amax(data):
    46             lims=[npy.amin(data)-0.5,npy.amax(data)+0.5]
     45        if np.amin(data)==np.amax(data):
     46            lims=[np.amin(data)-0.5,np.amax(data)+0.5]
    4747        else:
    48             lims=[npy.amin(data),npy.amax(data)]
     48            lims=[np.amin(data),np.amax(data)]
    4949    norm = mpl.colors.Normalize(vmin=lims[0], vmax=lims[1])
    5050    if datatype==1:
    … …  
    7070        vy=data[:,1]
    7171        #TODO write plot_quiver.py to handle this here
    72         color=npy.sqrt(vx**2+vy**2)
     72        color=np.sqrt(vx**2+vy**2)
    7373        scale=options.getfieldvalue('scale',1000)
    74         width=options.getfieldvalue('width',0.005*(npy.amax(x)-npy.amin(y)))
     74        width=options.getfieldvalue('width',0.005*(np.amax(x)-np.amin(y)))
    7575        headwidth=options.getfieldvalue('headwidth',3)
    7676        headlength=options.getfieldvalue('headlength',5)

  • issm/trunk-jpl/src/py3/plot/plotmodel.py

    r19895 r21255  
    1 import numpy as npy
     1import numpy as np
    22from plotoptions import plotoptions
    33
    … …  
    3535                nr=True
    3636        else:
    37                 nrows=npy.ceil(numberofplots/subplotwidth)
     37                nrows=np.ceil(numberofplots/subplotwidth)
    3838                nr=False
    3939       

  • issm/trunk-jpl/src/py3/plot/processdata.py

    r19895 r21255  
    11from math import isnan
    2 import numpy as npy
     2import numpy as np
    33
    44def processdata(md,data,options):
    … …  
    2626        numberofelements2d=md.mesh.numberofelements2d
    2727    else:
    28         numberofvertices2d=npy.nan
    29         numberofelements2d=npy.nan
     28        numberofvertices2d=np.nan
     29        numberofelements2d=np.nan
    3030   
    31     procdata=npy.copy(data)
     31    procdata=np.copy(data)
    3232   
    3333    #process patch
    … …  
    3737   
    3838    #get datasize
    39     if npy.ndim(procdata)==1:
    40         datasize=npy.array([len(procdata),1])
     39    if np.ndim(procdata)==1:
     40        datasize=np.array([len(procdata),1])
    4141    else:
    42         datasize=npy.shape(procdata)
     42        datasize=np.shape(procdata)
    4343        if len(datasize)>2:
    4444            raise ValueError('data passed to plotmodel has more than 2 dimensions; check that column vectors are rank-1')
    … …  
    4646    #process NaN's if any
    4747    nanfill=options.getfieldvalue('nan',-9999)
    48     if npy.any(npy.isnan(procdata)):
    49         lb=npy.min(data[~npy.isnan(data)])
    50         ub=npy.max(data[~npy.isnan(data)])
     48    if np.any(np.isnan(procdata)):
     49        lb=np.min(data[~np.isnan(data)])
     50        ub=np.max(data[~np.isnan(data)])
    5151        if lb==ub:
    5252            lb=lb-0.5
    5353            ub=ub+0.5
    5454            nanfill=lb-1
    55         procdata[npy.isnan(procdata)]=nanfill
     55        procdata[np.isnan(procdata)]=nanfill
    5656        options.addfielddefault('clim',[lb,ub])
    5757        options.addfielddefault('cmap_set_under','1')
    … …  
    9999   
    100100    #convert rank-2 array to rank-1
    101     if npy.ndim(procdata)==2 and npy.shape(procdata)[1]==1:
     101    if np.ndim(procdata)==2 and np.shape(procdata)[1]==1:
    102102        procdata=procdata.reshape(-1,)
    103103   

  • issm/trunk-jpl/src/py3/plot/processmesh.py

    r19895 r21255  
    11from math import isnan
    22import MatlabFuncs as m
    3 import numpy as npy
     3import numpy as np
    44
    55def processmesh(md,data,options):
    … …  
    3333                        z=md.mesh.z
    3434                else:
    35                         z=npy.zeros_like(md.mesh.x)
     35                        z=np.zeros_like(md.mesh.x)
    3636               
    3737                if 'elements2d' in dir(md.mesh):
Note: See TracChangeset for help on using the changeset viewer.