include ray path density calculation

src/VoroTomo_2D_MPI.py

@@ -39,14 +39,14 @@
 ##########################################
 
 # define the absolute path of data and project
-rootpath = '/Users/chengxinjiang/Documents/Github/VoroTomo_SW/example/real_data'
-data = glob.glob(os.path.join(rootpath,'dispersion/LVC_snr8_2wl_8s.dat'))
+rootpath = '/Users/chengxin/Documents/ANU/St_Helens/Voro_Tomo/synthetic/dispersion/Love'
+data = glob.glob(os.path.join(rootpath,'syn_selection_*.dat'))
 
 # useful parameters for location and inversion
-ncell = 300                                                                 # number of Voronoi cells for the target region
-nsets = 300                                                                 # number of realizations
-latmin,dlat,nlat = 35.5,0.05,110                                            # latitude range of the target region
-lonmin,dlon,nlon = -122.5,0.05,100                                          # longitude range of the target region
+ncell = 350                                                                 # number of Voronoi cells for the target region
+nsets = 50                                                                 # number of realizations
+latmin,dlat,nlat = 44.5,0.05,76                                             # latitude range of the target region
+lonmin,dlon,nlon = -124.8,0.05,124                                          # longitude range of the target region
 velall           = np.zeros(((nlon+1)*(nlat+1),nsets))                      # vector of the final model averaged over (nsets) times of realizations
 geogrid = {'latmin':latmin,'lonmin':lonmin,                                 # assemble geographic info into a dict
             'dlat':dlat,'dlon':dlon,
@@ -60,8 +60,6 @@
 
 if rank==0:
     # make sub-directory to save intermediate files
-    if not os.path.isdir(os.path.join(rootpath,'figures')):
-        os.mkdir(os.path.join(rootpath,'figures'))
     if not os.path.isdir(os.path.join(rootpath,'iterations')):
         os.mkdir(os.path.join(rootpath,'iterations'))
 
@@ -86,10 +84,11 @@
 
         # define the size of subsampling
         obs_all = len(dfile)
-        obs_sub = int(np.floor(obs_all*0.7))
+        obs_sub = int(np.floor(obs_all*0.8))
 
         # read traveltime info
-        ave     = float("{:.3f}".format(np.mean(dfile['vel'])))                     # averaged velocity from the observation
+        #ave     = float("{:.3f}".format(np.mean(dfile['vel'])))                     # averaged velocity from the observation
+        ave     = 3
         sphgrid = subf.construct_grid(geogrid,ave)                                  # spherical grid data and initial velocity model
 
         # generate parameters
@@ -167,7 +166,7 @@
         np.savez_compressed(tname,vel=vel)
 
         # plot each inversion results
-        subf.plot_tomo(geogrid,vel,rootpath,iset)
+        #subf.plot_tomo(geogrid,vel,rootpath,iset)
 
 ############################################################
 ######### RUN construct_final_model_residuals NEXT #########
@@ -176,4 +175,4 @@
 # syn time to exit
 comm.barrier()
 if rank == 0:
-    sys.exit()
+    sys.exit()

src/construct_final_model_residula.py

@@ -1,5 +1,7 @@
 import os 
 import glob 
+import sys
+import scipy
 import numpy as np 
 import pandas as pd 
 import subfunction as subf 
@@ -19,17 +21,20 @@
 ###################################
 
 # define absolute path
-rootpath = '/Users/chengxinjiang/Documents/ANU/Voro_Tomo/synthetic'
+#rootpath = '/Users/chengxin/Documents/ANU/St_Helens/Voro_Tomo/synthetic/Oct4_WithNoise_0.5s_spacing2/Love'
+rootpath = "/Users/chengxin/Documents/ANU/St_Helens/Voro_Tomo/real_data_inversion_Cascadia/Love_phase/roundII"
 
 # control parameter
-cal_resid  = True                                                           # calculate the residual of initial/final model
-plot_final = False                                                          # plot the average model
-plot_each  = False                                                          # plot tomo at each iteration 
-real_data  = False                                                          # whether this is for real data or not
+cal_raypath = True                                                           # calculate the ray path of the averaged model
+cal_resid  = False                                                           # calculate the residual of initial/final model
+plot_final = False                                                           # plot the average model
+plot_each  = False                                                           # plot tomo at each iteration 
+write_each = False                                                           # write each model into txt file
+real_data  = False                                                           # whether this is for real data or not
 
 # several useful parameters
-pper    = [2,3,4,5,6,8,10,12,14,16,18,20,24,28]                                   # all period range
-nset    = 300                                                               # number of inversions
+pper    = [3,4,5,6,8,10,12,14,16,18,20,24,28,32,36,40]                                    # all period range
+nset    = 70                                                                 # number of inversions
 
 # check whether dir exists or not
 if not os.path.isdir(os.path.join(rootpath,'models')):
@@ -40,24 +45,38 @@
     os.mkdir(os.path.join(rootpath,'figures'))
 
 # define geographic information
-latmin,dlat,nlat = 35.5,0.02,100                                            # latitude range of the target region
-lonmin,dlon,nlon = -120.0,0.02,100                                          # longitude range of the target region
+latmin,dlat,nlat = 24.5,0.2,123                                            # latitude range of the target region
+lonmin,dlon,nlon = -124.8,0.2,289                                          # longitude range of the target region
+
+latmin,dlat,nlat = 44.5,0.02,190                                            # latitude range of the target region
+lonmin,dlon,nlon = -124.8,0.02,310                                          # longitude range of the target region
 geogrid = {'latmin':latmin,'lonmin':lonmin,                                 # assemble geographic info into a dict
             'dlat':dlat,'dlon':dlon,
             'nlat':nlat,'nlon':nlon}
 
+# grid for output files
+latgrid = np.linspace(latmin+nlat*dlat,latmin,nlat+1) 
+longrid = np.linspace(lonmin,lonmin+nlon*dlon,nlon+1)
+lon,lat = np.meshgrid(longrid,latgrid)
 
 #######################################
 ######## LOOPS THROUGH PERIODS ########
 #######################################
 
 # loop through all periods
 for iper in range(len(pper)):
+    tdlat,tnlat = 0.1,38
+    tdlon,tnlon = 0.1,62
+    tlatgrid = np.linspace(latmin+tnlat*tdlat,latmin,tnlat+1)
+    tlongrid = np.linspace(lonmin,lonmin+tnlon*tdlon,tnlon+1)
+    tlon,tlat = np.meshgrid(tlongrid,tlatgrid)
+    tcount = np.zeros(tlon.size,dtype=np.int16)
     per = pper[iper]
     velall = np.zeros(((nlon+1)*(nlat+1),nset))                             # vector of the final model averaged over (nsets) times of realizations
 
     # find the observational data
-    dfile = os.path.join(rootpath,'dispersion/data_snr8_'+str(per)+'s.dat')
+    #dfile = os.path.join(rootpath,'dispersion/Love_phase/new_L{0:02d}_USANT15.txt'.format(per))
+    dfile = os.path.join(rootpath,'../dispersion/selection_'+str(per)+'s.dat')
     if not os.path.isfile(dfile):
         raise ValueError('double check! cannot find %s'%dfile)
 
@@ -69,65 +88,112 @@
         ave = 3.0
 
     # load the velocity model
-    allfiles = glob.glob(os.path.join(rootpath,'iterations/'+str(per)+'s_*.npz'))
+    allfiles = sorted(glob.glob(os.path.join(rootpath,'iterations/{0:d}s_*.npz'.format(per))))
     if not len(allfiles):
         raise ValueError('no model generated! double check')
 
     # load all models into a big matrix
     ii = 0
-    for ifile in allfiles:
+    for ifile in allfiles[:nset]:
         try:
             tvel = np.load(ifile)['vel']
             if plot_each:
                 subf.plot_tomo(geogrid,-tvel*ave*ave+ave,rootpath,ii)
+            if write_each:
+                tvel1 = -tvel*ave*ave+ave
+                # format the data column
+                fout0 = open(rootpath+'/models/new_tomo_'+str(per)+'s_'+str(ii)+'m.txt','w')
+                for jj in range(lon.size):
+                    fout0.write('%6d %8.3f %8.3f %6.3f %6.3f 0.0\n'%(jj,lon.flatten()[jj],lat.flatten()[jj],tvel[jj],tvel1[jj]))
+                fout0.close()
             velall[:,ii]  = tvel
             ii+=1
         except Exception as err:
             print('cannot load %s because %s'%(ifile,err))
 
+    print('total # of models %d'%ii)
     # remove empty traces
     velall = velall[:,:ii]
 
     # averaged model and get stad
-    tvel = np.mean(velall,axis=1)
-    vel_std = np.std(velall,axis=1)
-
-    # calculate the absolute vel and vel pertb
-    vel_abs = -tvel*ave*ave+ave
+    tvel = -velall*ave*ave+ave
+    vel_abs = np.mean(tvel,axis=1)
+    vel_std = np.std(tvel,axis=1)
     vel_per = (vel_abs-ave)/ave*100
-    #print('old and new averages are %6.3f %6.3f'%(ave,nave))
 
-    # output into txt file
-    latgrid = np.linspace(latmin+nlat*dlat,latmin,nlat+1) 
-    longrid = np.linspace(lonmin,lonmin+nlon*dlon,nlon+1)
-    lon,lat = np.meshgrid(longrid,latgrid)
+    #######################################
+    ######## CALCULATE RAYPATPH ###########
+    #######################################
 
+    if cal_raypath:
+        sphgrid = subf.construct_grid(geogrid,ave)                                  # spherical grid data and initial velocity model
+        xygrid = np.zeros(shape=(tlon.size,2),dtype=np.float32)
+        ttheta,tphi = subf.geo2sph(tlat.flatten(),tlon.flatten())
+        xygrid[:,0],xygrid[:,1] = subf.sph2xyz(ttheta,tphi)
+        sphgrid['vel1'] = np.reshape(vel_abs,(1,sphgrid['ntheta'],sphgrid['nphi']))
+
+        #######################################
+        ######## CALCULATE RESIDUALS ##########
+        #######################################
+
+        if cal_resid:
+            fname = os.path.join(rootpath,'residuals/tomo_'+str(per)+'s.txt')                     
+            res_initial,res_final = subf.calculate_residuals(data,geogrid,sphgrid,fname)
+
+            # plot the histogram
+            nbin = 30
+            plt.close('all')
+            plt.hist(res_initial,nbin,density=True,facecolor='g',alpha=0.6,label='initial')
+            plt.hist(res_final,nbin,density=True,facecolor='b',alpha=0.6,label='final')
+            plt.legend(loc='upper right')
+            plt.savefig(rootpath+'/residuals/hist'+str(per)+'s.pdf',format='pdf',mpi=300)
+
+        srcs    = data.evt_id.unique()
+        # loop through each virtual source
+        for isc in srcs:
+            print(len(srcs),isc)
+
+            # source location to initialize the traveltime field
+            arr4rc = data[data.evt_id.values==isc]
+            lat_s,lon_s   = arr4rc['evt_lat'].iloc[0],arr4rc['evt_lon'].iloc[0]
+
+            # generate parameters
+            solver = subf.solve_source_region(geogrid,sphgrid,lat_s,lon_s,1)
+
+            # loop through each receiver 
+            for irc in range(len(arr4rc)):
+                try:
+
+                    # receiver location
+                    lat_r,lon_r = arr4rc['sta_lat'].iloc[irc],arr4rc['sta_lon'].iloc[irc]
+                    theta_r,phi_r = subf.geo2sph(lat_r,lon_r)
+                    src_pos = (RR,theta_r,phi_r)
+
+                    # do ray tracing here
+                    ray = solver.trace_ray(np.array(src_pos))
+                    ray = pd.DataFrame(ray,columns=('rho','theta','phi'))
+                    ray['x'] = ray['rho']*np.sin(ray['theta'])*np.cos(ray['phi'])
+                    ray['y'] = ray['rho']*np.sin(ray['theta'])*np.sin(ray['phi'])
+                    dist = scipy.spatial.distance.cdist(ray[['x','y']].values,xygrid)
+                    argmin = np.argmin(dist,axis=1)
+                    uargmin = np.unique(argmin)
+                    tcount[uargmin] += 1
+
+                except Exception as err:
+                    print('ERROR for',irc,isc,err)               
+
+    fout1 = open('{0:s}/ray_{1:02d}s.txt'.format(rootpath,per),'w')
+    ttlon,ttlat = tlon.flatten(),tlat.flatten()
+    for ii in range(tlon.size):
+       fout1.write('%8.3f %8.3f %d\n'%(ttlon[ii],ttlat[ii],tcount[ii]))
+    fout1.close()
     # format the data column
-    fout = open(rootpath+'/models/tomo_'+str(per)+'s.txt','w')
+    fout = open('{0:s}/models/L_tomo_{1:02d}s.txt'.format(rootpath,per),'w')
     for ii in range(lon.size):
         fout.write('%6d %8.3f %8.3f %6.3f %6.3f %8.4f\n'%(ii,lon.flatten()[ii],lat.flatten()[ii],vel_per[ii],\
             vel_abs[ii],vel_std[ii]))
     fout.close()
 
-
-    #######################################
-    ######## CALCULATE RESIDUALS ##########
-    #######################################
-
-    if cal_resid:
-        sphgrid = subf.construct_grid(geogrid,ave)  
-        sphgrid['vel1'] = np.reshape(vel_abs,(1,sphgrid['ntheta'],sphgrid['nphi']))
-        fname = os.path.join(rootpath,'residuals/tomo_'+str(per)+'s.txt')                     
-        res_initial,res_final = subf.calculate_residuals(data,geogrid,sphgrid,fname)
-
-        # plot the histogram
-        nbin = 30
-        plt.close('all')
-        plt.hist(res_initial,nbin,density=True,facecolor='g',alpha=0.6,label='initial')
-        plt.hist(res_final,nbin,density=True,facecolor='b',alpha=0.6,label='final')
-        plt.legend(loc='upper right')
-        plt.savefig(rootpath+'/residuals/hist'+str(per)+'s.pdf',format='pdf',mpi=300)
-
     # plot the model
     if plot_final:
-        subf.plot_tomo(geogrid,vel_abs,rootpath,-1)
+        subf.plot_tomo(geogrid,vel_abs,rootpath,-1)

src/make_syn_2D.py

@@ -18,21 +18,22 @@
 syn_option = 'checkerboard'
 
 # absoluate path for traveltime data and for output
-rootpath = '/Users/chengxinjiang/Documents/Github/VoroTomo_SW/example/synthetic'
-dfiles   = glob.glob(os.path.join(rootpath,'dispersion/LVC_snr8_2wl_8s.dat'))
+rootpath = '/Users/chengxin/Documents/ANU/St_Helens/Voro_Tomo/synthetic/Oct4_WithNoise_0.5s_spacing2/Love'
+dfiles   = glob.glob(os.path.join(rootpath,'selection_*s.dat'))
 
 # full 2D research area
-latmin,dlat,nlat = 35.5,0.05,110                                            # latitude range of the target region
-lonmin,dlon,nlon = -122.5,0.05,100                                          # longitude range of the target region
+latmin,dlat,nlat = 44.5,0.02,190                                            # latitude range of the target region
+lonmin,dlon,nlon = -124.8,0.02,310                                          # longitude range of the target region
 geogrid = {'latmin':latmin,'lonmin':lonmin,                                 # assemble geographic info into a dict 
             'dlat':dlat,'dlon':dlon,
             'nlat':nlat,'nlon':nlon}
 
 if syn_option=='checkerboard':
     # target anomaly information
     anomaly = {'avs':3,'amp':0.3,                                           # amplitude in percentage
-                'size':15,'spacing': True,                                  # size of anomaly and whether of spacing
-                'noise':0.3}
+                'size':15,'spacing': 2,                                      # size of anomaly and whether of spacing
+                'noise_ave':0,
+                'noise_std':0.5}
 
     # make sphgrid for traveltime prediction
     sphgrid = subf.construct_grid(geogrid,anomaly['avs'])                       # spherical grid data and initial velocity model
@@ -46,7 +47,7 @@
     anomaly = {'avs':3,'amp':-0.3,                                               # amplitude in percentage
                 'latmin':46.1,'latmax':46.4,                                     # size of anomaly and whether of spacing
                 'lonmin':-122.4,'lonmax':-122.1,
-                'noise':0.3}
+                'noise':0}
 
     # make sphgrid for traveltime prediction
     sphgrid = subf.construct_grid(geogrid,anomaly['avs'])                       # spherical grid data and initial velocity model
@@ -97,8 +98,9 @@
                 theta_r,phi_r = subf.geo2sph(lat_r,lon_r)
                 src_pos  = (RR,theta_r,phi_r)
                 syn_data = solver.traveltime.value(np.array(src_pos))
-                if anomaly['noise']:
-                    syn_data += np.random.normal(loc=anomaly['noise'],scale=0.05)
+                if anomaly['noise_std']:
+                    syn_data += np.random.normal(loc=anomaly['noise_ave'],scale=anomaly['noise_std'])
+
                 fout.write('%6d,%6d,%8.3f,%8.3f,%6d,%8.3f,%8.3f,%8.3f,%8.3f\n'%(arr4rc['index'].iloc[irc],\
                     isc,lon_s,lat_s,arr4rc['sta_id'].iloc[irc],lon_r,lat_r,arr4rc['dist'].iloc[irc]/syn_data,\
                     arr4rc['dist'].iloc[irc]))
@@ -107,4 +109,4 @@
                 print('ERROR for',irc,isc,err)
 
     fout.close()
-    print('synthetic done for %s'%dfiles[iper])
+    print('synthetic done for %s'%dfiles[iper])

src/remove_data_residuals.py

@@ -16,8 +16,8 @@
 '''
 
 # define rootpath and find all files
-rootpath = '/Users/chengxinjiang/Documents/ANU/Voro_Tomo/real_data'
-pper = [3,6,8,12,14,18,20,24,28]
+rootpath = '/Users/chengxin/Documents/ANU/St_Helens/Voro_Tomo/real_data_inversion_Cascadia/Love_group'
+pper = [2,3,4,5,6,8,10,12,14,16,18,20,24,28,32,36,40]
 
 # loop through each period
 for ii in range(len(pper)):
@@ -34,8 +34,8 @@
     dist = np.array(data['dist'])
     vel  = np.array(data['vel'])
     tt   = np.array(dist/vel)
-    res1 = np.array(data['res_initial'])
-    res2 = np.array(data['res_final'])
+    res1 = np.array(data['res_initial']).astype(np.float32)
+    res2 = np.array(data['res_final']).astype(np.float32)
 
     # get statistics
     ave1 = np.mean(res1)
@@ -45,7 +45,7 @@
 
     # make selections
     indx = np.where((res2>=ave2-2*std2)&(res2<=ave2+2*std2))[0]
-    fout = open(rootpath+'/dispersion/selection2_'+str(per)+'s.dat','w')
+    fout = open(rootpath+'/dispersion/selection_'+str(per)+'s.dat','w')
     fout.write('index,evt_id,evt_lon,evt_lat,sta_id,sta_lon,sta_lat,vel,dist\n')
     for jj in range(len(indx)):
         tindx = indx[jj]
@@ -81,6 +81,6 @@
     plt.text(15+np.min(dist),0.85*np.max(tt),text)
     plt.title(str(per)+'s final model')
     plt.colorbar()
-    outname = rootpath+'/dispersion/obs2_residual_'+str(per)+'s.pdf'
+    outname = rootpath+'/dispersion/obs_residual_'+str(per)+'s.pdf'
     plt.savefig(outname,format='pdf',dpi=300)
     #plt.show()