fixed a bug in to_egf()

xtyangpsp · Apr 5, 2022 · 0be3a46 · 0be3a46
1 parent f9aff4d
commit 0be3a46
Show file tree

Hide file tree

Showing 2 changed files with 9 additions and 6 deletions.
diff --git a/changes.txt b/changes.txt
@@ -10,6 +10,7 @@ TYPES:
 1. Added CorrData.shaping() to shape the data with wavelet.
 2. Added CorrData.save() to wrap saving functions.
 3. In CorrData.to_asdf(), save stack_method.
+4. Fixed a bug in CorrData.to_egf() where the zero lag was not handled correctly. The negative side was wrong.
 
 UTILS:
 1. Added gaussian() and ricker() as the shaping wavelets.

diff --git a/seisgo/types.py b/seisgo/types.py
@@ -909,27 +909,29 @@ def to_egf(self,taper_frac=0.01,taper_maxlen=10,verbose=False):
             egf=np.zeros(self.data.shape,dtype=self.data.dtype)
             if self.substack:
                 if side.lower()=="a":
-                    nhalfpoint=np.int(self.data.shape[1]/2)
+                    nhalfpoint=int(self.data.shape[1]/2)
                     #positive side
                     egf[:,nhalfpoint:]=utils.taper(-1.0*np.gradient(self.data[:,nhalfpoint:],axis=1)/dt,
                                                     fraction=taper_frac,maxlen=taper_maxlen)
                     #negative side
-                    egf[:,:nhalfpoint+1]=np.flip(utils.taper(np.gradient(np.flip(self.data[:,:nhalfpoint+1],axis=1),
+                    egf[:,:nhalfpoint+1]=np.flip(utils.taper(-1.0*np.gradient(np.flip(self.data[:,:nhalfpoint+1],axis=1),
                                                     axis=1)/dt,fraction=taper_frac,maxlen=taper_maxlen),axis=1)
-                    egf[:,[0,nhalfpoint,-1]]=0
+                    egf[:,[0,-1]]=0
+                    egf[:,nhalfpoint]=np.mean(egf[:,nhalfpoint-1:nhalfpoint+1],axis=1)
                 else:
                     egf=utils.taper(-1.0*np.gradient(self.data,axis=1)/dt,
                                                     fraction=taper_frac,maxlen=taper_maxlen)
             else:
                 if side.lower()=="a":
-                    nhalfpoint=np.int(self.data.shape[0]/2)
+                    nhalfpoint=int(self.data.shape[0]/2)
                     #positive side
                     egf[nhalfpoint:]=utils.taper(-1.0*np.gradient(self.data[nhalfpoint:])/dt,
                                                 fraction=taper_frac,maxlen=taper_maxlen)
                     #negative side
-                    egf[:nhalfpoint+1]=np.flip(utils.taper(np.gradient(np.flip(self.data[:nhalfpoint+1]))/dt,
+                    egf[:nhalfpoint+1]=np.flip(utils.taper(-1.0*np.gradient(np.flip(self.data[:nhalfpoint+1]))/dt,
                                                 fraction=taper_frac,maxlen=taper_maxlen))
-                    egf[[0,nhalfpoint,-1]]=0
+                    egf[0,-1]=0
+                    egf[nhalfpoint]=np.mean(egf[nhalfpoint-1:nhalfpoint+1])
                 else:
                     egf=utils.taper(-1.0*np.gradient(self.data)/dt,
                                                 fraction=taper_frac,maxlen=taper_maxlen)