<jupyter_start><jupyter_text>PART 2: Remove trend and seasonality by least-squares modeling This script aims to:- Least squares model to remove the tidal and seasonal signals as well as co-seismic jumps- Earthquakes selected are obtained from wilber 3 web interface- Co-seismic jumps from selected earthquakes (>5.1 Mw)  Required:- python- pandas- jupyter- notebook- matplotlib- seabornthis should be easy to set up in a conda env: conda create -c conda-forge -n dtwclustering python=3.7 pandas numpy jupyter notebook matplotlib seaborn__Author: <NAME> @Institute of Earth Sciences, Academia Sinica__  Import Libraries<jupyter_code>import pandas as pd
import numpy as np
from scipy.optimize import least_squares as ls
from numpy import linalg as LA
import matplotlib.pyplot as plt
from datetime import datetime, timedelta
from dtwhaclustering.analysis_support import toYearFraction as tyf
from dtwhaclustering.leastSquareModeling import lsqmodeling

import matplotlib
from matplotlib import style
import scipy.io as sio
import os, glob
from IPython.display import Image, display
from random import randrange 
## default matplotlib parameters
import matplotlib
font = {'family' : 'Times',
        'weight' : 'bold',
        'size'   : 22}

matplotlib.rc('font', **font)
plt.rcParams["figure.figsize"] = (12, 6)
plt.style.use('ggplot')
matplotlib.rcParams['pdf.fonttype'] = 42<jupyter_output><empty_output><jupyter_text>Load data<jupyter_code>#load pickle data
dataloc = "pickleFiles"
dUU=pd.read_pickle(os.path.join(dataloc,"dU_data.pickle"))
dNN=pd.read_pickle(os.path.join(dataloc,"dN_data.pickle"))
dEE=pd.read_pickle(os.path.join(dataloc,"dE_data.pickle"))
dUU.head()<jupyter_output><empty_output><jupyter_text>Least square modeling for all the stations for three-components  Remove seasonlity only (keep the trend and co-seismic jumps)<jupyter_code>if not os.path.exists(os.path.join(dataloc,"dU_wo_seasn.pickle")) or not os.path.exists(os.path.join(dataloc,"dN_wo_seasn.pickle")) or not os.path.exists(os.path.join(dataloc,"dE_wo_seasn.pickle")):
    final_dU, final_dN, final_dE = lsqmodeling(dUU, dNN, dEE,stnlocfile="helper_files/stn_loc.txt",  plot_results=True, remove_trend=False, remove_seasonality=True, remove_jumps=False)
else:
    final_dU=pd.read_pickle(os.path.join(dataloc,"dU_wo_seasn.pickle"))
    final_dN=pd.read_pickle(os.path.join(dataloc,"dN_wo_seasn.pickle"))
    final_dE=pd.read_pickle(os.path.join(dataloc,"dE_wo_seasn.pickle"))
    final_dU.head()<jupyter_output>Plotting LSQ (U): PKGM: : 115it [02:54,  1.46it/s]
Plotting LSQ (N): SHJU: : 115it [03:19,  1.22it/s]
Plotting LSQ (E): GS16: : 115it [03:39,  1.20it/s]<jupyter_text>Check the least square modeling results<jupyter_code>all_U_png = glob.glob(os.path.join("LSQOut","*_U.png"))
all_N_png = glob.glob(os.path.join("LSQOut","*_N.png"))
all_E_png = glob.glob(os.path.join("LSQOut","*_E.png"))

# all_U_png = glob.glob(os.path.join("LSQOut","*SLNP_U.png"))
# all_N_png = glob.glob(os.path.join("LSQOut","*SLNP_N.png"))
# all_E_png = glob.glob(os.path.join("LSQOut","*SLNP_E.png"))


if all_U_png:
    display(Image(filename=all_U_png[randrange(0, len(all_U_png))]))   
else:
    print("No LS figures available")
if all_N_png:
    display(Image(filename=all_N_png[randrange(0, len(all_N_png))])) 
if all_E_png:
    display(Image(filename=all_E_png[randrange(0, len(all_E_png))]))<jupyter_output><empty_output><jupyter_text>Save lsq modeled result in pickle and mat format<jupyter_code>if not os.path.exists(os.path.join(dataloc,"dU_wo_seasn.pickle")) or not os.path.exists(os.path.join(dataloc,"dN_wo_seasn.pickle")) or not os.path.exists(os.path.join(dataloc,"dE_wo_seasn.pickle")):

    final_dU.to_pickle(os.path.join(dataloc,"dU_wo_seasn.pickle"))
    final_dN.to_pickle(os.path.join(dataloc,"dN_wo_seasn.pickle"))
    final_dE.to_pickle(os.path.join(dataloc,"dE_wo_seasn.pickle"))

# read station information
station_loc_file="helper_files/stn_loc.txt"
stnloc = pd.read_csv(station_loc_file, header=None,
                          sep='\s+', names=['stn', 'lon', 'lat'])
stnloc.set_index('stn', inplace=True)
new_stn_info_list = []
for stn in final_dU.columns.values:
    stn = stn.split("_")[0]
    _dict={}
    _dict['stn'] = stn
    _dict['lon'] = stnloc.loc[stn,'lon']
    _dict['lat'] = stnloc.loc[stn,'lat']
    new_stn_info_list.append(_dict)
    
stn_info_df = pd.DataFrame(new_stn_info_list)
stn_info_df.head()
# if not os.path.exists(os.path.join(dataloc,'all_data_wo_seasn.mat')):
#     sio.savemat(os.path.join(dataloc,'all_data_wo_seasn.mat'), {'slat': stn_info_df['lat'].values, 'slon': stn_info_df['lon'].values, 'tdata': final_dU.index.values,
#                                         'stns': stn_info_df['stn'].values, 'dU': np.array(final_dU), 'dN': np.array(final_dN), 'dE': np.array(final_dE)})

# # Saving the slope info in the mat file
# stn_slope = pd.read_csv('stn_slope_res_U.txt', sep='\s+',
#                         header=None, names=['stn', 'lon', 'lat', 'slope'])
# sio.savemat('slope_info.mat', {'slat': stn_slope['lat'].values, 'slon': stn_slope['lon'].values,
#                                'stns': stn_slope['stn'].values, 'slope': stn_slope['slope'].values})<jupyter_output><empty_output><jupyter_text>Compute dN, dE and dU without seasonality, trend and jumps- to compute the CME (common mode error)<jupyter_code>if not os.path.exists(os.path.join(dataloc,"dU_wo_seasn_trend_jump.pickle")) or not os.path.exists(os.path.join(dataloc,"dN_wo_seasn_trend_jump.pickle")) or not os.path.exists(os.path.join(dataloc,"dE_wo_seasn_trend_jump.pickle")):
    final_dU_2, final_dN_2, final_dE_2 = lsqmodeling(dUU, dNN, dEE,stnlocfile="helper_files/stn_loc.txt",  plot_results=False, remove_trend=True, remove_seasonality=True, remove_jumps=True)
else:
    final_dU_2=pd.read_pickle(os.path.join(dataloc,"dU_wo_seasn_trend_jump.pickle"))
    final_dN_2=pd.read_pickle(os.path.join(dataloc,"dN_wo_seasn_trend_jump.pickle"))
    final_dE_2=pd.read_pickle(os.path.join(dataloc,"dE_wo_seasn_trend_jump.pickle"))
    final_dU_2.head()
if not os.path.exists(os.path.join(dataloc,"dU_wo_seasn_trend_jump.pickle")) or not os.path.exists(os.path.join(dataloc,"dN_wo_seasn_trend_jump.pickle")) or not os.path.exists(os.path.join(dataloc,"dE_wo_seasn_trend_jump.pickle")):
    final_dU_2.to_pickle(os.path.join(dataloc,"dU_wo_seasn_trend_jump.pickle"))
    final_dN_2.to_pickle(os.path.join(dataloc,"dN_wo_seasn_trend_jump.pickle"))
    final_dE_2.to_pickle(os.path.join(dataloc,"dE_wo_seasn_trend_jump.pickle"))

# read station information
station_loc_file="helper_files/stn_loc.txt"
stnloc = pd.read_csv(station_loc_file, header=None,
                          sep='\s+', names=['stn', 'lon', 'lat'])
stnloc.set_index('stn', inplace=True)
new_stn_info_list = []

# east_coast_stations = ['PEPU', 'DAJN', 'NDHU', 'CHUN', 'SHUL', 'TUNH', 'DAWU', 'CHGO', 'YENL', 'SHAN', 'SOFN', 'TAPE', 'ERPN', 'CHEN', 'TAPO', 'SINL', 'LONT', 'JULI', 'JSUI', 'TTUN', 'NAAO', 'SPAO', 'MOTN', 'SLNP', 'WARO', 'SLIN', 'WULU']
for stn in final_dU.columns.values:
    stn = stn.split("_")[0]
    
    _dict={}
    _dict['stn'] = stn
    _dict['lon'] = stnloc.loc[stn,'lon']
    _dict['lat'] = stnloc.loc[stn,'lat']
    new_stn_info_list.append(_dict)
    
stn_info_df = pd.DataFrame(new_stn_info_list)
stn_info_df.head()
stn_info_df.to_csv('helper_files/selected_stations_info.txt', index=False)
# if not os.path.exists(os.path.join(dataloc,'all_data_wo_seasn_trend_jump.mat')):
#     sio.savemat(os.path.join(dataloc,'all_data_wo_seasn_trend_jump.mat'), {'slat': stn_info_df['lat'].values, 'slon': stn_info_df['lon'].values, 'tdata': final_dU_2.index.values,
#                                         'stns': stn_info_df['stn'].values, 'dU': np.array(final_dU_2), 'dN': np.array(final_dN_2), 'dE': np.array(final_dE_2)})<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15528901365318362, 'fraction_numerical': 0.008166390200331759, 'mean_word_length': 5.6199324324324325, 'pattern_counts': {'":': 0, '<': 25, '<?xml version=': 0, '>': 26, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '641269', 'license': 'Apache-2.0', 'n_tokens_mistral': 3342, 'n_tokens_neox': 3252, 'n_words': 488, 'path': 'example_notebooks/.ipynb_checkpoints/part2_remove_trend_seasonality_least_square_modeling-checkpoint.ipynb', 'repo_name': 'earthinversion/Dynamic-Time-Warping-based-Hierarchical-Agglomerative-Clustering', 'size': '7842'}