Source code for create_2D_migrated_profile_DK2

"""
Location: Chavannes-pres-renens, CH
Date: April 2022
Author: Konstantinos Michailos
"""
import rfmpy.utils.migration_mtz as mtz
import numpy as np
import os
from obspy.geodetics import degrees2kilometers, kilometers2degrees
import rfmpy.utils.migration_plots_spher as plot_migration_sphr

# Define paths


[docs]
work_dir = os.getcwd()





[docs]
path2grid = '/home/' + work_dir.split('/')[2] + '/Desktop/mtz_example/'



# Read the 3D grid (epcrust.npy) of stacked migrated RF amplitudes.
with open(path2grid + 'example.npy', 'rb') as f:
    mObs = np.load(f)

## Define MIGRATION parameters
# Ray-tracing parameters


[docs]
inc = 2  # km




[docs]
zmax = 800 # km


# Determine study area (x -> perpendicular to the profile)


[docs]
minx = -13.0 # degrees 2.optional:2 or -4




[docs]
maxx = 46.0 # degrees 2.optional:30 or 38




[docs]
pasx = 0.26 # degrees oldest 0.38




[docs]
miny = 30.0 # degrees 2.optional:41 or 38




[docs]
maxy = 64.0 # degrees 2.optional:51 or 54




[docs]
pasy = 0.18 # degrees oldest 0.27




[docs]
minz = -5 # km


# maxz needs to be >= zmax


[docs]
maxz = 800 # km




[docs]
pasz = 2 # km


# Pass all the migration parameters in a dictionary to use them in functions called below


[docs]
m_params = {'minx': minx, 'maxx': maxx,
            'pasx': pasx, 'pasy': pasy, 'miny': miny, 'maxy': maxy,
            'minz': minz, 'maxz': maxz, 'pasz': pasz, 'inc': inc, 'zmax': zmax}



# Define path to RFs


[docs]
path = '/home/' + work_dir.split('/')[2] + '/Desktop/mtz_example/RF/'



# Read station coordinates from the rfs (sac files) in a pandas dataframe


[docs]
sta = mtz.read_stations_from_sac(path2rfs=path)





[docs]
profile_A = np.array([[8, 50.5], [30, 45.2]])




[docs]
prof_name = 'Cross-section_A_and_A'



G2_, sta_, xx, zz = plot_migration_sphr.create_2d_profile(mObs, m_params, profile_A, sta, swath=300, plot=True)
mObs = mtz.ccp_smooth(G2_, m_params)


[docs]
mObs = mtz.ccpFilter(mObs)



# File for creating cross-sections with GMT
for i, x in enumerate(xx):
    for j, z in enumerate(zz):
        print(kilometers2degrees(x), z, mObs[i,j])
        with open(path2grid + prof_name + '.txt', 'a') as of:
            of.write('{} {} {} \n'.
                     format(kilometers2degrees(x), z, mObs[i, j]))