用python重寫ncl的wrf_user_intrp3d函數

不多說了,,直接把代碼貼上來了

如果對大家有幫助的話，發PAPER時請帶上我的名字HI.KONG哈哈哈，謝謝啦

__author__ = 'kong'
import numpy as np

##--use to interpolation--##

def interp3d(dat0,mat0,pos0,direction):
  sized=np.shape(dat0)
  sizem=np.shape(mat0)
  #y=0
  if direction<>0:##not top(big) bottom(small)
    dat0=dat0[::-1,...]
    mat0=mat0[::-1,...]
  if sized==sizem:
    #del(y)
    mysize=np.zeros(np.shape(sizem),int)
    mysize[0]=sizem[0]-1
    mysize[1:]=sizem[1:]
    meow=np.all([mat0[0:-1,...]>pos0,mat0[1:,...]<=pos0],0)
    mt=np.zeros(mysize)
    mb=np.zeros(mysize)
    dt=np.zeros(mysize)
    db=np.zeros(mysize)
    mt[meow]=mat0[0:-1,...][meow]
    mb[meow]=mat0[1:,...][meow]
    dt[meow]=dat0[0:-1,...][meow]
    db[meow]=dat0[1:,...][meow]
    mt=np.sum(mt,0)
    mb=np.sum(mb,0)
    dt=np.sum(dt,0)
    db=np.sum(db,0)
    dt=np.ma.array(dt,mask=(np.sum(meow,0)==0))
    y=(dt-db)/(mt-mb)*(pos0-mb)+db
  return(y)

__author__ = 'kong'
import numpy as np

def interpcs(data,i1,j1,i2,j2):
  meow=0
  if np.ndim(data)==2:
      meow=1
      data=np.ma.reshape(data,(1,np.ma.shape(data)[0],np.ma.shape(data)[1]))

  if np.abs(i2-i1)>=np.abs(j2-j1):
    if i2>=i1:
      i_index=np.arange(i1+1,i2,1.0)
    else:
      i_index=np.arange(i1-1,i2,-1.0)
    outdata=np.ma.zeros((np.shape(data)[0],np.abs(i2-i1)+1))
    outdata[:,0]=data[:,i1,j1]
    outdata[:,-1]=data[:,i2,j2]
    out_i=0
    for i in i_index:
      out_i=out_i+1
      j=(i-i1)*(j2-j1)/(i2-i1)+j1
      if np.floor(j)<>j:
        jjs=np.floor(j)
        jjb=np.ceil(j)
        ds=data[:,int(i),int(jjs)]
        db=data[:,int(i),int(jjb)]
        outdata[:,out_i]=(db-ds)/(jjb-jjs)*(j-jjs)+ds
      else:
        outdata[:,out_i]=data[:,int(i),int(j)]
  else:
    if j2>=j1:
      j_index=np.arange(j1+1,j2,1.0)
    else:
      j_index=np.arange(j1-1,j2,-1.0)
    outdata=np.ma.zeros((np.shape(data)[0],np.abs(j2-j1)+1))
    outdata[:,0]=data[:,i1,j1]
    outdata[:,-1]=data[:,i2,j2]
    out_j=0
    for j in j_index:
      out_j=out_j+1
      i=(j-j1)*(i2-i1)/(j2-j1)+i1
      if np.floor(i)<>i:
        iis=np.floor(i)
        iib=np.ceil(i)
        ds=data[:,int(iis),int(j)]
        db=data[:,int(iib),int(j)]
        outdata[:,out_j]=(db-ds)/(iib-iis)*(i-iis)+ds
      else:
        outdata[:,out_j]=data[:,int(i),int(j)]
  if meow==1:
      outdata=np.ma.reshape(outdata,(np.ma.shape(outdata)[1]))
  return(outdata)