From c0dd1e224ec81bd3c2ce44cc764b9510feb5c381 Mon Sep 17 00:00:00 2001 From: Diego Melgar Date: Thu, 19 Aug 2021 08:48:44 -0700 Subject: [PATCH] Added teleseismic modeling to fakequakes --- src/python/mudpy/fakequakes.py | 4 +- src/python/mudpy/forward.py | 163 +++++++++++++++++ .../mudpy/generate_ruptures_parallel.py | 4 +- src/python/mudpy/gmttools.py | 5 +- src/python/mudpy/parallel.py | 168 ++++++++++++++++++ src/python/mudpy/runslip.py | 113 +++++++++++- src/python/mudpy/view.py | 6 +- src/python/mudpy/viewFQ.py | 24 +-- 8 files changed, 469 insertions(+), 18 deletions(-) diff --git a/src/python/mudpy/fakequakes.py b/src/python/mudpy/fakequakes.py index f285758..04892c8 100644 --- a/src/python/mudpy/fakequakes.py +++ b/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] diff --git a/src/python/mudpy/forward.py b/src/python/mudpy/forward.py index dd4261e..867e7d9 100644 --- a/src/python/mudpy/forward.py +++ b/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() diff --git a/src/python/mudpy/generate_ruptures_parallel.py b/src/python/mudpy/generate_ruptures_parallel.py index 18c3227..435f76d 100755 --- a/src/python/mudpy/generate_ruptures_parallel.py +++ b/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 diff --git a/src/python/mudpy/gmttools.py b/src/python/mudpy/gmttools.py index 41e0f00..3c01c7d 100644 --- a/src/python/mudpy/gmttools.py +++ b/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]) diff --git a/src/python/mudpy/parallel.py b/src/python/mudpy/parallel.py index 92c8ee2..71a0a6b 100644 --- a/src/python/mudpy/parallel.py +++ b/src/python/mudpy/parallel.py @@ -579,6 +579,161 @@ def run_parallel_synthetics(home,project_name,station_file,model_name,integrate, DSdf.to_csv(green_path+'subfault'+num+'.DS.static.neu',sep='\t',index=False,header=False) + + + + +def run_parallel_teleseismics_green(home,project_name,time_epi,station_file,model_name,teleseismic_vel_mod,endtime,rank,size): + ''' + Use green functions and compute synthetics at stations for a single source + and multiple stations. This code makes an external system call to syn.c first it + will make the external call for the strike-slip component then a second externall + call will be made for the dip-slip component. The unit amount of moment is 1e15 + which corresponds to Mw=3.9333... + + IN: + source: 1-row numpy array containig informaiton aboutt he source, lat, lon, depth, etc... + station_file: File name with the station coordinates + green_path: Directopry where GFs are stored + model_file: File containing the Earth velocity structure + integrate: =0 if youw ant velocity waveforms, =1 if you want displacements + static: =0 if computing full waveforms, =1 if computing only the static field + subfault: String indicating the subfault being worked on + coord_type: =0 if problem is in cartesian coordinates, =1 if problem is in lat/lon + + OUT: + log: Sysytem standard output and standard error for log + ''' + + import os + import requests + from mudpy.forward import get_mu + from numpy import array,genfromtxt + from obspy import read + from mudpy.green import src2sta + + #What parameters are we using? + if rank==0: + out='''Running all processes with: + home = %s + project_name = %s + station_file = %s + model_name = %s + time_epi = %s + end_time = %s + ''' %(home,project_name,station_file,model_name,str(time_epi),str(endtime)) + print(out) + + #url for web request + url='http://service.iris.edu/irisws/syngine/1/query' + + #Read your corresponding source file + mpi_source=genfromtxt(home+project_name+'/data/model_info/mpi_source.'+str(rank)+'.fault') + + #Constant parameters + rakeDS=90 #90 is thrust, -90 is normal + rakeSS=0 #0 is left lateral, 180 is right lateral + + #Load structure + model_file=home+project_name+'/structure/'+model_name + structure=genfromtxt(model_file) + + for ksource in range(len(mpi_source)): + + source=mpi_source[ksource,:] + + #Parse the soruce information + xs=source[1] + ys=source[2] + zs=source[3] + zs_in_meters=int(zs*1000) + strike=source[4] + dip=source[5] + + #check longitude because iris only allows +/-180 + if xs > 180: + xs -= 360 + + ss_length=source[8] + ds_length=source[9] + + strdepth='%.4f' % zs + subfault=str(int(source[0])).rjust(4,'0') + + #Where to save dynamic waveforms + green_path=home+project_name+'/GFs/dynamic/'+model_name+"_"+strdepth+".sub"+subfault+"/" + + #check if folder exists if not make it + if os.path.isdir(green_path) == False: + os.mkdir(green_path) + + staname=genfromtxt(station_file,dtype="U",usecols=0) + if staname.shape==(): #Single staiton file + staname=array([staname]) + + #Compute distances and azimuths + d,az,lon_sta,lat_sta=src2sta(station_file,source,output_coordinates=True) + + #Get moment corresponding to 1 meter of slip on subfault + mu=get_mu(structure,zs) + Mo=mu*ss_length*ds_length*1.0 + + #Move to output folder + os.chdir(green_path) + print('Processor '+str(rank)+' is working on subfault '+str(int(source[0]))+' and '+str(len(d))+' stations ') + + #This is looping over "sites" + for ksta in range(len(d)): + + #cehck longitude for valid range + if lon_sta[ksta] > 180: + lon_sta[ksta] -= 360 + + #Form web request for SS syntehtic + parameters = {'model': teleseismic_vel_mod, + 'sourcelatitude':str(ys), + 'sourcelongitude':str(xs), + 'sourcedepthinmeters':str(zs_in_meters), + 'sourcedoublecouple':str(strike)+','+str(dip)+','+str(rakeSS)+','+str(Mo), + 'receiverlatitude':str(lat_sta[ksta]), + 'receiverlongitude':str(lon_sta[ksta]), + 'format':'miniseed', + 'origintime':str(time_epi), + 'endtime':str(endtime)} + + web_request = requests.get(url, params = parameters) + + #Filename for temporary write + out=green_path+staname[ksta]+'.subfault'+subfault+'.SS.mseed' + f=open(out,'wb') + f.write(web_request.content) + + + #Form web request for DS syntehtic + parameters = {'model': teleseismic_vel_mod, + 'sourcelatitude':str(ys), + 'sourcelongitude':str(xs), + 'sourcedepthinmeters':str(zs_in_meters), + 'sourcedoublecouple':str(strike)+','+str(dip)+','+str(rakeDS)+','+str(Mo), + 'receiverlatitude':str(lat_sta[ksta]), + 'receiverlongitude':str(lon_sta[ksta]), + 'format':'miniseed', + 'origintime':str(time_epi), + 'endtime':str(endtime)} + + web_request = requests.get(url, params = parameters) + + #Filename for temporary write + out=green_path+staname[ksta]+'.subfault'+subfault+'.DS.mseed' + f=open(out,'wb') + f.write(web_request.content) + + + + + + + def run_parallel_synthetics_mt3d(home,project_name,station_file,model_name,forceMT,mt,insar,rank,size): @@ -996,6 +1151,7 @@ def run_parallel_synthetics_mt3d(home,project_name,station_file,model_name,force elif insar=='False': insar=False run_parallel_green(home,project_name,station_file,model_name,dt,NFFT,static,dk,pmin,pmax,kmax,tsunami,insar,rank,size) + elif sys.argv[1]=='run_parallel_synthetics': #Parse command line arguments home=sys.argv[2] @@ -1026,6 +1182,18 @@ def run_parallel_synthetics_mt3d(home,project_name,station_file,model_name,force mu_okada=float(sys.argv[15]) run_parallel_synthetics(home,project_name,station_file,model_name,integrate,static,tsunami,time_epi,beta,custom_stf,impulse,rank,size,insar,okada,mu_okada) + + elif sys.argv[1]=='run_parallel_teleseismics_green': + home=sys.argv[2] + project_name=sys.argv[3] + time_epi=UTCDateTime(sys.argv[4]) + station_file=sys.argv[5] + model_name=sys.argv[6] + teleseismic_vel_mod=sys.argv[7] + endtime=sys.argv[8] + + run_parallel_teleseismics_green(home,project_name,time_epi,station_file,model_name,teleseismic_vel_mod,endtime,rank,size) + elif sys.argv[1]=='run_parallel_synthetics_mt3d': #Parse command line arguments home=sys.argv[2] diff --git a/src/python/mudpy/runslip.py b/src/python/mudpy/runslip.py index fc1c103..f3f73ef 100644 --- a/src/python/mudpy/runslip.py +++ b/src/python/mudpy/runslip.py @@ -245,6 +245,74 @@ def make_parallel_green(home,project_name,station_file,fault_name,model_name,dt, mpi=split(mpi) p=subprocess.Popen(mpi) p.communicate() + + + + + +def make_parallel_teleseismics_green(home,project_name,station_file,fault_name,model_name,teleseismic_vel_mod,time_epi,endtime,ncpus,hot_start=0): + ''' + This routine set's up the computation of GFs for each subfault to all stations. + The GFs are impulse sources, they don't yet depend on strike and dip. + + IN: + home: Home directory + project_name: Name fo the problem + station_file: File with coordinates of stations + fault_name: Name of fault file + model_Name: Name of Earth structure model file + dt: Desired sampling itnerval for waveforms + NFFT: No. of samples requested in waveform (must be power of 2) + static: =0 if computing full waveforms, =1 if computing only the static field + hot_start: =k if you want to start computations at k-th subfault, =0 to compute all + + + OUT: + Nothing + ''' + from numpy import arange,savetxt,genfromtxt + from os import path,makedirs,environ + from shlex import split + import subprocess + + + green_path=home+project_name+'/GFs/' + station_file=home+project_name+'/data/station_info/'+station_file + fault_file=home+project_name+'/data/model_info/'+fault_name + #Load source model for station-event distance computations + source=genfromtxt(fault_file) + + #Create all output folders + for k in range(len(source)): + + strdepth='%.4f' % source[k,3] + subfault=str(k+1).rjust(4,'0') + + subfault_folder=green_path+'dynamic/'+model_name+'_'+strdepth+'.sub'+subfault + if path.exists(subfault_folder)==False: + #It doesn't, make it, don't be lazy + makedirs(subfault_folder) + + #Create individual source files + for k in range(ncpus): + i=arange(k+hot_start,len(source),ncpus) + mpi_source=source[i,:] + fmt='%d\t%10.6f\t%10.6f\t%.8f\t%10.6f\t%10.6f\t%10.6f\t%10.6f\t%10.6f\t%10.6f' + savetxt(home+project_name+'/data/model_info/mpi_source.'+str(k)+'.fault',mpi_source,fmt=fmt) + + #Make mpi system call + print("MPI: Starting GFs computation on", ncpus, "CPUs\n") + mud_source=environ['MUD']+'/src/python/mudpy/' + + + mpi='mpiexec -n '+str(ncpus)+' python '+mud_source+'parallel.py run_parallel_teleseismics_green '+home+' '+project_name+' '+str(time_epi)+' '+station_file+' '+model_name+' '+teleseismic_vel_mod+' '+str(endtime) + print(mpi) + mpi=split(mpi) + p=subprocess.Popen(mpi) + p.communicate() + + + @@ -523,7 +591,50 @@ def inversionGFs(home,project_name,GF_list,tgf_file,fault_name,model_name, make_parallel_synthetics(home,project_name,station_file,fault_name,model_name,integrate,static,tsunami,beta,hot_start,time_epi,ncpus,custom_stf,impulse,insar) - + + + + +#Compute GFs for the ivenrse problem +def teleseismicGFs(home,project_name,GF_list_teleseismic,fault_name,model_name,teleseismic_vel_mod,time_epi,endtime,ncpus): + ''' + This routine will read a .gflist file and compute the required GF type for each station + ''' + from numpy import genfromtxt,shape,floor + from os import remove + + #Read in GFlist and decide what to compute + gf_file=home+project_name+'/data/station_info/'+GF_list_teleseismic + stations=genfromtxt(gf_file,usecols=0,dtype='U') + lonlat=genfromtxt(gf_file,usecols=[1,2,]) + fault_file=home+project_name+'/data/model_info/'+fault_name + source=genfromtxt(fault_file) + num_faults=shape(source)[0] + + if num_faults/ncpus < 2: + ncpus=int(floor(num_faults/2.)) + print('Cutting back to ' + str(ncpus) + ' cpus for ' + str(num_faults) + ' subfaults') + # GFs can be computed all at the same time + station_file='temp.sta' + try: + remove(home+project_name+'/data/station_info/'+station_file) #Cleanup + except: + pass + + print('Teleseismic GFs requested...') + f=open(home+project_name+'/data/station_info/'+station_file,'w') + for k in range(len(stations)): #Write temp .sta file + out=stations[k]+'\t'+repr(lonlat[k,0])+'\t'+repr(lonlat[k,1])+'\n' + f.write(out) + f.close() + + make_parallel_teleseismics_green(home,project_name,station_file,fault_name,model_name,teleseismic_vel_mod,time_epi,endtime,ncpus) + + + + + + def run_inversion(home,project_name,run_name,fault_name,model_name,GF_list,G_from_file,G_name,epicenter, rupture_speed,num_windows,reg_spatial,reg_temporal,nfaults,beta,decimate,bandpass, diff --git a/src/python/mudpy/view.py b/src/python/mudpy/view.py index 4882e6c..e599b06 100644 --- a/src/python/mudpy/view.py +++ b/src/python/mudpy/view.py @@ -2721,7 +2721,7 @@ def dtopo_slices(dtopo_file,fault_file,out): plt.close("all") -def plot_grd(grdfile,zlims,cmap,flip_lon=False,return_data=False,z_variable=None): +def plot_grd(grdfile,zlims,cmap,flip_lon=False,return_data=False,z_variable=None,contours=None): ''' Quick plot of any GMT grd file, this will only work for NETCDF4 files, i.e. if you use GMT5. If you are outdated and use NETCDF3 you can edit this @@ -2754,6 +2754,10 @@ def plot_grd(grdfile,zlims,cmap,flip_lon=False,return_data=False,z_variable=None plt.title(grdfile) plt.pcolormesh(X,Y,z,vmin=zlims[0],vmax=zlims[1],cmap=cmap) plt.colorbar() + + if contours != None: + plt.contour(X,Y,z,levels=coqntours) + plt.xlabel('Longitude') plt.ylabel('Latitude') plt.show() diff --git a/src/python/mudpy/viewFQ.py b/src/python/mudpy/viewFQ.py index 0cb7296..dce71ab 100644 --- a/src/python/mudpy/viewFQ.py +++ b/src/python/mudpy/viewFQ.py @@ -445,7 +445,7 @@ def plot_LW_scaling(home,project_name,run_name): return L,W,Mw -def analyze_sources(home,project_name,run_name,Mw_lims=[7.75,9.35],return_values=False): +def analyze_sources(rupture_folder,Mw_lims=[7.75,9.35],return_values=False): ''' Basic parameters of the sources ''' @@ -456,8 +456,8 @@ def analyze_sources(home,project_name,run_name,Mw_lims=[7.75,9.35],return_values from matplotlib.ticker import MultipleLocator plt.rcParams.update({'font.size': 12}) - logs=sort(glob(home+project_name+'/output/ruptures/'+run_name+'*.log')) - ruptures=sort(glob(home+project_name+'/output/ruptures/'+run_name+'*.rupt')) + logs=sort(glob(rupture_folder+'*.log')) + ruptures=sort(glob(rupture_folder+'*.rupt')) L=zeros(len(logs)) W=zeros(len(logs)) @@ -498,7 +498,7 @@ def analyze_sources(home,project_name,run_name,Mw_lims=[7.75,9.35],return_values plt.subplot(331) Mw_synth=arange(7.8,9.3,0.1) - plt.plot(Mw_synth,10**(-2.37+0.57*Mw_synth),c='k',lw=2) + # plt.plot(Mw_synth,10**(-2.37+0.57*Mw_synth),c='k',lw=2) plt.scatter(Mw_actual,L,marker='+') plt.xlim(Mw_lims) plt.ylim([70,1100]) @@ -521,7 +521,7 @@ def analyze_sources(home,project_name,run_name,Mw_lims=[7.75,9.35],return_values plt.subplot(332) - plt.plot(Mw_synth,10**(-1.86+0.46*Mw_synth),c='k',lw=2) + # plt.plot(Mw_synth,10**(-1.86+0.46*Mw_synth),c='k',lw=2) plt.scatter(Mw_actual,W,marker='+') plt.xlim(Mw_lims) plt.ylim([9,300]) @@ -574,11 +574,11 @@ def analyze_sources(home,project_name,run_name,Mw_lims=[7.75,9.35],return_values ax.tick_params(which='major',length=7,width=1) ax.tick_params(which='minor',length=4,width=1) #Allen Hayes scaling lines - plt.plot([7.0,9.6],[0.4,21.1],c='k',lw=2) - plt.plot([7.0,9.6],[0.4,7.87],c='r',lw=2) - plt.plot([7.0,9.6],[0.55,10.71],c='orange',lw=2) - plt.plot([7.0,9.6],[1.20,23.33],c='g',lw=2) - plt.plot([7.0,9.27],[1.08,39.16],c='violet',lw=2) + # plt.plot([7.0,9.6],[0.4,21.1],c='k',lw=2) + # plt.plot([7.0,9.6],[0.4,7.87],c='r',lw=2) + # plt.plot([7.0,9.6],[0.55,10.71],c='orange',lw=2) + # plt.plot([7.0,9.6],[1.20,23.33],c='g',lw=2) + # plt.plot([7.0,9.27],[1.08,39.16],c='violet',lw=2) plt.subplot(335) plt.scatter(Mw_actual,slip_max,marker='+') @@ -598,8 +598,8 @@ def analyze_sources(home,project_name,run_name,Mw_lims=[7.75,9.35],return_values ax.tick_params(which='major',length=7,width=1) ax.tick_params(which='minor',length=4,width=1) #Allen Hayes scaling lines - plt.plot([7.0,9.6],[1.07,75],c='k',lw=2) - plt.plot([7.0,9.6],[3.21,225],c='k',lw=2) + # plt.plot([7.0,9.6],[1.07,75],c='k',lw=2) + # plt.plot([7.0,9.6],[3.21,225],c='k',lw=2) plt.subplot(336) plt.scatter(Mw_actual,slip_stdev,marker='+')