How to use the mudpy.forward.lowpass function in mudpy

Eds[ktrace].data=hfilt(Eds[ktrace].data,BP[0],fsample,2)
                            Nds[ktrace].data=hfilt(Nds[ktrace].data,BP[0],fsample,2)
                            Zds[ktrace].data=hfilt(Zds[ktrace].data,BP[0],fsample,2)                        
                        else: #A low pass or bandpass filter has been requested
                            if kfault==0:
                                print('.... lowpass')
                            #Ess[ktrace].data=lfilt(Ess[ktrace].data,BP[0],fsample,2)
                            #Nss[ktrace].data=lfilt(Nss[ktrace].data,BP[0],fsample,2)
                            #Zss[ktrace].data=lfilt(Zss[ktrace].data,BP[0],fsample,2)
                            #Eds[ktrace].data=lfilt(Eds[ktrace].data,BP[0],fsample,2)
                            #Nds[ktrace].data=lfilt(Nds[ktrace].data,BP[0],fsample,2)
                            #Zds[ktrace].data=lfilt(Zds[ktrace].data,BP[0],fsample,2)
                            Ess[ktrace].data=lfilt(Ess[ktrace].data,BP,fsample,2)
                            Nss[ktrace].data=lfilt(Nss[ktrace].data,BP,fsample,2)
                            Zss[ktrace].data=lfilt(Zss[ktrace].data,BP,fsample,2)
                            Eds[ktrace].data=lfilt(Eds[ktrace].data,BP,fsample,2)
                            Nds[ktrace].data=lfilt(Nds[ktrace].data,BP,fsample,2)
                            Zds[ktrace].data=lfilt(Zds[ktrace].data,BP,fsample,2)
                        #bandpass=None
                    
                    ### HAAAAAACCCCKKKK!!!!
                    #if vord=='vel': #Apply low pass filter to data (** DIRTY HACK!**)
                    #    if kfault==0:
                    #        print(vord
                    #        print('Bandpassing hack on velocity...'
                    #    bandpass=[1./20,1./2]
                    #    fsample=1./Ess[0].stats.delta
                    #    Ess[ktrace].data=lfilt(Ess[ktrace].data,bandpass,fsample,2)
                    #    Nss[ktrace].data=lfilt(Nss[ktrace].data,bandpass,fsample,2)
                    #    Zss[ktrace].data=lfilt(Zss[ktrace].data,bandpass,fsample,2)
                    #    Eds[ktrace].data=lfilt(Eds[ktrace].data,bandpass,fsample,2)
                    #    Nds[ktrace].data=lfilt(Nds[ktrace].data,bandpass,fsample,2)

Zds[ktrace].data=hfilt(Zds[ktrace].data,BP[0],fsample,2)                        
                        else: #A low pass or bandpass filter has been requested
                            if kfault==0:
                                print('.... lowpass')
                            #Ess[ktrace].data=lfilt(Ess[ktrace].data,BP[0],fsample,2)
                            #Nss[ktrace].data=lfilt(Nss[ktrace].data,BP[0],fsample,2)
                            #Zss[ktrace].data=lfilt(Zss[ktrace].data,BP[0],fsample,2)
                            #Eds[ktrace].data=lfilt(Eds[ktrace].data,BP[0],fsample,2)
                            #Nds[ktrace].data=lfilt(Nds[ktrace].data,BP[0],fsample,2)
                            #Zds[ktrace].data=lfilt(Zds[ktrace].data,BP[0],fsample,2)
                            Ess[ktrace].data=lfilt(Ess[ktrace].data,BP,fsample,2)
                            Nss[ktrace].data=lfilt(Nss[ktrace].data,BP,fsample,2)
                            Zss[ktrace].data=lfilt(Zss[ktrace].data,BP,fsample,2)
                            Eds[ktrace].data=lfilt(Eds[ktrace].data,BP,fsample,2)
                            Nds[ktrace].data=lfilt(Nds[ktrace].data,BP,fsample,2)
                            Zds[ktrace].data=lfilt(Zds[ktrace].data,BP,fsample,2)
                        #bandpass=None
                    
                    ### HAAAAAACCCCKKKK!!!!
                    #if vord=='vel': #Apply low pass filter to data (** DIRTY HACK!**)
                    #    if kfault==0:
                    #        print(vord
                    #        print('Bandpassing hack on velocity...'
                    #    bandpass=[1./20,1./2]
                    #    fsample=1./Ess[0].stats.delta
                    #    Ess[ktrace].data=lfilt(Ess[ktrace].data,bandpass,fsample,2)
                    #    Nss[ktrace].data=lfilt(Nss[ktrace].data,bandpass,fsample,2)
                    #    Zss[ktrace].data=lfilt(Zss[ktrace].data,bandpass,fsample,2)
                    #    Eds[ktrace].data=lfilt(Eds[ktrace].data,bandpass,fsample,2)
                    #    Nds[ktrace].data=lfilt(Nds[ktrace].data,bandpass,fsample,2)
                    #    Zds[ktrace].data=lfilt(Zds[ktrace].data,bandpass,fsample,2)
                    #    #bandpass=None

def make_shakemap_slice(home,project_name,run_name,time_epi,GF_list,dt,tmax):
    '''
    Make xyz files with current ground velocity
    '''
    from numpy import genfromtxt,arange,sqrt,zeros,where,c_,savetxt
    from obspy import read
    from mudpy.forward import lowpass
    
    t=arange(0,tmax,dt)
    fcorner=0.5
    sta=genfromtxt(home+project_name+'/data/station_info/'+GF_list,usecols=0,dtype='S')
    lonlat=genfromtxt(home+project_name+'/data/station_info/'+GF_list,usecols=[1,2])
    for ksta in range(len(sta)):
        n=read(home+project_name+'/output/forward_models/'+run_name+'.'+sta[ksta]+'.vel.n')
        e=read(home+project_name+'/output/forward_models/'+run_name+'.'+sta[ksta]+'.vel.n')
        n[0].data=lowpass(n[0].data,fcorner,1./n[0].stats.delta,2)
        e[0].data=lowpass(n[0].data,fcorner,1./e[0].stats.delta,2)
        if ksta==0:
            h=n.copy()
        else:
            h+=n[0].copy()
        h[ksta].data=sqrt(n[0].data**2+e[0].data**2)
        h[ksta].trim(starttime=time_epi)
        print(h[ksta].stats.starttime)
    vout=zeros(len(lonlat))
    maxv=0
    for kt in range(len(t)):
        for ksta in range(len(sta)):
            i=where(h[ksta].times()==t[kt])[0]
            vout[ksta]=h[ksta].data[i]
        if vout.max()>maxv:
            maxv=vout.max()

if sort.lower()=='lon':
        j=argsort(lon[i])[::-1]
        i=i[j]
    elif sort.lower()=='lat':
        j=argsort(lat[i])[::-1] 
        i=i[j]
    nsta=len(i)
    fig, axarr = plt.subplots(nsta, 3)  
    for k in range(len(i)):
        n=read(datapath+sta[i[k]]+'.'+datasuffix+'.n')
        e=read(datapath+sta[i[k]]+'.'+datasuffix+'.e')
        u=read(datapath+sta[i[k]]+'.'+datasuffix+'.u')
        if lowpass!=None:
            fsample=1./e[0].stats.delta
            e[0].data=lfilter(e[0].data,lowpass,fsample,2)
            n[0].data=lfilter(n[0].data,lowpass,fsample,2)
            u[0].data=lfilter(u[0].data,lowpass,fsample,2)
        if decimate!=None:
            n[0]=stdecimate(n[0],decimate)
            e[0]=stdecimate(e[0],decimate)
            u[0]=stdecimate(u[0],decimate)
        if scale!=None:
            n[0].data=n[0].data/scale
            e[0].data=e[0].data/scale
            u[0].data=u[0].data/scale
        #Make plot
        if nsta>1:
            axn=axarr[k,0]
            axe=axarr[k,1]
            axu=axarr[k,2]
        else:
            axn=axarr[0]

e[0].stats.starttime=round_time(e[0].stats.starttime,dt)
        n[0].stats.starttime=round_time(n[0].stats.starttime,dt)
        u[0].stats.starttime=round_time(u[0].stats.starttime,dt)
        dvel=append(dvel,r_[n[0].data,e[0].data,u[0].data])
    
    #Tsunami
    kgf=3
    dtsun=array([])
    i=where(GF[:,kgf]==1)[0]
    for ksta in range(len(i)):
        if quiet==False:  
            print('Assembling tsunami waveforms from '+stations[i[ksta]]+' into data vector.')
        tsun=read(GFfiles[i[ksta],kgf])
        if tsunami_bandpass is not None: #Apply low pass filter to data
            fsample=1./tsun[0].stats.delta
            tsun[0].data=lfilt(tsun[0].data,tsunami_bandpass,fsample,2)
            
        dtsun=append(dtsun,tsun)
    
    #InSAR LOS
    kgf=4
    dlos=array([])  
    i=where(GF[:,kgf]==1)[0]
    for ksta in range(len(i)):
        if quiet==False:  
            print('Assembling InSAR LOS offsets from '+stations[i[ksta]]+' into data vector.')
        dtemp=genfromtxt(GFfiles[i[ksta],kgf])
        los=dtemp[0]
        dlos=append(dlos,los)          
    #Done, concatenate all, convert to column vector and exit
    d=concatenate([dx for dx in [dstatic,ddisp,dvel,dtsun,dlos] if dx.size > 0])
    D=zeros((d.shape[0],1))

datasuffix='kdisp'
        synthsuffix='disp'
    elif v_or_d.lower()=='v':
        kgf=1 #disp
        datasuffix='kvel'
        synthsuffix='vel'
    i=where(gf[:,kgf]==1)[0]
    for k in range(len(i)):
        print('Working on '+sta[i[k]])
        #Read data
        n=read(datapath+sta[i[k]]+'.'+datasuffix+'.n')
        e=read(datapath+sta[i[k]]+'.'+datasuffix+'.e')
        u=read(datapath+sta[i[k]]+'.'+datasuffix+'.u')
        if lowpass!=None:
            fsample=1./e[0].stats.delta
            e[0].data=lfilter(e[0].data,lowpass,fsample,10)
            n[0].data=lfilter(n[0].data,lowpass,fsample,10)
            u[0].data=lfilter(u[0].data,lowpass,fsample,10)
        if decimate!=None:
            n[0]=stdecimate(n[0],decimate)
            e[0]=stdecimate(e[0],decimate)
            u[0]=stdecimate(u[0],decimate)
        #Read synthetics
        nsyn=read(synthpath+run_name+'.'+run_number+'.'+sta[i[k]]+'.'+synthsuffix+'.n.sac')
        esyn=read(synthpath+run_name+'.'+run_number+'.'+sta[i[k]]+'.'+synthsuffix+'.e.sac')
        usyn=read(synthpath+run_name+'.'+run_number+'.'+sta[i[k]]+'.'+synthsuffix+'.u.sac')
        #What's the sampling rate correction?
        Fs=1./n[0].stats.delta
        Fsc=Fs/(2.*pi)
        #Compute coherence
        data=c_[n[0].data,nsyn[0].data].T
        fn,cn=tsa.cohere.coherence_regularized(data,epsilon=eps,alpha=al,

for k in range(len(i)):
        n=read(datapath+stations[i[k]]+'.'+datasuffix+'.n')
        e=read(datapath+stations[i[k]]+'.'+datasuffix+'.e')
        u=read(datapath+stations[i[k]]+'.'+datasuffix+'.u')
        ns=read(synthpath+run_name+'.'+run_number+'.'+stations[i[k]]+'.'+synthsuffix+'.n.sac')
        es=read(synthpath+run_name+'.'+run_number+'.'+stations[i[k]]+'.'+synthsuffix+'.e.sac')
        us=read(synthpath+run_name+'.'+run_number+'.'+stations[i[k]]+'.'+synthsuffix+'.u.sac')
        
        #Lowpass
        fsample=1./e[0].stats.delta
        e[0].data=lfilter(e[0].data,displacement_bandpass,fsample,2)
        n[0].data=lfilter(n[0].data,displacement_bandpass,fsample,2)
        u[0].data=lfilter(u[0].data,displacement_bandpass,fsample,2)
        es[0].data=lfilter(es[0].data,displacement_bandpass,fsample,2)
        ns[0].data=lfilter(ns[0].data,displacement_bandpass,fsample,2)
        us[0].data=lfilter(us[0].data,displacement_bandpass,fsample,2)
        
        if use_weights==True:
            nw=ones(n[0].stats.npts)*nweight[k]
            ew=ones(n[0].stats.npts)*eweight[k]
            uw=ones(n[0].stats.npts)*uweight[k]
            dtemp=r_[n[0].data/nw,e[0].data/ew,u[0].data/uw]
            dstemp=r_[ns[0].data/nw,es[0].data/ew,us[0].data/uw]  
        else:
            dtemp=r_[n[0].data,e[0].data,u[0].data]
            dstemp=r_[ns[0].data,es[0].data,us[0].data]
        
        if k==0:
            d_disp=dtemp.copy()
            ds_disp=dstemp.copy()
        else:

ddisp=append(ddisp,r_[n[0].data,e[0].data,u[0].data])
    
    #Velocities
    kgf=2
    dvel=array([])
    i=where(GF[:,kgf]==1)[0]
    for ksta in range(len(i)):
        if quiet==False:  
            print('Assembling acceleration waveforms from '+stations[i[ksta]]+' into data vector.')
        n=read(GFfiles[i[ksta],kgf]+'.n')
        e=read(GFfiles[i[ksta],kgf]+'.e')
        u=read(GFfiles[i[ksta],kgf]+'.u')
        if velocity_bandpass is not None: #Apply low pass filter to data
            fsample=1./n[0].stats.delta
            n[0].data=lfilt(n[0].data,velocity_bandpass,fsample,2)
            e[0].data=lfilt(e[0].data,velocity_bandpass,fsample,2)
            u[0].data=lfilt(u[0].data,velocity_bandpass,fsample,2)
        #Decimate
        if decimate!=None:
            n[0]=stdecimate(n[0],decimate)
            e[0]=stdecimate(e[0],decimate)
            u[0]=stdecimate(u[0],decimate)
        #Make sure they are rounded to a dt interval
        dt=e[0].stats.delta
        e[0].stats.starttime=round_time(e[0].stats.starttime,dt)
        n[0].stats.starttime=round_time(n[0].stats.starttime,dt)
        u[0].stats.starttime=round_time(u[0].stats.starttime,dt)
        dvel=append(dvel,r_[n[0].data,e[0].data,u[0].data])
    
    #Tsunami
    kgf=3
    dtsun=array([])

#Displacements
    kgf=1
    ddisp=array([])
    i=where(GF[:,kgf]==1)[0]
    for ksta in range(len(i)):
        if quiet==False:  
            print('Assembling displacement waveforms from '+str(stations[i[ksta]])+' into data vector.')
        #print(str(GFfiles[i[ksta],kgf]))
        n=read(GFfiles[i[ksta],kgf]+'.n')
        e=read(GFfiles[i[ksta],kgf]+'.e')
        u=read(GFfiles[i[ksta],kgf]+'.u')
        if displacement_bandpass is not None: #Apply low pass filter to data
            fsample=1./n[0].stats.delta
            n[0].data=lfilt(n[0].data,displacement_bandpass,fsample,2)
            e[0].data=lfilt(e[0].data,displacement_bandpass,fsample,2)
            u[0].data=lfilt(u[0].data,displacement_bandpass,fsample,2)
        #Decimate
        if decimate!=None:
            n[0]=stdecimate(n[0],decimate)
            e[0]=stdecimate(e[0],decimate)
            u[0]=stdecimate(u[0],decimate)
        #Make sure they are rounded to a dt interval
        dt=e[0].stats.delta
        e[0].stats.starttime=round_time(e[0].stats.starttime,dt)
        n[0].stats.starttime=round_time(n[0].stats.starttime,dt)
        u[0].stats.starttime=round_time(u[0].stats.starttime,dt)
        ddisp=append(ddisp,r_[n[0].data,e[0].data,u[0].data])
    
    #Velocities
    kgf=2
    dvel=array([])