source: issm/trunk-jpl/src/py3/miscellaneous/prctile_issm.py@ 23677

import numpy as np
from scipy.interpolate import interp1d

#
#  wrapper for prctile to avoid using the matlab statistics toolbox.
#
def prctile_issm(x,p,dim):

        try:
                raise RuntimeError('hello world')

                #numpy has no interpolation method that matches matlab's percentile function
                #y = np.percentile(x,p,dim,interpolation='higher')
        except:
                if len(np.shape(x)) > 2:
                        raise RuntimeError('Number of dimensions #d not implemented.'+str(len(np.shape(x))))

                # presumably at least 1 input value has been given
                #       np.shape(integer) -> (), must be at least (1,)
                psize=np.shape(p) or (1,)
                if len(psize) > 1 and np.size(p,1)>1:
                        p=p.T
                
                xsize=np.shape(x) or (1,)
                if dim==2:
                        x=x.T

                #  check for any NaN in any columns
                if not np.isnan(x).any():
                        x=np.sort(x,axis=0)
                        n=np.size(x,0)

                        #  branch based on number of elements
                        if n>1:
                                #  set up percent values and interpolate
                                xi=[((i+0.5)*100/n) for i in range(n)]
                                # scipy's interp1d returns a function
                                y=interp1d(xi,x,axis=dim,bounds_error=False)
                                y=y(p)

                                #  fill in high and low values outside of interp range
                                if p>xi[n-1]:
                                        y=np.tile(x[n-1,:],1)
                                if p<xi[0]:
                                        y=np.tile(x[0,:],1)

                        #  if one value, just copy it
                        elif n==1:
                                y=np.tile(x[0,:],(len(p),1))

                        #  if no values, use NaN
                        else:
                                y=np.tile(float('NaN'),(size(p,0),size(x,0)))
                else:
                        #  must loop over columns, since number of elements could be different
                        y=np.zeros((np.size(p,0),np.size(x,1)))
                        for j in range(np.size(x,1)):
                        #  remove any NaN and recursively call column
                                y[:,j]=prctile_issm(x[np.where(not np.isnan(x[:,j]),j)],p)

                if (np.min(xsize)==1 and len(xsize) > 1 and xsize[dim]>1 and len(p) > 1 and psize[1]>1) or (np.min(xsize)> 1 and dim==2):
                        y=y.T

        return y