Added teleseismic modeling to fakequakes

src/python/mudpy/fakequakes.py

@@ -800,6 +800,8 @@ def get_rise_times(M0,slip,fault_array,rise_time_depths,stoc_rake,rise_time_std=
 
 def get_rupture_onset(home,project_name,slip,fault_array,model_name,hypocenter,
         rise_time_depths,M0,velmod,sigma_rise_time=0.2,shear_wave_fraction_shallow=0.49,shear_wave_fraction_deep=0.8):
+
+    home,project_name,slip,fault_array,model_name,hypocenter,rise_time_depths,M0,velmod,shear_wave_fraction_shallow,shear_wave_fraction_deep
     '''
     Using a custom built tvel file ray trace from hypocenter to determine rupture
     onset times
@@ -834,6 +836,7 @@ def get_rupture_onset(home,project_name,slip,fault_array,model_name,hypocenter,
 #    rupture_multiplier[i]=slope*depth_to_top[i]+intercept
 
 
+
     #Get rupture speed shear-wave multipliers
     rupture_multiplier=zeros(len(vel))
     # Shallow 
@@ -849,7 +852,6 @@ def get_rupture_onset(home,project_name,slip,fault_array,model_name,hypocenter,
     rupture_multiplier[i]=slope*depth_to_top[i]+intercept
 
 
-
     #Perturb depths of the hypocenter so that faults at the same depth are not zero onset
     delta=0.00001
     i_same_as_hypo=where(fault_array[:,3]==hypocenter[2])[0]

src/python/mudpy/forward.py

@@ -324,6 +324,57 @@ def waveforms_fakequakes(home,project_name,fault_name,rupture_list,GF_list,
         #Write output
         write_fakequakes_waveforms(home,project_name,rupture_name,waveforms,GF_list,NFFT,time_epi,dt)
 
+
+
+
+def tele_waveforms(home,project_name,fault_name,rupture_list,GF_list_teleseismic,
+                model_name,run_name,G_from_file,G_name,source_time_function='dreger',zeta=0.2,
+                stf_falloff_rate=4.0,rupture_name=None,epicenter=None,time_epi=None,
+                hot_start=0,decimation_factor=1):
+    '''
+
+    '''
+    from numpy import genfromtxt,array
+    import datetime
+
+    print('Solving for kinematic problem(s)')
+    #Time for log file
+    now=datetime.datetime.now()
+    now=now.strftime('%b-%d-%H%M')
+
+    #load source names
+    if rupture_list==None: 
+        all_sources=array([rupture_name])  
+    else:
+        all_sources=genfromtxt(home+project_name+'/data/'+rupture_list,dtype='U')
+
+    #Load all synthetics
+    print('... loading all synthetics into memory')
+    Nss,Ess,Zss,Nds,Eds,Zds=load_fakequakes_tele_synthetics(home,project_name,fault_name,model_name,GF_list_teleseismic,G_from_file,G_name,decimation_factor)
+    print('... ... done')
+
+    Npoints=Nss[0].stats.npts
+    dt=Nss[0].stats.delta
+    #Now loop over rupture models
+    for ksource in range(hot_start,len(all_sources)):
+        print('... solving for source '+str(ksource)+' of '+str(len(all_sources)))
+        rupture_name=all_sources[ksource]
+        print(rupture_name)
+
+        if rupture_list!=None:
+            #Get epicentral time
+            epicenter,time_epi=read_fakequakes_hypo_time(home,project_name,rupture_name)
+            forward=False
+        else:
+            forward=True #This controls where we look for the rupture file
+
+        # Put in matrix
+        m,G=get_fakequakes_G_and_m(Nss,Ess,Zss,Nds,Eds,Zds,home,project_name,rupture_name,time_epi,GF_list_teleseismic,epicenter,Npoints,source_time_function,stf_falloff_rate,zeta,forward=forward)
+        # Solve
+        waveforms=G.dot(m)
+        #Write output
+        write_fakequakes_waveforms(home,project_name,rupture_name,waveforms,GF_list_teleseismic,Npoints,time_epi,dt)
+
 
 
 def waveforms_fakequakes_dynGF(home,project_name,fault_name,rupture_list,GF_list,dynamic_GFlist,dist_threshold,
@@ -880,6 +931,115 @@ def load_fakequakes_synthetics(home,project_name,fault_name,model_name,GF_list,G
 
 
 
+def load_fakequakes_tele_synthetics(home,project_name,fault_name,model_name,GF_list_teleseismic,G_from_file,G_name,decimation_factor):
+    ''''
+    Load the miniseed files with all the synthetics
+    '''
+    from numpy import genfromtxt,loadtxt
+    from obspy import read,Stream
+
+    if G_from_file==True: #load from file
+        Eds=read(home+project_name+'/GFs/matrices/'+G_name+'.Eds.tele.mseed')
+        Nds=read(home+project_name+'/GFs/matrices/'+G_name+'.Nds.tele.mseed')
+        Zds=read(home+project_name+'/GFs/matrices/'+G_name+'.Zds.tele.mseed')
+        Ess=read(home+project_name+'/GFs/matrices/'+G_name+'.Ess.tele.mseed')
+        Nss=read(home+project_name+'/GFs/matrices/'+G_name+'.Nss.tele.mseed')
+        Zss=read(home+project_name+'/GFs/matrices/'+G_name+'.Zss.tele.mseed')
+    else: #assemble G one data type at a time, just displacememnt right now
+        #Load station info
+        station_file=home+project_name+'/data/station_info/'+GF_list_teleseismic
+        staname=genfromtxt(station_file,dtype="U",usecols=0)
+        Nsta=len(staname)
+        #Load fault model
+        source=loadtxt(home+project_name+'/data/model_info/'+fault_name,ndmin=2)
+        Nfaults=source.shape[0] #Number of subfaults
+        kindex=0
+
+        for ksta in range(Nsta):
+
+            print('Reading green functions for station #'+str(ksta+1)+' of '+str(Nsta))
+
+            for kfault in range(Nfaults):
+
+                #Get subfault GF directory
+                nsub='sub'+str(int(source[kfault,0])).rjust(4,'0')
+                nfault='subfault'+str(int(source[kfault,0])).rjust(4,'0')
+                strdepth='%.4f' % source[kfault,3]
+                syn_path=home+project_name+'/GFs/dynamic/'+model_name+'_'+strdepth+'.'+nsub+'/'
+
+                #Get synthetics
+                if kfault==0 and ksta==0: #It's the first one, initalize stream object
+                    SS=read(syn_path+staname[ksta]+'.'+nfault+'.SS.mseed')
+                    DS=read(syn_path+staname[ksta]+'.'+nfault+'.DS.mseed')
+
+                    for trace in SS:
+
+                        if decimation_factor >1:
+                            trace.decimate(factor=decimation_factor)    
+
+
+                        if trace.stats.channel=='BXE':
+                            Ess=Stream(trace)
+                        if trace.stats.channel=='BXN':
+                            Nss=Stream(trace)
+                        if trace.stats.channel=='BXZ':
+                            Zss=Stream(trace)
+
+                    for trace in DS:
+
+                        if decimation_factor >1:
+                            trace.decimate(factor=decimation_factor)    
+
+                        if trace.stats.channel=='BXE':
+                            Eds=Stream(trace)
+                        if trace.stats.channel=='BXN':
+                            Nds=Stream(trace)
+                        if trace.stats.channel=='BXZ':
+                            Zds=Stream(trace)
+
+
+                else: #Just add to stream object
+
+                    SS=read(syn_path+staname[ksta]+'.'+nfault+'.SS.mseed')
+                    DS=read(syn_path+staname[ksta]+'.'+nfault+'.DS.mseed')
+
+                    for trace in SS:
+
+                        if decimation_factor >1:
+                            trace.decimate(factor=decimation_factor)  
+
+                        if trace.stats.channel=='BXE':
+                            Ess+=trace
+                        if trace.stats.channel=='BXN':
+                            Nss+=trace
+                        if trace.stats.channel=='BXZ':
+                            Zss+=trace
+
+                    for trace in DS:
+
+                        if decimation_factor >1:
+                            trace.decimate(factor=decimation_factor)  
+
+                        if trace.stats.channel=='BXE':
+                            Eds+=trace
+                        if trace.stats.channel=='BXN':
+                            Nds+=trace
+                        if trace.stats.channel=='BXZ':
+                            Zds+=trace
+                kindex+=1
+        print('Writting synthetics to miniSEED, hang on this might take a minute or two.')
+        Ess.write(home+project_name+'/GFs/matrices/'+G_name+'.Ess.tele.mseed',format='MSEED')
+        Nss.write(home+project_name+'/GFs/matrices/'+G_name+'.Nss.tele.mseed',format='MSEED')
+        Zss.write(home+project_name+'/GFs/matrices/'+G_name+'.Zss.tele.mseed',format='MSEED')
+        Eds.write(home+project_name+'/GFs/matrices/'+G_name+'.Eds.tele.mseed',format='MSEED')
+        Nds.write(home+project_name+'/GFs/matrices/'+G_name+'.Nds.tele.mseed',format='MSEED')
+        Zds.write(home+project_name+'/GFs/matrices/'+G_name+'.Zds.tele.mseed',format='MSEED')
+    return Nss,Ess,Zss,Nds,Eds,Zds
+
+
+
+
+
 def get_fakequakes_G_and_m(Nss,Ess,Zss,Nds,Eds,Zds,home,project_name,rupture_name,time_epi,GF_list,epicenter,NFFT,
                 source_time_function,stf_falloff_rate,zeta=0.2,forward=False):
     '''
@@ -938,6 +1098,9 @@ def get_fakequakes_G_and_m(Nss,Ess,Zss,Nds,Eds,Zds,home,project_name,rupture_nam
 
         for ksource in range(len(i_non_zero)):
 
+            if ksource % 50 == 0:
+                print('... ... ... subfault %d of %d' % (ksource,len(i_non_zero)))
+
             #Get synthetics
             nss=Nss[read_start+i_non_zero[ksource]].copy()
             ess=Ess[read_start+i_non_zero[ksource]].copy()

src/python/mudpy/generate_ruptures_parallel.py

@@ -226,7 +226,9 @@ def run_parallel_generate_ruptures(home,project_name,run_name,fault_name,slab_na
             else:
                 hypocenter=whole_fault[shypo,1:4]
 
-            t_onset=fakequakes.get_rupture_onset(home,project_name,slip,fault_array,model_name,hypocenter,rise_time_depths,M0,velmod,shear_wave_fraction_shallow,shear_wave_fraction_deep)
+            t_onset=fakequakes.get_rupture_onset(home,project_name,slip,fault_array,model_name,hypocenter,rise_time_depths,
+                                                 M0,velmod,shear_wave_fraction_shallow=shear_wave_fraction_shallow,
+                                                 shear_wave_fraction_deep=shear_wave_fraction_deep)
             fault_out[:,12]=0
             fault_out[ifaults,12]=t_onset
 

src/python/mudpy/gmttools.py

@@ -718,7 +718,7 @@ def read_neic_param(fault_file):
 
 
 
-def triangular_rupt_2_gmt(meshfile,slipfile,outfile,kinematic_out_folder=None):
+def triangular_rupt_2_gmt(meshfile,slipfile,outfile,kinematic_out_folder=None,percentage=0):
 
     '''
     DM Note: Modified from Brendan's script because he refused to do a pull request :)
@@ -815,11 +815,12 @@ def triangular_rupt_2_gmt(meshfile,slipfile,outfile,kinematic_out_folder=None):
     #Total slip model
     moment = 0
     fso = open(outfile,'w')
+    slip_threshold=(percentage/100)*TOTS.max()
     for i in range(0, numpy.amax(MESN)):
         a1 = numpy.where(MESN[i] == INVN)[0]
         totslip = numpy.sum(TOTS[a1])
 #        print (i+1,totslip*100)
-        if (totslip >= 0.0):
+        if (totslip >= slip_threshold):
             moment = moment+FA[i]*1000*1000*numpy.mean(RIG[a1])*totslip
             lon1 = "{0:.4f}".format(meshlon1[i])
             lon2 = "{0:.4f}".format(meshlon2[i])