Merge pull request #237 from SeismoFelix/master

SPECFEM3D_SAC client support depth gridsearch
uafgeotools · Nov 9, 2023 · c759d16 · c759d16
2 parents 0641750 + ec96aea
commit c759d16
Show file tree

Hide file tree

Showing 3 changed files with 29 additions and 8 deletions.
diff --git a/data/examples/SPECFEM3D_SAC.tgz b/data/examples/SPECFEM3D_SAC.tgz
diff --git a/mtuq/io/clients/SPECFEM3D_SAC.py b/mtuq/io/clients/SPECFEM3D_SAC.py
@@ -4,6 +4,7 @@
 import numpy as np
 
 from obspy.core import Stream
+from os.path import join
 from mtuq.greens_tensor.SPECFEM3D import GreensTensor 
 from mtuq.io.clients.base import Client as ClientBase
 from mtuq.util.signal import resample
@@ -68,7 +69,8 @@ class Client(ClientBase):
     """
 
     def __init__(self, path_or_url=None, model=None, 
-                 include_mt=True, include_force=False):
+                 include_mt=True, include_force=False,
+                 units='km'):
 
         self.path = path_or_url
 
@@ -79,6 +81,13 @@ def __init__(self, path_or_url=None, model=None,
         self.include_mt = include_mt
         self.include_force = include_force
 
+        if not units:
+            pass
+        elif units.lower() not in ['','m','km']:
+            raise ValueError
+            units = units.lower()
+        self.units = units
+
 
     def get_greens_tensors(self, stations=[], origins=[], verbose=False):
         """ Extracts Green's tensors
@@ -105,29 +114,40 @@ def _get_greens_tensor(self, station=None, origin=None):
         stream = Stream()
         stream.id = station.id
 
+        # optionally, append depth subdirectory to path
+        path = self.path
+        if hasattr(origin, 'depth_in_m'):
+           if self.units=='m':
+               subdir = str(int(np.ceil(origin.depth_in_m)))
+           elif self.units=='km':
+               subdir = str(int(np.ceil(origin.depth_in_m/1000.)))
+           else:
+               raise ValueError
+           path = join(path, subdir)
+
         # check if Green's functions exist for given station code
-        if _exists(self.path+'/'+station.id+'.*.sac'):
+        if _exists(path+'/'+station.id+'.*.sac'):
             prefix = station.id
 
         else:
             print('No Green\'s functions found for "%s"\n\n'
                   'Trying other codes instead:'
                   % station.id)
 
-            prefix = _try_wildcards(self.path, station)
-            
+            prefix = _try_wildcards(path, station)
+
         if self.include_mt:
             for suffix in EXT_MT:
                 trace = obspy.read(
-                    self.path+'/'+prefix+'.'+suffix+'.sac', format='sac')[0]
+                    path+'/'+prefix+'.'+suffix+'.sac', format='sac')[0]
                 trace.stats.channel = suffix
                 trace.stats._component = suffix[0]
                 stream += trace
 
         if self.include_force:
             for suffix in EXT_FORCE:
                 trace = obspy.read(
-                    self.path+'/'+prefix+'.'+suffix+'.sac', format='sac')[0]
+                    path+'/'+prefix+'.'+suffix+'.sac', format='sac')[0]
                 trace.stats.channel = suffix
                 trace.stats._component = suffix[0]
                 stream += trace
@@ -178,7 +198,7 @@ def _try_wildcards(path, station):
 
     wildcards = [ 
         station.network+'.'+station.station+'.*',
-        '*'            +'.'+station.station+'.*',
+        '*'+'.'+station.station+'.*',
         ] 
 
     for wildcard in wildcards:

diff --git a/tests/test_greens_SPECFEM3D_SAC.py b/tests/test_greens_SPECFEM3D_SAC.py
@@ -94,7 +94,7 @@
         'time': '2009-04-07T20:12:55.000000Z',
         'latitude': 61.454200744628906,
         'longitude': -149.7427978515625,
-        'depth_in_m': 33033.599853515625,
+        'depth_in_m': 33000.0,
         })
 
 
@@ -148,6 +148,7 @@
     best_mt = grid.get(idx)
     lune_dict = grid.get_dict(idx)
 
+
     if run_figures:
 
         plot_data_greens2(event_id+'DC_waveforms.png',