How to use the green.stdecimate function in green

e[0].data=ds[kinsert+npts:kinsert+2*npts]
            u[0].data=ds[kinsert+2*npts:kinsert+3*npts]
            kinsert+=3*npts
            n.write(home+project_name+'/output/inverse_models/waveforms/'+run_name+'.'+num+'.'+sta+'.disp.n.sac',format='SAC')
            e.write(home+project_name+'/output/inverse_models/waveforms/'+run_name+'.'+num+'.'+sta+'.disp.e.sac',format='SAC')
            u.write(home+project_name+'/output/inverse_models/waveforms/'+run_name+'.'+num+'.'+sta+'.disp.u.sac',format='SAC')
    #Velocity
    kgf=2
    i=where(GF[:,kgf]==1)[0]
    if len(i)>0:
        for ksta in range(len(i)):
            sta=stations[i[ksta]]
            n=read(GFfiles[i[ksta],kgf]+'.n')
            e=read(GFfiles[i[ksta],kgf]+'.e')
            u=read(GFfiles[i[ksta],kgf]+'.u')
            n=stdecimate(n,decimate)
            e=stdecimate(e,decimate)
            u=stdecimate(u,decimate)
            npts=n[0].stats.npts
            n[0].data=ds[kinsert:kinsert+npts]
            e[0].data=ds[kinsert+npts:kinsert+2*npts]
            u[0].data=ds[kinsert+2*npts:kinsert+3*npts]
            kinsert+=3*npts
            n.write(home+project_name+'/output/inverse_models/waveforms/'+run_name+'.'+num+'.'+sta+'.vel.n.sac',format='SAC')
            e.write(home+project_name+'/output/inverse_models/waveforms/'+run_name+'.'+num+'.'+sta+'.vel.e.sac',format='SAC')
            u.write(home+project_name+'/output/inverse_models/waveforms/'+run_name+'.'+num+'.'+sta+'.vel.u.sac',format='SAC')

u[0].data=ds[kinsert+2*npts:kinsert+3*npts]
            kinsert+=3*npts
            n.write(home+project_name+'/output/inverse_models/waveforms/'+run_name+'.'+num+'.'+sta+'.disp.n.sac',format='SAC')
            e.write(home+project_name+'/output/inverse_models/waveforms/'+run_name+'.'+num+'.'+sta+'.disp.e.sac',format='SAC')
            u.write(home+project_name+'/output/inverse_models/waveforms/'+run_name+'.'+num+'.'+sta+'.disp.u.sac',format='SAC')
    #Velocity
    kgf=2
    i=where(GF[:,kgf]==1)[0]
    if len(i)>0:
        for ksta in range(len(i)):
            sta=stations[i[ksta]]
            n=read(GFfiles[i[ksta],kgf]+'.n')
            e=read(GFfiles[i[ksta],kgf]+'.e')
            u=read(GFfiles[i[ksta],kgf]+'.u')
            n=stdecimate(n,decimate)
            e=stdecimate(e,decimate)
            u=stdecimate(u,decimate)
            npts=n[0].stats.npts
            n[0].data=ds[kinsert:kinsert+npts]
            e[0].data=ds[kinsert+npts:kinsert+2*npts]
            u[0].data=ds[kinsert+2*npts:kinsert+3*npts]
            kinsert+=3*npts
            n.write(home+project_name+'/output/inverse_models/waveforms/'+run_name+'.'+num+'.'+sta+'.vel.n.sac',format='SAC')
            e.write(home+project_name+'/output/inverse_models/waveforms/'+run_name+'.'+num+'.'+sta+'.vel.e.sac',format='SAC')
            u.write(home+project_name+'/output/inverse_models/waveforms/'+run_name+'.'+num+'.'+sta+'.vel.u.sac',format='SAC')

ddisp=array([])
    i=where(GF[:,kgf]==1)[0]
    for ksta in range(len(i)):
        print 'Assembling displacement 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 lowpass!=None: #Apply low pass filter to data
            fsample=1./n[0].stats.delta
            n[0].data=lfilt(n[0].data,lowpass,fsample,9)
            e[0].data=lfilt(e[0].data,lowpass,fsample,9)
            u[0].data=lfilt(u[0].data,lowpass,fsample,9)
        #Decimate
        if decimate!=None:
            n=stdecimate(n,decimate)
            e=stdecimate(e,decimate)
            u=stdecimate(u,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([])
    i=where(GF[:,kgf]==1)[0]
    for ksta in range(len(i)):
        print 'Assembling velocity 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')

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,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=stdecimate(n,decimate)
            e=stdecimate(e,decimate)
            u=stdecimate(u,decimate)
        ns=read(synthpath+run_name+'.'+run_number+'.'+sta[i[k]]+'.'+synthsuffix+'.n.sac')
        es=read(synthpath+run_name+'.'+run_number+'.'+sta[i[k]]+'.'+synthsuffix+'.e.sac')
        us=read(synthpath+run_name+'.'+run_number+'.'+sta[i[k]]+'.'+synthsuffix+'.u.sac')
        #Make plot
        axn=axarr[k,0]
        axe=axarr[k,1]
        axu=axarr[k,2]
        axn.plot(n[0].times(),n[0].data,'k',ns[0].times(),ns[0].data,'r')
        axn.grid(which='both')
        axe.plot(e[0].times(),e[0].data,'k',es[0].times(),es[0].data,'r')
        axe.grid(which='both')
        axu.plot(u[0].times(),u[0].data,'k',us[0].times(),us[0].data,'r')
        axu.grid(which='both')
        axe.yaxis.set_ticklabels([])
        axu.yaxis.set_ticklabels([])
        axe.set_xlim(t_lim)

eds=tshift(eds,tdelay[kfault])
                    nds=tshift(nds,tdelay[kfault])
                    zds=tshift(zds,tdelay[kfault])
                    if decimate!=None: 
                        ess=stdecimate(ess,decimate)
                        nss=stdecimate(nss,decimate)
                        zss=stdecimate(zss,decimate)    
                        eds=stdecimate(eds,decimate)
                        nds=stdecimate(nds,decimate)
                        zds=stdecimate(zds,decimate)
                    if decimate!=None and kfault==0:#Only need to do this once
                        #Load raw data, this will be used to trim the GFs
                        edata=read(datafiles[ksta]+'.e')
                        ndata=read(datafiles[ksta]+'.n')
                        udata=read(datafiles[ksta]+'.u')
                        edata=stdecimate(edata,decimate)
                        ndata=stdecimate(ndata,decimate)
                        udata=stdecimate(udata,decimate)
                        #How many points left in the tiem series
                        npts=edata[0].stats.npts
                        print "... "+str(npts)+" data points left over after decimation"
                    #Now time align stuff (This is where we might have some timing issues, check later, consider re-sampling to data time vector)                                       
                    ess=resample_to_data(ess,edata)
                    ess=prep_synth(ess,edata)
                    nss=resample_to_data(nss,ndata)
                    nss=prep_synth(nss,ndata)
                    zss=resample_to_data(zss,udata)
                    zss=prep_synth(zss,udata)
                    eds=resample_to_data(eds,edata)
                    eds=prep_synth(eds,edata)
                    nds=resample_to_data(nds,ndata)
                    nds=prep_synth(nds,ndata)

'''
    
    from obspy import read
    from numpy import zeros
    from green import stdecimate
    
    npts=0
    if decimate==None:
        decimate=1
    for k in range(len(datafiles)):
        e=read(datafiles[k]+'.e')
        n=read(datafiles[k]+'.n')
        u=read(datafiles[k]+'.u')
        if e[0].stats.npts==n[0].stats.npts==u[0].stats.npts:
            if decimate!=None:
                e=stdecimate(e,decimate)
            npts+=e[0].stats.npts
        else:
            print str(e[0].stats.npts)+' pts in east component'
            print str(n[0].stats.npts)+' pts in north component'
            print str(u[0].stats.npts)+' pts in up component'
            print 'ERROR: The 3 components of data are not the same length'
            return 'Error in forming G'
    G=zeros([3*npts,nfaults*2])
    return G            

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,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=stdecimate(n,decimate)
            e=stdecimate(e,decimate)
            u=stdecimate(u,decimate)
        ns=read(synthpath+run_name+'.'+run_number+'.'+sta[i[k]]+'.'+synthsuffix+'.n.sac')
        es=read(synthpath+run_name+'.'+run_number+'.'+sta[i[k]]+'.'+synthsuffix+'.e.sac')
        us=read(synthpath+run_name+'.'+run_number+'.'+sta[i[k]]+'.'+synthsuffix+'.u.sac')
        #Make plot
        axn=axarr[k,0]
        axe=axarr[k,1]
        axu=axarr[k,2]
        axn.plot(n[0].times(),n[0].data,'k',ns[0].times(),ns[0].data,'r')
        axn.grid(which='both')
        axe.plot(e[0].times(),e[0].data,'k',es[0].times(),es[0].data,'r')
        axe.grid(which='both')
        axu.plot(u[0].times(),u[0].data,'k',us[0].times(),us[0].data,'r')
        axu.grid(which='both')
        axe.yaxis.set_ticklabels([])

kgf=1
    ddisp=array([])
    i=where(GF[:,kgf]==1)[0]
    for ksta in range(len(i)):
        print 'Assembling displacement 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 lowpass!=None: #Apply low pass filter to data
            fsample=1./n[0].stats.delta
            n[0].data=lfilt(n[0].data,lowpass,fsample,9)
            e[0].data=lfilt(e[0].data,lowpass,fsample,9)
            u[0].data=lfilt(u[0].data,lowpass,fsample,9)
        #Decimate
        if decimate!=None:
            n=stdecimate(n,decimate)
            e=stdecimate(e,decimate)
            u=stdecimate(u,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([])
    i=where(GF[:,kgf]==1)[0]
    for ksta in range(len(i)):
        print 'Assembling velocity waveforms from '+stations[i[ksta]]+' into data vector.'
        n=read(GFfiles[i[ksta],kgf]+'.n')
        e=read(GFfiles[i[ksta],kgf]+'.e')

ess[0].data=lfilt(ess[0].data,lowpass,fsample,9)
                        nss[0].data=lfilt(nss[0].data,lowpass,fsample,9)
                        zss[0].data=lfilt(zss[0].data,lowpass,fsample,9)
                        eds[0].data=lfilt(eds[0].data,lowpass,fsample,9)
                        nds[0].data=lfilt(nds[0].data,lowpass,fsample,9)
                        zds[0].data=lfilt(zds[0].data,lowpass,fsample,9)
                    #Time shift them according to subfault rupture time, zero pad, round to dt interval,decimate
                    #and extend to maximum time
                    ess=tshift(ess,tdelay[kfault])
                    nss=tshift(nss,tdelay[kfault])
                    zss=tshift(zss,tdelay[kfault])
                    eds=tshift(eds,tdelay[kfault])
                    nds=tshift(nds,tdelay[kfault])
                    zds=tshift(zds,tdelay[kfault])
                    if decimate!=None: 
                        ess=stdecimate(ess,decimate)
                        nss=stdecimate(nss,decimate)
                        zss=stdecimate(zss,decimate)    
                        eds=stdecimate(eds,decimate)
                        nds=stdecimate(nds,decimate)
                        zds=stdecimate(zds,decimate)
                    if decimate!=None and kfault==0:#Only need to do this once
                        #Load raw data, this will be used to trim the GFs
                        edata=read(datafiles[ksta]+'.e')
                        ndata=read(datafiles[ksta]+'.n')
                        udata=read(datafiles[ksta]+'.u')
                        edata=stdecimate(edata,decimate)
                        ndata=stdecimate(ndata,decimate)
                        udata=stdecimate(udata,decimate)
                        #How many points left in the tiem series
                        npts=edata[0].stats.npts
                        print "... "+str(npts)+" data points left over after decimation"