Body Waves Synthetics Cross-correlation Functions¶

## Plot archive
plt.figure()
plt.pcolormesh(time, distance, archive,cmap='seismic', vmin= -10, vmax = 10)
plt.colorbar()
plt.xlabel('time')
plt.ylabel('distance')
plt.xlim(0, 3600)
plt.show()

## Taper
file_name_taper = '../../data/NOISE_vertforce_dirac_0-ak135f_1.s_3600_tapered_P_PP.h5'

if os.path.isfile(file_name_taper):
    print('Tapered archive already exists')
    ## Open tapered archive
    h5_file = h5py.File(file_name_taper, 'r')
    tapered_archive = h5_file['WAVEFORMS'][:]
    
else:
    tapered_archive = taper_axisem_archive_body_waves(time, distance, archive_name=path_file_axisem, phase1='P', phase2='PP')

## Plot Tapered archive 

plt.figure(figsize=(10, 5))
plt.pcolormesh(time, distance, tapered_archive, cmap='seismic', vmin=-10, vmax = 10)
plt.colorbar()
plt.xlabel('time')
plt.ylabel('distance')
plt.xlim(0, 3600)
plt.savefig('GF_archive.png', dpi=300, bbox_inches='tight')

file_name_spectrum = '../../data/spectrum_archive_tapered_P_PP.h5'

if os.path.isfile(file_name_spectrum):
    print('Spectrum archive already exists')
    ## Open Spectrum archive
    h5_file = h5py.File(file_name_spectrum, 'r')
    spectrum_archive = h5_file['SPECTRUM'][:]
    frequency = h5_file['frequency'][:]
    distance = h5_file['distance'][:]
    print(spectrum_archive.shape)
else:
    ## Create Spectrum archive from the tapered one
    create_spectrum_archive(time, distance, tapered_archive, h5_name_spectrum = file_name_spectrum)

## Station Pairs Metadata
df = pd.read_csv('../../data/stations_pair_bodywaves.txt', sep=',', header=None, index_col=False, names=['idA', 'idB', 'latA', 'latB', 'lonA', 'lonB'], dtype={'idA':str, 'idB':str, 'latA':float, 'latB':float, 'lonA':float, 'lonB':float}, engine='python')
# ## Plot
fig = plt.figure(figsize=(10,10))
ax = fig.add_subplot(111, projection=ccrs.Mollweide())
#ax.set_extent(extent, crs=ccrs.PlateCarree())
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.BORDERS, linestyle='-', alpha=1)
ax.add_feature(cfeature.LAKES, alpha=0.5)
ax.add_feature(cfeature.RIVERS)
ax.add_feature(cfeature.STATES)
ax.add_feature(cfeature.LAND, color='linen', zorder=0)
ax.add_feature(cfeature.OCEAN, color='lightsteelblue', zorder=0)
for i in range(len(df)):
    ax.annotate(df.idA[i], (df.lonA[i]-0.1, df.latA[i]+0.8), transform = ccrs.PlateCarree(), family = 'sans-serif')
    ax.annotate(df.idB[i], (df.lonB[i]+0.1, df.latB[i]-0.8), transform = ccrs.PlateCarree(), family = 'sans-serif')
    ax.plot((df.lonA[i], df.lonB[i]),(df.latA[i], df.latB[i]), transform = ccrs.PlateCarree(), color='k', alpha=0.5, linewidth=3)
df.plot.scatter(x='lonA', y='latA', s=30, ax=ax, c='r', transform = ccrs.PlateCarree())
df.plot.scatter(x='lonB', y='latB', s=30, ax=ax, c='r', transform = ccrs.PlateCarree())
plt.show()

## Compute CCF and save it
### Station Pair Info
id_A = df.idA[0]
id_B = df.idB[0]
coords_staA = (df.lonA[0], df.latA[0])
coords_staB = (df.lonB[0], df.latB[0])

#### Keep Correlations
corr_3h = np.zeros((len(dates_vector), N_time))

### Date Info
for i, date_vect in enumerate(dates_vector):
    ### CCF
    filename_ccf = './CCF/ccf_%s_%s_%d%02d%02d%02d.nc'%(id_A, id_B, date_vect[0],date_vect[1],date_vect[2],date_vect[3]) # File Name Understandable

    if os.path.exists(filename_ccf):  # Check if file already exists
        print('%s already computed'%filename_ccf)
        corr = xr.open_dataarray(filename_ccf)
    else:  # if not, compute CCF
        t1 = ttime.perf_counter()
        corr, time_corr = ccf_computation(coords_staA, coords_staB, path_model, date_vect, spectrum_synth,extent=extent, fe=fe, N=N_time, comp='Z')
        t2 = ttime.perf_counter()
        print(" Time for station pair: ", id_A, id_B, " elapsed: ", t2-t1)
        try:
            corr.to_netcdf('./CCF/ccf_%s_%s_%d%02d%02d%02d.nc'%(id_A, id_B, date_vect[0],date_vect[1],date_vect[2],date_vect[3]))
        except:
            print("make directory")
            os.mkdir('./CCF')
            corr.to_netcdf('./CCF/ccf_%s_%s_%d%02d%02d%02d.nc'%(id_A, id_B, date_vect[0],date_vect[1],date_vect[2],date_vect[3]))
    corr_3h[i, :] = corr.values
    # TAKES 12 MIN #

## Plot stack 

plt.figure()
plt.plot(corr.time, np.mean(corr_3h, axis=0), 'k')
plt.xlim(400, 600)
plt.xlabel('Lag Time (s)')
plt.title('%s - %s stack'%(id_A, id_B))
plt.show()

## Plot CCF
normalize_corr = np.nanmax(abs(corr_3h))
dates_corr = dates

fig, ax = plt.subplots(figsize=(10,8))
im = plt.pcolormesh(corr.time, np.arange(len(dates_vector)), corr_3h/normalize_corr, cmap='seismic', vmin= -0.8, vmax = 0.8)
#ax.plot(corr.time, corr_3h.T/normalize_corr + np.arange(len(dates_vector)), 'k') 
plt.colorbar(im)
ax.set_xlabel('Lag Time (s)')
ax.set_ylabel('Date')
ax.set_title('%s-%s CCF'%(id_A, id_B))
ax.set_xlim(400, 600)
ax.set_ylim(0, len(dates_vector))
ax.set_yticks(np.arange(len(dates_vector))[::8])
ax.set_yticklabels(dates[::8].strftime('%Y-%m-%dT%H'))

plt.savefig('ccf_%s_%s_synth.png'%(id_A, id_B), dpi=300, bbox_inches='tight')