plot_tool.py

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib.gridspec import GridSpec
from obspy import UTCDateTime
import matplotlib as mpl
from pyproj import Geod
import seaborn as sns
import os
import glob
import matplotlib.dates as mdates
import datetime
import datetime
import os


def calculate_true_distance(dist_horizontal_m, source_depth_km, inversion_dep_km):
    source_depth_m = source_depth_km *1000
    inversion_dep_m = inversion_dep_km*1000
    dist_vertical_m = source_depth_m - inversion_dep_m 
    dist_true_m = np.sqrt(dist_horizontal_m**2 + dist_vertical_m**2)
    return dist_true_m

def cal_misfit(src_lon,src_lat,src_dep, loc_lon ,loc_lat ,loc_dep):
    geod = Geod(ellps='WGS84')
    #震中距
    dist_horizontal_loc_m = geod.inv(src_lon, src_lat,loc_lon, loc_lat)[2]
    # print(src_lon,src_lat,src_dep, loc_lon ,loc_lat ,loc_dep , dist_horizontal_loc_m)
    #考虑高程差
    return calculate_true_distance(dist_horizontal_loc_m,src_dep,loc_dep)/1000

def split_file(file_path, split_size,name):
    with open(file_path, 'r') as f:
        line_count = 0
        file_number = 1
        out_dir = "/media/wzm/HLPBook/DiTingProject/data/synthesis_data/azimyths/real/split/20250521/"
        if not os.path.exists(f'{out_dir}/split{file_number}'):
            os.makedirs(f'{out_dir}/split{file_number}')
        current_file = open(f'{out_dir}/split{file_number}/{name}', 'w')
        for line in f:
            line = line[:-3]+line[-1]
            current_file.write(line)
            line_count += 1
            
            if line_count == split_size:
                current_file.close()
                file_number += 1
                if not os.path.exists(f'{out_dir}/split{file_number}'):
                    os.makedirs(f'{out_dir}/split{file_number}')
                    
                current_file =open(f'{out_dir}/split{file_number}/{name}', 'w')
                line_count = 0
        current_file.close()


def merge_real_output(combined_file , real_dir ):
    # 合并real输出
    # combined_file = open(dir+'/'+model+'/'+test+'.csv', "w")
    # 逐个将文本文件追加到新文件中
    for catalog in glob.glob(real_dir+"/*.catalog_sel.txt"):
        with open( catalog, "r") as file:
            content = file.read()
            combined_file.write(content)
    # 关闭新文件
    combined_file.close()
    print("合并完成！")


def get_loc_df(model,test,dir= "/media/wzm/HLPBook/DiTingProject/results/synthesis_data/"):

    if model == 'REAL' :

        REAL_df = pd.read_csv(dir+'/'+model+'/'+test+'.csv',  delim_whitespace=True,header=None,usecols=[1,2,3,4,7, 8, 9],names=['year','month','day','date','loc_lat', 'loc_lon', 'loc_dep'])
        # REAL_df = pd.read_csv('/media/wzm/HLPBook/DiTingProject/results/synthesis_data/archive/azimyths/real_1min_sep/merge_catalog.csv',  
        #                      delim_whitespace=True,header=None,usecols=[1,2,3,4,7, 8, 9],names=['year','month','day','date','loc_lat', 'loc_lon', 'loc_dep'])

        REAL_df['time'] = pd.to_datetime(REAL_df['year'].map(str) +'-'+ REAL_df['month'].map(str)+'-'+ (REAL_df['day'].map(str)) +' '+ REAL_df['date'] ,utc =True)
        return REAL_df


    # # pyocto 结果
    if model == 'pyocto' :
        pyocto_df = pd.read_csv(dir+'/'+model+'/'+test+'/associate_cat.csv',skiprows=1, sep=',',header=None, usecols=[1, 6, 7, 8], names=['time','loc_lat', 'loc_lon', 'loc_dep'])
        # pyocto_df = pd.read_csv('/media/wzm/HLPBook/DiTingProject/results/synthesis_data/archive/azimyths/real_1min_sep/pyocto_result1.csv', 
        # sep=',',skiprows=1, header=None, usecols=[1, 6, 7, 8], names=['time','loc_lat', 'loc_lon', 'loc_dep'])
        pyocto_df['time'] = pd.to_datetime(pyocto_df['time'],utc =True)
        return pyocto_df

    # #gamma
    if model == 'gamma':
        gamma_df = pd.read_csv(dir+'/'+model+'/'+test+'/gamma_events_no_duplicate.csv')
        gamma_df.rename(columns={"longitude": "loc_lon", "latitude": "loc_lat", "depth_km": "loc_dep"}, inplace=True)
        gamma_df['time'] = pd.to_datetime(gamma_df['time'],utc =True)
        return gamma_df


def plot_grid_misfit(loc_df , model, test  ,title_suffixes = "",station_csv = '/media/wzm/HLPBook/DiTingProject/results/synthesis_data/single_source/station.csv'  ,reference_time= UTCDateTime("2025-05-21T00:00:00.000")):
    # 读取通道信息
    df = pd.read_csv(station_csv, sep='\s+')
    df_downsampled = df.iloc[::4, :].reset_index(drop=True)  
    das_lon = df_downsampled['longitude'].values
    das_lat = df_downsampled['latitude'].values
    # 创建输入震源网格
    elat_trial_arr = np.linspace(23.35, 24.35, 11)
    elon_trial_arr = np.linspace(114.0, 115.0, 11)
    # elat_trial_arr = np.linspace(23.75, 24.00, 11)
    # elon_trial_arr = np.linspace(114.40, 114.60, 11)

    elat_grid, elon_grid = np.meshgrid(elat_trial_arr, elon_trial_arr)  #经纬度网格
    dep_grid = np.linspace(2, 20, 10)  # 深度网格
    # 展平网格并创建输入震源数组
    src_lon = np.repeat(elon_grid.flatten(), len(dep_grid))
    src_lat = np.repeat(elat_grid.flatten(), len(dep_grid))
    src_dep = np.tile(dep_grid, len(elat_grid.flatten()))

    #时间
    # reference_time = UTCDateTime("2025-05-21T00:00:00.000")

    # 展平网格并创建震源DataFrame
    sources = []
    for i, (lat, lon) in enumerate(zip(elat_grid.flatten(), elon_grid.flatten())):
        for j, dep in enumerate(dep_grid): 
            sources.append({
                #'source_id': f"src_{i:04d}_dep{dep:.1f}km",
                'src_lon': lon.round(2),
                'src_lat': lat.round(2),
                'src_dep_km': dep,
                'time': (reference_time + 60 *60* (i * len(dep_grid) + j)).datetime
            })
    source_df = pd.DataFrame(sources)
    pyocto_df = loc_df
    #根据发震时间 目录对齐
    source_df['time'] = pd.to_datetime(source_df['time'],utc =True)
    pyocto_df['time'] = pd.to_datetime(pyocto_df['time'],utc =True)
    pyocto_df = pyocto_df.sort_values(by=['time'])
    misfit_df = pd.merge_asof( source_df , pyocto_df,on="time", tolerance=pd.Timedelta("30min") , direction='nearest')
    # print(misfit_df.isnull().sum())
    #计算misfit
    misfit_df["misfit"]=misfit_df.apply(lambda x:cal_misfit(x['src_lon'],x['src_lat'] ,x['src_dep_km'],x['loc_lon'], x['loc_lat'], x['loc_dep']),axis=1)
    # print(misfit_df.isnull().sum())
    misfit_df["misfit"].fillna(1000, inplace=True)

    # 设置colorbar
    cmap = mpl.cm.viridis_r
    bounds = [0.5, 1, 2, 4, 8, 12, 16]
    # bounds = [0.5, 1, 1.5, 2, 2.5, 3, 4]
    norm = mpl.colors.BoundaryNorm(bounds, cmap.N, extend='both')

    #用seaborn画图是直接给dataframe, 是最准确的，但是画不出阵列分布
    plt.figure(figsize=(10, 10))
    plt.suptitle(model+title_suffixes, y=1 ,  fontsize='x-large')
    for i ,dep in enumerate(dep_grid[1:9:2]):
        plot_misfit_df = misfit_df[misfit_df['src_dep_km'] == dep].pivot(index='src_lat', columns='src_lon', values='misfit').sort_index(axis=0, ascending=False)
        ax = plt.subplot(220+i+1)
        plt.title("Misfit Depth = " +str(dep) )
        ax.set_aspect(1)
        ax = sns.heatmap(data=plot_misfit_df,square=True,norm=norm, cmap='viridis_r' , shading='auto',alpha=0.7 , annot=True  ,zorder=-10) 
        sns.scatterplot(df_downsampled, x='longitude',y='latitude', marker='^', color='red', 
                label='DAS Channels', alpha=0.7 ,  ax=ax , zorder=4)
        # plt.scatter(das_lon, das_lat, marker='.', color='black', 
        #             label='DAS Channels', s=10, alpha=0.7,zorder=4 )
        ax.set_rasterization_zorder(0)
        ax.get_legend().remove()

    plt.tight_layout()
    plt.savefig('/media/wzm/HLPBook/DiTingProject/results/synthesis_data/'+model+'/'+test+"/misfit_with_num"+title_suffixes+".png")
    #用imshow可以画出阵列分布，但是打上数字会没那么好看
    plt.figure(figsize=(10, 10))
    plt.suptitle(model+title_suffixes, y=1 ,  fontsize='x-large')
    for i ,dep in enumerate(dep_grid[1:9:2]):
        plot_misfit_df = misfit_df[misfit_df['src_dep_km'] == dep].pivot(index='src_lat', columns='src_lon', values='misfit').sort_index(axis=0, ascending=False)
        # print(plot_misfit_df)
        ax = plt.subplot(220+i+1)
        plt.title("Misfit Depth = " +str(dep) )
        plt.imshow(plot_misfit_df, norm=norm, cmap='viridis_r' , alpha=0.7 , extent =[113.95, 115.05,23.30, 24.4])  # 大网格
        # plt.imshow(plot_misfit_df, norm=norm, cmap='viridis_r' , alpha=0.7 , extent =[114.39, 114.61,23.74, 24.01])  #小网格
        plt.colorbar()
        ax.set_aspect(1)
        sns.scatterplot(df_downsampled, x='longitude',y='latitude', marker='.', color='red', 
                label='DAS Channels', alpha=0.7 ,  ax=ax , zorder=4)
        ax.set_rasterization_zorder(0)
        # for index, row in misfit_df[misfit_df['src_dep_km'] == dep].iterrows():#先lat+，再lon+
        #     src_lon, src_lat , misfit  = row['src_lon'], row['src_lat'], row['misfit']
        #     plt.text(src_lon, src_lat , round(misfit,1) ,ha="center", va="center", color="black" , zorder=7)
        # 设置图例位置
        ax.get_legend().remove()
    plt.tight_layout()
    plt.savefig('/media/wzm/HLPBook/DiTingProject/results/synthesis_data/'+model+'/'+test+"/misfit"+title_suffixes+".png")
    return misfit_df

def plot_polar_misfit( loc_df , model, test , title_suffixes = "",station_csv = '/media/wzm/HLPBook/DiTingProject/results/synthesis_data/single_source/station.csv'   ,reference_time= UTCDateTime("2025-05-21T00:00:00.000")):
    # 初始化地理计算工具
    geod = Geod(ellps='WGS84')

    # 读取通道信息
    df = pd.read_csv(station_csv, sep='\s+')
    df_downsampled = df.iloc[::4, :].reset_index(drop=True)  
    das_lon = df_downsampled['longitude'].values
    das_lat = df_downsampled['latitude'].values
    center_lon = np.mean(das_lon)
    center_lat = np.mean(das_lat)
    das_distance = [ geod.inv(center_lon, center_lat, das_lon[i], das_lat[i])[2]/1000 for i in range(len(das_lat))]
    forward_azimuth= [ geod.inv(center_lon, center_lat, das_lon[i], das_lat[i])[0] for i in range(len(das_lat))]
    back_azimuth = [ geod.inv(center_lon, center_lat, das_lon[i], das_lat[i])[1] for i in range(len(das_lat))]
    # 设置极坐标网格参数
    distances = np.arange(2, 30, 2)  # 1-30km，间隔2km
    azimuths = np.linspace(0, 360, 16, endpoint=False)  # 8个方位角（0-360°）
    depths = np.linspace(2, 20, 10)  # 深度网格 2-20km
    print(distances)

    # 创建震源
    sources = []
    # reference_time = UTCDateTime("2025-05-21T00:00:00.000")

    for i, dist in enumerate(distances):
        for j, az in enumerate(azimuths):
            # 计算每个距离和方位角对应的经纬度
            lon, lat, _ = geod.fwd(center_lon, center_lat, az, dist * 1000)  # 距离转换为米
            
            for k, dep in enumerate(depths):
                sources.append({
                    'src_lon': lon,
                    'src_lat': lat,
                    'src_dep_km': dep,
                    'time': (reference_time + 60*60* (i * len(azimuths) * len(depths) 
                                                + j * len(depths) 
                                                + k)).datetime,
                    'distance_km': dist,
                    'azimuth': az
                })
    source_df = pd.DataFrame(sources)
    # source_df.rename(columns={"src_lon":"longitude" , "src_lat": "latitude","src_dep_km":"depth_km"}, inplace=True)
    # source_df.to_csv('/media/wzm/HLPBook/DiTingProject/data/synthesis_data/'+test+'_standard_catalog.csv', index=True)
    pyocto_df = loc_df
    source_df['time'] = pd.to_datetime(source_df['time'],utc =True)
    # print(source_df.head(5))
    #根据发震时间 目录对齐
    pyocto_df = pyocto_df.sort_values(by=['time'])
    misfit_df = pd.merge_asof( source_df , pyocto_df,on="time", tolerance=pd.Timedelta("70min") , direction='nearest')
    # print(misfit_df.isnull().sum())
    #计算misfit
    misfit_df["misfit"]=misfit_df.apply(lambda x:cal_misfit(x['src_lon'],x['src_lat'] ,x['src_dep_km'],x['loc_lon'], x['loc_lat'], x['loc_dep']),axis=1)
    print(misfit_df.isnull().sum())
    misfit_df["misfit"].fillna(100, inplace=True)
    #设置colorbar
    cmap = mpl.cm.viridis_r
    bounds = [0.5, 1, 2, 4, 8, 12, 16]
    norm = mpl.colors.BoundaryNorm(bounds, cmap.N, extend='both')
    #plot 
    fig = plt.figure(figsize=(20, 5))
    plt.suptitle(model+title_suffixes, y=1 ,  fontsize='x-large')
    for i ,dep in enumerate(depths[1:9:2]):
        plot_misfit_df = misfit_df[misfit_df['src_dep_km'] == dep].pivot(index='distance_km', columns='azimuth', values='misfit')
        p_azimuths= plot_misfit_df.columns.tolist()
        p_distances= plot_misfit_df.index.tolist()
        theta = np.deg2rad(p_azimuths)
        ax = plt.subplot(140+i+1 , polar = True)
        plt.title("Misfit Depth = " +str(dep) )
        # ax = sns.heatmap(data=plot_misfit_df,square=False,norm=norm, cmap='viridis_r' , shading='auto',alpha=0.7 , annot=True  ,zorder=-10) 
        a = plt.pcolormesh(theta, p_distances, plot_misfit_df.values.tolist(),norm=norm, cmap='viridis_r' , shading='auto',alpha=0.8)
        plt.grid(False)
        ax.set_theta_zero_location('N')
        ax.set_theta_direction(-1)
        ax.set_rlabel_position(90)
        plt.scatter(np.deg2rad(forward_azimuth), das_distance, marker='.', color='red', 
                    label='DAS Channels', s=10, alpha=0.7)
        # plt.colorbar()
        for index, row in misfit_df[misfit_df['src_dep_km'] == dep].iterrows():#先lat+，再lon+
            disi, azi , misfit  = row['distance_km'], row['azimuth'], row['misfit']
            if misfit == 100 :
                plt.scatter(np.deg2rad(azi), disi , marker='^', color='gray', 
                    label='no loc', s=10, alpha=0.7)
                # plt.text(np.deg2rad(azi), disi , round(misfit,1) ,ha="center", va="center", color="black" , zorder=7)
    fig.subplots_adjust(right=0.9)
    cbar_ax = fig.add_axes([0.92, 0.12, 0.015, 0.76])
    cb = plt.colorbar(a, cax=cbar_ax ,label='Average Misfit Value(km)', pad=0.1)
    # plt.tight_layout()

    plt.savefig('/media/wzm/HLPBook/DiTingProject/results/synthesis_data/'+model+'/'+test+"/misfit"+title_suffixes+".png")
    return misfit_df

def compare_catalog(loc_df ,model, test , title_suffixes=""):
    catalog_label="Standard"
    pyocto_events = loc_df
    figure_dir = lambda x: "/media/wzm/HLPBook/DiTingProject/results/synthesis_data/"+model+f"/{test}/"+x
    model_label = model+title_suffixes
    if os.path.exists("/media/wzm/HLPBook/DiTingProject/data/synthesis_data/"+test+"_standard_catalog.csv"):
        standard_catalog = pd.read_csv("/media/wzm/HLPBook/DiTingProject/data/synthesis_data/"+test+"_standard_catalog.csv")
        starttime = standard_catalog["time"].min()
        endtime = standard_catalog["time"].max()
    else:
        standard_catalog = None
        starttime = pyocto_events["time"].min()
        endtime = pyocto_events["time"].max()
    print(standard_catalog["longitude"].max() , standard_catalog["longitude"].min())
    print(standard_catalog["latitude"].max() , standard_catalog["latitude"].min())

    plt.figure()
    plt.hist(pyocto_events["time"], range=(starttime, endtime), bins=24, edgecolor="k", alpha=1.0, linewidth=0.5, label=f"{model_label}: {len(pyocto_events['time'])}")
    if standard_catalog is not None:
        plt.hist(standard_catalog["time"], range=(starttime, endtime), bins=24, edgecolor="k", alpha=0.6, linewidth=0.5, label=f"{catalog_label}: {len(standard_catalog['time'])}")
    plt.ylabel("Frequency")
    plt.xlabel("Date")
    plt.gca().autoscale(enable=True, axis='x', tight=True)
    plt.gca().xaxis.set_major_formatter(mdates.DateFormatter('%m-%d:%H'))
    plt.gcf().autofmt_xdate()
    plt.legend()
    plt.savefig(figure_dir("earthquake_number"+title_suffixes+".png"), bbox_inches="tight", dpi=300)
    plt.savefig(figure_dir("earthquake_number"+title_suffixes+".pdf"), bbox_inches="tight")
    plt.show()


def plot_source_compare(loc_df,model, test , title_suffixes=""):
    pyocto_events = loc_df
    pyocto_events.rename(columns={"loc_lon":"longitude" , "loc_lat": "latitude","loc_dep":"depth"}, inplace=True)
    catalog_label="Standard"
    figure_dir = lambda x: "/media/wzm/HLPBook/DiTingProject/results/synthesis_data/"+model+f"/{test}/"+x
    model_label = model+title_suffixes

    stations = pd.read_csv("/media/wzm/HLPBook/DiTingProject/data/synthesis_data/pyocto/station.csv", sep='\s+')
    stations.rename(columns={"elevation": "elevation_m",}, inplace=True)
    standard_catalog = pd.read_csv("/media/wzm/HLPBook/DiTingProject/data/synthesis_data/"+test+"_standard_catalog.csv")
    config = {  "ylim_degree":(23.6, 24.15),
        "xlim_degree":(114.2, 114.8),
        "z(km)":(0, 22),
    }

    fig = plt.figure(figsize=plt.rcParams["figure.figsize"]*np.array([1.5,1]))
    box = dict(boxstyle='round', facecolor='white', alpha=1)
    text_loc = [0.05, 0.92]
    grd = fig.add_gridspec(ncols=2, nrows=2, width_ratios=[1.5, 1], height_ratios=[1,1])
    fig.add_subplot(grd[:, 0])
    plt.plot(pyocto_events["longitude"], pyocto_events["latitude"], '^',markersize=2, alpha=1.0)
    if standard_catalog is not None:
        plt.plot(standard_catalog["longitude"], standard_catalog["latitude"], '.', markersize=2, alpha=0.6)
    plt.axis("scaled")
    # plt.xlim(np.array(config["xlim_degree"]))
    # plt.ylim(np.array(config["ylim_degree"]))
    # plt.ylim(pyocto_events["latitude"].min()-0.2,pyocto_events["latitude"].max()+0.2)
    # plt.xlim(pyocto_events["longitude"].min()-0.2,pyocto_events["longitude"].max()+0.2)
    plt.xlabel("Latitude")
    plt.ylabel("Longitude")
    plt.gca().set_prop_cycle(None)
    plt.plot([], [], '.', markersize=10, label=f"{model_label}", rasterized=True)
    plt.plot([], [], '.', markersize=10, label=f"{catalog_label}", rasterized=True)
    plt.plot(stations["longitude"], stations["latitude"], 'k', markersize=1, alpha=0.7, label="Stations")
    plt.legend(loc="lower right")
    plt.text(text_loc[0], text_loc[1], '(i)', horizontalalignment='left', verticalalignment="top", 
            transform=plt.gca().transAxes, fontsize="large", fontweight="normal", bbox=box)

    fig.add_subplot(grd[0, 1])
    plt.plot(pyocto_events["longitude"], pyocto_events["depth"], '.', markersize=2, alpha=1.0, rasterized=True)
    if standard_catalog is not None:
        plt.plot(standard_catalog["longitude"], standard_catalog["depth_km"], '.', markersize=2, alpha=0.6, rasterized=True)
    # plt.xlim(np.array(config["xlim_degree"])+np.array([0.2,-0.27]))
    plt.ylim(config["z(km)"])
    plt.gca().invert_yaxis()
    plt.xlabel("Longitude")
    plt.ylabel("Depth (km)")
    plt.gca().set_prop_cycle(None)
    plt.plot([], [], '.', markersize=10, label=f"{model_label}")
    plt.plot([], [], '.', markersize=10, label=f"{catalog_label}")
    plt.legend(loc="lower right")
    plt.text(text_loc[0], text_loc[1], '(ii)', horizontalalignment='left', verticalalignment="top", 
            transform=plt.gca().transAxes, fontsize="large", fontweight="normal", bbox=box)

    fig.add_subplot(grd[1, 1])
    plt.plot(pyocto_events["latitude"], pyocto_events["depth"], '.', markersize=2, alpha=1.0, rasterized=True)
    if standard_catalog is not None:
        plt.plot(standard_catalog["latitude"], standard_catalog["depth_km"], '.', markersize=2, alpha=0.6, rasterized=True)
    # plt.xlim(np.array(config["ylim_degree"])+np.array([0.2,-0.27]))
    plt.ylim(config["z(km)"])
    plt.gca().invert_yaxis()
    plt.xlabel("Latitude")
    plt.ylabel("Depth (km)")
    plt.gca().set_prop_cycle(None)
    plt.plot([], [], '.', markersize=10, label=f"{model_label}")
    plt.plot([], [], '.', markersize=10, label=f"{catalog_label}")
    plt.legend(loc="lower right")
    plt.tight_layout()
    plt.text(text_loc[0], text_loc[1], '(iii)', horizontalalignment='left', verticalalignment="top", 
            transform=plt.gca().transAxes, fontsize="large", fontweight="normal", bbox=box)
    plt.savefig(figure_dir("earthquake_location"+title_suffixes+".png"), bbox_inches="tight", dpi=300)
    plt.savefig(figure_dir("earthquake_location"+title_suffixes+".pdf"), bbox_inches="tight", dpi=300)
    plt.show()

def source_link(loc_df, misfit_df,model, test, title_suffixes=""):
    pyocto_events = loc_df
    
    catalog_label="Standard"
    figure_dir = lambda x: "/media/wzm/HLPBook/DiTingProject/results/synthesis_data/"+model+f"/{test}/"+x
    model_label = model+title_suffixes
    pyocto_events.rename(columns={"loc_lon":"longitude" , "loc_lat": "latitude","loc_dep":"depth"}, inplace=True)
    stations = pd.read_csv("/media/wzm/HLPBook/DiTingProject/data/synthesis_data/pyocto/station.csv", sep='\s+')
    stations.rename(columns={"elevation": "elevation_m",}, inplace=True)
    standard_catalog = pd.read_csv("/media/wzm/HLPBook/DiTingProject/data/synthesis_data/"+test+"_standard_catalog.csv")

    config = {  "ylim_degree":(23.6, 24.15),
        "xlim_degree":(114.2, 114.8),
        "z(km)":(0, 22),
    }

    fig = plt.figure(figsize=plt.rcParams["figure.figsize"]*np.array([1.5,1]))
    box = dict(boxstyle='round', facecolor='white', alpha=1)
    text_loc = [0.05, 0.92]
    grd = fig.add_gridspec(ncols=2, nrows=2, width_ratios=[1.5, 1], height_ratios=[1,1])
    fig.add_subplot(grd[:, 0])
    for index ,row in misfit_df.iterrows():
        plt.plot([row["src_lon"], row["loc_lon"]], [row["src_lat"], row["loc_lat"]], color='green',alpha=0.4)
    plt.plot(pyocto_events["longitude"], pyocto_events["latitude"], '^',markersize=2, alpha=1.0)
    if standard_catalog is not None:
        plt.plot(standard_catalog["longitude"], standard_catalog["latitude"], '.', markersize=2, alpha=0.6,color='orange')
    plt.axis("scaled")

    # plt.xlim(np.array(config["xlim_degree"]))
    # plt.ylim(np.array(config["ylim_degree"]))
    # plt.ylim(pyocto_events["latitude"].min()-0.2,pyocto_events["latitude"].max()+0.2)
    # plt.xlim(pyocto_events["longitude"].min()-0.2,pyocto_events["longitude"].max()+0.2)
    plt.xlabel("Latitude")
    plt.ylabel("Longitude")
    plt.gca().set_prop_cycle(None)
    plt.plot([], [], '.', markersize=10, label=f"{model_label}", rasterized=True)
    plt.plot([], [], '.', markersize=10, label=f"{catalog_label}", rasterized=True)
    plt.plot(stations["longitude"], stations["latitude"], 'k', markersize=1, alpha=0.7, label="Stations")
    plt.legend(loc="lower right")
    plt.text(text_loc[0], text_loc[1], '(i)', horizontalalignment='left', verticalalignment="top", 
            transform=plt.gca().transAxes, fontsize="large", fontweight="normal", bbox=box)

    fig.add_subplot(grd[0, 1])
    for index ,row in misfit_df.iterrows():
        if row["src_dep_km"] %10 ==0:
            plt.plot([row["src_lon"], row["loc_lon"]], [row["src_dep_km"], row["loc_dep"]], color='green',alpha=0.4)
    plt.plot(pyocto_events["longitude"], pyocto_events["depth"], '.', markersize=2, alpha=1.0, rasterized=True)
    if standard_catalog is not None:
        plt.plot(standard_catalog["longitude"], standard_catalog["depth_km"], '.', markersize=2, alpha=0.6, rasterized=True , color='orange')
    # plt.xlim(np.array(config["xlim_degree"])+np.array([0.2,-0.27]))
    plt.ylim(config["z(km)"])
    plt.gca().invert_yaxis()
    plt.xlabel("Longitude")
    plt.ylabel("Depth (km)")
    plt.gca().set_prop_cycle(None)
    plt.plot([], [], '.', markersize=10, label=f"{model_label}")
    plt.plot([], [], '.', markersize=10, label=f"{catalog_label}")

    plt.legend(loc="lower right")
    plt.text(text_loc[0], text_loc[1], '(ii)', horizontalalignment='left', verticalalignment="top", 
            transform=plt.gca().transAxes, fontsize="large", fontweight="normal", bbox=box)

    fig.add_subplot(grd[1, 1])
    plt.plot(pyocto_events["latitude"], pyocto_events["depth"], '.', markersize=2, alpha=1.0, rasterized=True)
    if standard_catalog is not None:
        plt.plot(standard_catalog["latitude"], standard_catalog["depth_km"], '.', markersize=2, alpha=0.6, rasterized=True)
    # plt.xlim(np.array(config["ylim_degree"])+np.array([0.2,-0.27]))
    plt.ylim(config["z(km)"])
    plt.gca().invert_yaxis()
    plt.xlabel("Latitude")
    plt.ylabel("Depth (km)")
    plt.gca().set_prop_cycle(None)
    plt.plot([], [], '.', markersize=10, label=f"{model_label}")
    plt.plot([], [], '.', markersize=10, label=f"{catalog_label}")
    plt.legend(loc="lower right")
    plt.tight_layout()
    plt.text(text_loc[0], text_loc[1], '(iii)', horizontalalignment='left', verticalalignment="top", 
            transform=plt.gca().transAxes, fontsize="large", fontweight="normal", bbox=box)
    plt.savefig(figure_dir("earthquake_location_link"+title_suffixes+".png"), bbox_inches="tight", dpi=300)
    plt.savefig(figure_dir("earthquake_location_link"+title_suffixes+".pdf"), bbox_inches="tight", dpi=300)
    plt.show()



def plot_all(events,model,test , title_suffixes):
    compare_catalog(events,model,test , title_suffixes)
    plot_source_compare(events,model,test , title_suffixes)
    events.rename(columns={"longitude": "loc_lon",  "latitude":"loc_lat","depth":"loc_dep"}, inplace=True)
    misfit_df = plot_polar_misfit(events,model,test , title_suffixes)
    source_link(events,misfit_df,model,test , title_suffixes)