<jupyter_start><jupyter_text>Analysis on correlation between earthquakes and surface mass loading<jupyter_code>import xarray as xr
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
from global_land_mask import globe
import scipy.stats as stats
import os
import importlib
import geopandas as gpd

# file containing functions:
import loadquakes
importlib.reload(loadquakes)<jupyter_output><empty_output><jupyter_text>Load datasets (GRACE and declustered catalog)<jupyter_code>dataset_path = '/data/fast1/surface-load/'
ds = xr.open_dataset(dataset_path + "GRCTellus.JPL.200204_202201.GLO.RL06M.MSCNv02CRI.nc")

# Events from 2002-04-16 to 2022-1-16 with magnitude >= 5.4
catalog = pd.read_csv('ind_and_bg_probs.csv')
catalog["time"] = catalog["time"].astype("datetime64")
# This is the catalog with all details from the USGS Earthquake Catalog

original_catalog = pd.read_csv('19970101-20220116.csv')
original_catalog = original_catalog[original_catalog['time']>'2002-04-16']
original_catalog = original_catalog[original_catalog['time']<'2022-01-16']
cat = original_catalog.reset_index()
depth = cat['depth']
catalog['depth'] = depth

catalog<jupyter_output><empty_output><jupyter_text>Filter events that that have less than 90% chance of being triggered The mean probability of being triggered (P_triggered) for the catalog is 97.97% with a standard deviation of 9.25% . Here, we examine events that have less 90% probability of being triggered by another earthquake<jupyter_code>mainshock = catalog[catalog['P_triggered']<.9]
mainshock.sort_values("magnitude")<jupyter_output><empty_output><jupyter_text>Add surface load to the earthquake dataframes<jupyter_code>def get_mass(eq):
    mass_change = ds.sel(lat=(eq["latitude"]),lon=(eq["longitude"]),time=(eq["time"]), method="nearest")
    mass_change = mass_change["lwe_thickness"]
    return mass_change.data
  
mainshock['dmass'] = mainshock.apply(lambda row: get_mass(row), axis=1)
catalog['dmass'] = catalog.apply(lambda row: get_mass(row), axis=1)<jupyter_output><ipython-input-6-62d8ffe61e8d>:6: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  mainshock['dmass'] = mainshock.apply(lambda row: get_mass(row), axis=1)<jupyter_text>Are the mean/meadian values different between: 1. The whole dataset, 2. times with earthquakes, and 3. times with earthquakes that have low probability of being triggered?The answer appears to be: no, not really. The Bayesian analysis explains why! It's because mainshocks are favored at large positive AND negative surface loads. The symmetry implies that the means don't change much.<jupyter_code>print('Mainshock mean:   %f'%mainshock['dmass'].mean())
print('Mainshock median: %f'%mainshock['dmass'].median())

print('Earthquake mean:   %f'%catalog['dmass'].mean())
print('Earthquake median:   %f'%catalog['dmass'].median())

print(ds.lwe_thickness.mean())
print(ds.lwe_thickness.median())<jupyter_output>Mainshock mean:   0.441004
Mainshock median: 0.284740
Earthquake mean:   0.827411
Earthquake median:   1.054229
<xarray.DataArray 'lwe_thickness' ()>
array(-0.65387518)
<xarray.DataArray 'lwe_thickness' ()>
array(0.7491707)<jupyter_text>Correlation between all events and surface mass loading<jupyter_code># takes a long time to run (has already been run for 2002-04-16 to 2022-01-16)
# loadquakes.calc_stats(np.array(catalog.dmass),np.array(ds['lwe_thickness'].data.flatten()))
importlib.reload(loadquakes)

def set_of_figures_load(all_time, earthquake_only,bayes_title,method):

    fig,(ax1,ax2,ax3) = plt.subplots(3,1, figsize=(7,14))

    plt.style.use('fivethirtyeight')
    loadquakes.plot_hist(all_time, earthquake_only, ax1, ax2, 
              'a. Cumulative Distribution', 'b. Probability Density', method)

    loadquakes.plot_bayes(all_time, earthquake_only, ax3, bayes_title,
                         method)

    fig.tight_layout()
loadquakes.freedman_diaconis(ds['lwe_thickness'].data.flatten(), returnas = "bins")
set_of_figures_load(ds['lwe_thickness'].data.flatten(), catalog.dmass,
                    'c. Bayesian Probability of Event (n=12573)','fd')
set_of_figures_load(ds['lwe_thickness'].data.flatten(), catalog.dmass,
                    'c. Bayesian Probability of Event (n=12573)','Sturge')<jupyter_output>32
32
31<jupyter_text>Correlation between mainshocks and surface mass loading  Both statistical tests produce p-values < 0.05, which indicate that the two samples are not drawn from the same distribution. This means that the distribution of surface mass load during earthquake occurrence is significantly different from background.<jupyter_code># takes a long time to run (has already been run for 2002-04-16 to 2022-01-16)
# loadquakes.calc_stats(np.array(mainshock.dmass),np.array(ds['lwe_thickness'].data.flatten()))
importlib.reload(loadquakes)

set_of_figures_load(ds['lwe_thickness'].data.flatten(), mainshock.dmass,
                    'c. Bayesian Probability of Event (n=537)','fd')<jupyter_output>/home/amand4/.local/lib/python3.9/site-packages/scipy/stats/stats.py:2836: VisibleDeprecationWarning: use of scale='raw' is deprecated, use scale=1.0 instead
  warnings.warn(<jupyter_text>The conditional probability will change quantitatively depending on the bin size. However, the same qualitative distribution remains, where there is a higher relative conditional probability during periods of large loading and unloading.  Here, we calculate the Bayesian probability using two different bin sizes which are calculated via two methods: a. Sturge's Rule b. Freedman-Diaconis Rule<jupyter_code>importlib.reload(loadquakes)

fig,(ax1,ax2) = plt.subplots(1,2, figsize=(17,6))

plt.style.use('fivethirtyeight')
    
loadquakes.plot_bayes(ds['lwe_thickness'].data.flatten(), mainshock.dmass, ax1,
                      "a. Bin size calculated via Sturge's Rule", method='Sturge')

loadquakes.plot_bayes(ds['lwe_thickness'].data.flatten(), mainshock.dmass, ax2,
                      "b. Bin size calculated via Freedman-Diaconis Rule", method='fd')

fig.tight_layout()<jupyter_output>21
20<jupyter_text>Mapping earthquakes with high conditional probability  First, calculate the conditional probability of each mainshock<jupyter_code>conditional_probability = loadquakes.get_cond_probability(
                            np.array(ds['lwe_thickness'].data.flatten()),
                            np.array(mainshock.dmass), np.array(mainshock.dmass), 'fd')
mainshock['cp'] = conditional_probability
fig,ax=plt.subplots()
plt.hist(np.array(mainshock['cp']),11)
plt.title('Histogram of mainshock conditional probabilities')
plt.yscale('log')
from scipy.stats import iqr
iqr(mainshock['cp'])
np.median(mainshock.cp)
importlib.reload(loadquakes)

def all_figs(all_loads,full_catalog,declustered_catalog,bayes_title1,bayes_title2,method):

    fig,((ax1,ax4),(ax2,ax5),(ax3,ax6)) = plt.subplots(3,2, figsize=(14,14))

    plt.style.use('fivethirtyeight')
    
    loadquakes.plot_hist(all_loads, full_catalog, ax1, ax2, 
              'a. Cumulative Distribution', 'b. Probability Density', method)

    loadquakes.plot_bayes(all_loads, full_catalog, ax3, bayes_title1,
                         method)
    
    loadquakes.plot_hist(all_loads, declustered_catalog, ax4, ax5, 
              'd. Cumulative Distribution (declustered catalog)', 
                         'e. Probability Density (declustered catalog)',
                         method)

    loadquakes.plot_bayes(all_loads, declustered_catalog, ax6, bayes_title2,
                         method)

    fig.tight_layout()
    
all_figs(ds['lwe_thickness'].data.flatten(),
         catalog.dmass,
         mainshock.dmass,
         'c. Bayesian Probability of Event (n=12573)',
         'f. Bayesian Probability of Event (n=537)','fd')
np.min(mainshock.dmass)<jupyter_output><empty_output><jupyter_text>The relative conditional probability for the declustered catalog has a median of is 1.06 and an IQR of 0.68. We use the addition of these values as the threshold to categorize an event as a SLIQ.<jupyter_code>sliq = mainshock[mainshock['cp']>1.74]
len(sliq)
sliq.sort_values('cp')
sliq.sort_values('magnitude')
from mpl_toolkits.axes_grid1 import make_axes_locatable
import matplotlib.cm as cm
# c = [0.5, 0.5, 0.5]
def probability_map_cb(events,color,label):

    gdf = gpd.GeoDataFrame(events,
                           geometry=gpd.points_from_xy(events.longitude, 
                                                       events.latitude))
    world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
    ax = world.plot(color='white', edgecolor='black', figsize=(15,10))
    divider = make_axes_locatable(ax)
    cax = divider.append_axes("bottom", size="5%", pad=0.6)
    ax.scatter(catalog.longitude, catalog.latitude, c="darkgrey", s=20,marker=".")
    cmap = cm.get_cmap('viridis', 11)    # 11 discrete colors
    gdf.plot(ax=ax,cax=cax,alpha=0.5,column=color,cmap=cmap,legend=True,
             edgecolor='k',
             markersize = 1e-5*(events.magnitude)**10,
             legend_kwds={'label': "Relative conditional probability of event",
                            'orientation': "horizontal"},
            vmax=5,
            vmin=1.74)
    gdf.plot(ax=ax,facecolor="None",
         edgecolor='k',
         markersize = 1e-5*(events.magnitude)**10)
    ax.set_xlabel('Longitude', fontsize = 15)
    ax.set_ylabel("Latitude", fontsize = 15)
    ax.set_title(label)
    plt.show()

probability_map_cb(sliq,sliq.cp,
                'SLIQs with relative conditional probability > 1.74')
montana = sliq[sliq['longitude']<-100]
montana<jupyter_output><empty_output><jupyter_text>Next, we examine the depth of a SLIQ compared to the depth of all events in the full catalog and declustered catalog.<jupyter_code># fig,(ax1,ax2) = plt.subplots(1,2, figsize=(15,5))

def depth_fig(ax,catalog,cumulative,label,title):

    bins_depth = loadquakes.calculate_bin_sizes(catalog.depth,'fd')

    ax.hist(sliq.depth,bins,density = True,cumulative=cumulative, histtype='step',
            label='SLIQs',linewidth=1.5)
    ax.hist(catalog.depth, 50,density = True, cumulative=cumulative,histtype='step',
            label=label,linewidth=1.5)
    yl = ax.get_ylim()
    ax.set_ylim((-0.01,1.4*yl[1]))
    xl = ax.get_xlim()
    ax.set_xlim(xl[0],xl[1])
    ax.legend()
    ax.set_xlabel('Depth (km)', fontsize = 17)
    ax.set_ylabel("Cumulative probability", fontsize = 17)
    ax.set_title(title, fontsize = 17)

fig,((ax1,ax2),(ax3,ax4)) = plt.subplots(2,2, figsize=(15,11))

depth_fig(ax1,mainshock,True,'Declustered catalog',"Cumulative Distribution")
depth_fig(ax2,mainshock,False,'Declustered catalog','Probability Density')
depth_fig(ax3,catalog,True,'Full catalog',"Cumulative Distribution")
depth_fig(ax4,catalog,False,'Full catalog','Probability Density')

fig.tight_layout()<jupyter_output>/home/amand4/.local/lib/python3.9/site-packages/scipy/stats/stats.py:2836: VisibleDeprecationWarning: use of scale='raw' is deprecated, use scale=1.0 instead
  warnings.warn(
/home/amand4/.local/lib/python3.9/site-packages/scipy/stats/stats.py:2836: VisibleDeprecationWarning: use of scale='raw' is deprecated, use scale=1.0 instead
  warnings.warn(
/home/amand4/.local/lib/python3.9/site-packages/scipy/stats/stats.py:2836: VisibleDeprecationWarning: use of scale='raw' is deprecated, use scale=1.0 instead
  warnings.warn(
/home/amand4/.local/lib/python3.9/site-packages/scipy/stats/stats.py:2836: VisibleDeprecationWarning: use of scale='raw' is deprecated, use scale=1.0 instead
  warnings.warn(<jupyter_text>Depth statistics:<jupyter_code>print('Catalog mean:   %f'%catalog['depth'].mean())
print('Catalog median:   %f'%catalog['depth'].median())

print('Declustered catalog mean:   %f'%mainshock['depth'].mean())
print('Declustered catalog median: %f'%mainshock['depth'].median())

print('SLIQ mean:   %f'%sliq['depth'].mean())
print('SLIQ median:   %f'%sliq['depth'].median())
sliq.to_csv('sliq.csv', index=False)<jupyter_output><empty_output><jupyter_text>Correlation between earthquakes and random data *mean and std deviation of random array same as GRACE data<jupyter_code>np.mean(ds['lwe_thickness'])
np.std(ds['lwe_thickness'])
random_load1 = np.random.normal(loc=-0.65387518, scale=25.49255185, size=(205, 360, 720))
ds_random = ds
ds_random = ds_random.assign(load=ds_random['lwe_thickness'] - ds_random['lwe_thickness'] + random_load1)
# earthquake catalog

mainshock_random = catalog[catalog['P_triggered']<.9]
def get_mass_random(eq):
    mass_change = ds_random.sel(lat=(eq["latitude"]),lon=(eq["longitude"]),time=(eq["time"]), method="nearest")
    mass_change = mass_change["load"]
    return mass_change.data
  
mainshock_random['dmass'] = mainshock_random.apply(lambda row: get_mass_random(row), axis=1)
set_of_figures_load(ds_random['load'].data.flatten(), mainshock_random.dmass,'c. Bayesian Probability of Mainshock (n=537)')
stats.ks_2samp(mainshock_random.dmass,ds_random['load'].data.flatten())<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12830988046603117, 'fraction_numerical': 0.03207747011650779, 'mean_word_length': 4.98415572657311, 'pattern_counts': {'":': 0, '<': 55, '<?xml version=': 0, '>': 54, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '484106', 'license': 'MIT', 'n_tokens_mistral': 4823, 'n_tokens_neox': 4408, 'n_words': 910, 'path': 'surface-load-induced-quakes.ipynb', 'repo_name': 'amandasyamsul/SLIM', 'size': '13218'}