(시도) Eearthquake
- imports
- load data and clean it
- define class
- analysis_df_global(2010~2015)
- analysis_df_global(2015~2020)
import tqdm
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import plotly.express as px
import warnings
warnings.simplefilter("ignore", np.ComplexWarning)
from haversine import haversine
from IPython.display import HTML
import rpy2
import rpy2.robjects as ro
from rpy2.robjects.vectors import FloatVector
from rpy2.robjects.packages import importr
- load
df= pd.read_csv('https://raw.githubusercontent.com/plotly/datasets/master/earthquakes-23k.csv')
df
df_korea= pd.read_csv('earthquake_korea2.csv').iloc[:,[1,2,5,6]].rename(columns={'규모':'Magnitude'})
df_global= pd.concat([pd.read_csv('00_05.csv'),pd.read_csv('05_10.csv'),pd.read_csv('10_15.csv'),pd.read_csv('15_20.csv')]).iloc[:,[0,1,2,4]].rename(columns={'latitude':'Latitude','longitude':'Longitude','mag':'Magnitude'}).reset_index().iloc[:,1:]
- cleaning
df.Date[df.Date == '1975-02-23T02:58:41.000Z']
df.iloc[3378,0] = '02/03/1975'
df.Date[df.Date == '1985-04-28T02:53:41.530Z']
df.iloc[7512,0] = '04/28/1985'
df.Date[df.Date == '2011-03-13T02:23:34.520Z']
df.iloc[20650,0] = '03/13/2011'
df= df.assign(Year=list(map(lambda x: x.split('/')[-1], df.Date))).iloc[:,1:]
df
df.Year = df.Year.astype(np.float64)
df_korea = df_korea.assign(Year=list(map(lambda x: x.split('/')[0], df_korea.발생시각))).iloc[:,1:]
df_korea = df_korea.assign(Latitude=list(map(lambda x: x.split(' ')[0], df_korea.위도))).iloc[:,[0,2,3,4]]
df_korea = df_korea.assign(Longitude=list(map(lambda x: x.split(' ')[0], df_korea.경도))).iloc[:,[0,2,3,4]]
df_global = df_global.assign(Year=list(map(lambda x: x.split('-')[0], df_global.time))).iloc[:,1:]
df_korea.Year = df_korea.Year.astype(np.float64)
df_korea.Latitude = df_korea.Latitude.astype(np.float64)
df_korea.Longitude = df_korea.Longitude.astype(np.float64)
df_global.Year = df_global.Year.astype(np.float64)
class MooYaHo:
def __init__(self,df):
self.df = df
self.f = df.Magnitude.to_numpy()
self.year = df.Year.to_numpy()
self.lat = df.Latitude.to_numpy()
self.long = df.Longitude.to_numpy()
self.n = len(self.f)
self.theta= None
def get_distance(self):
self.D = np.zeros([self.n,self.n])
locations = np.stack([self.lat, self.long],axis=1)
for i in tqdm.tqdm(range(self.n)):
for j in range(i,self.n):
self.D[i,j]=haversine(locations[i],locations[j])
self.D = self.D+self.D.T
def get_weightmatrix(self,theta=1,beta=0.5,kappa=4000):
self.theta = theta
dist = np.where(self.D<kappa,self.D,0)
self.W = np.exp(-(dist/self.theta)**2)
def _eigen(self):
d= self.W.sum(axis=1)
D= np.diag(d)
self.L = np.diag(1/np.sqrt(d)) @ (D-self.W) @ np.diag(1/np.sqrt(d))
self.lamb, self.Psi = np.linalg.eigh(self.L)
self.Lamb = np.diag(self.lamb)
def fit(self,m):
self._eigen()
self.fhat = self.Psi[:,0:m]@self.Psi[:,0:m].T@self.f
self.df = self.df.assign(MagnitudeHat = self.fhat)
self.df = self.df.assign(Residual = self.df.Magnitude- self.df.MagnitudeHat)
plt.plot(self.f,'.')
plt.plot(self.fhat,'x')
def vis(self,MagThresh=7,ResThresh=1):
fig = px.density_mapbox(self.df,
lat='Latitude',
lon='Longitude',
z='Magnitude',
radius=5,
center=dict(lat=37, lon=160),
zoom=1.5,
height=900,
opacity = 0.4,
mapbox_style="stamen-terrain",
range_color=[-7,7])
fig.update_layout(margin={"r":0,"t":0,"l":0,"b":0})
fig.add_scattermapbox(lat = self.df.query('Magnitude > @MagThresh')['Latitude'],
lon = self.df.query('Magnitude > @MagThresh')['Longitude'],
text = self.df.query('Magnitude > @MagThresh')['Magnitude'],
marker_size= 8,
marker_color= 'red',
opacity = 0.6
)
fig.add_scattermapbox(lat = self.df.query('Residual**2 > @ResThresh')['Latitude'],
lon = self.df.query('Residual**2 > @ResThresh')['Longitude'],
text = self.df.query('Magnitude > @ResThresh')['Magnitude'],
marker_size= 8,
marker_color= 'blue',
opacity = 0.5
)
return HTML(fig.to_html(include_mathjax=False, config=dict({'scrollZoom':False})))
def visf(self):
fig = px.density_mapbox(self.df,
lat='Latitude',
lon='Longitude',
z='Magnitude',
radius=5,
center=dict(lat=37, lon=160),
zoom=1.5,
height=900,
opacity = 0.7,
mapbox_style="stamen-terrain",
range_color=[-7,7])
fig.update_layout(margin={"r":0,"t":0,"l":0,"b":0})
return HTML(fig.to_html(include_mathjax=False, config=dict({'scrollZoom':False})))
def visfhat(self):
fig = px.density_mapbox(self.df,
lat='Latitude',
lon='Longitude',
z='MagnitudeHat',
radius=5,
center=dict(lat=37, lon=160),
zoom=1.5,
height=900,
opacity = 0.7,
mapbox_style="stamen-terrain",
range_color=[-7,7])
fig.update_layout(margin={"r":0,"t":0,"l":0,"b":0})
return HTML(fig.to_html(include_mathjax=False, config=dict({'scrollZoom':False})))
def visres(self,MagThresh=7,ResThresh=1):
fig = px.density_mapbox(self.df,
lat='Latitude',
lon='Longitude',
z=[0] * len(self.df),
radius=5,
center=dict(lat=37, lon=160),
zoom=1.5,
height=900,
opacity = 0.7,
mapbox_style="stamen-terrain",
range_color=[-7,7])
fig.update_layout(margin={"r":0,"t":0,"l":0,"b":0})
fig.add_scattermapbox(lat = self.df.query('Residual**2 > @ResThresh')['Latitude'],
lon = self.df.query('Residual**2 > @ResThresh')['Longitude'],
text = self.df.query('Magnitude > @ResThresh')['Magnitude'],
marker_size= 8,
marker_color= 'blue',
opacity = 0.7
)
return HTML(fig.to_html(include_mathjax=False, config=dict({'scrollZoom':False})))
class MooYaHo2(MooYaHo): # ebayesthresh 기능추가
def fit2(self,ref=0.5): # fit with ebayesthresh
self._eigen()
self.fbar = self.Psi.T @ self.f # fbar := graph fourier transform of f
self.power = self.fbar**2
ebayesthresh = importr('EbayesThresh').ebayesthresh
self.power_threshed=np.array(ebayesthresh(FloatVector(self.fbar**2)))
self.fbar_threshed = np.where(self.power_threshed>0,self.fbar,0)
self.fhat = self.Psi@self.fbar_threshed
self.df = self.df.assign(MagnitudeHat = self.fhat)
self.df = self.df.assign(Residual = self.df.Magnitude- self.df.MagnitudeHat)
self.con = np.where(self.df.Residual>0.7,1,0)
#plt.plot(self.f,'.')
#plt.plot(self.fhat,'x')
# fig, axs = plt.subplots(2,2,figsize=(16,10))
# axs[0,0].plot(self.f,'b')
# axs[0,0].set_title('Magnitude')
# axs[0,0].set_ylim([4.5,9])
# axs[0,1].plot(self.fhat,'k')
# axs[0,1].set_title('MagnitudeHat')
# axs[0,1].set_ylim([4.5,9])
# axs[1,0].plot(self.con,'r*')
# axs[1,0].set_title('Residual square')
# axs[1,1].plot(self.f,'b')
# axs[1,1].plot(self.fhat,'k')
# axs[1,1].plot(self.con,'r*')
# axs[1,1].set_title('Graph')
# axs[1,1].set_ylim([4.5,9])
# plt.tight_layout()
# plt.show()
class MooYaHo3(MooYaHo2):
def vis(self,MagThresh=7,ResThresh=1):
fig = px.density_mapbox(self.df,
lat='Latitude',
lon='Longitude',
z='Magnitude',
radius=5,
center=dict(lat=37, lon=126),
zoom=5.7,
height=900,
opacity = 0.3,
mapbox_style="stamen-terrain")
fig.update_layout(margin={"r":0,"t":0,"l":0,"b":0})
fig.add_scattermapbox(lat = self.df.query('Magnitude > @MagThresh')['Latitude'],
lon = self.df.query('Magnitude > @MagThresh')['Longitude'],
text = self.df.query('Magnitude > @MagThresh')['Magnitude'],
marker_size= 8,
marker_color= 'red',
opacity = 0.5
)
fig.add_scattermapbox(lat = self.df.query('Residual**2 > @ResThresh')['Latitude'],
lon = self.df.query('Residual**2 > @ResThresh')['Longitude'],
text = self.df.query('Magnitude > @ResThresh')['Magnitude'],
marker_size= 8,
marker_color= 'blue',
opacity = 0.5
)
return HTML(fig.to_html(include_mathjax=False, config=dict({'scrollZoom':False})))
ebayesthresh = importr('EbayesThresh').ebayesthresh
- make instance for analysis
moo_global=MooYaHo2(df_global.query("2010 <= Year < 2015"))
- get distance
moo_global.get_distance()
moo_global.D[moo_global.D>0].mean()
plt.hist(moo_global.D[moo_global.D>0])
- weight matrix
moo_global.get_weightmatrix(theta=(8810.865423093777),kappa=2500)
- fit
moo_global.fit2()
moo_global.df.sort_values("Residual",ascending=False).iloc[:40,:]
(2010~2014 시도)
- 21번째 Ouest Department, Haiti 아이티 지진 2010년 진도 7.0
- 24번쨰 Puchuncavi, Valparaíso, Chile 칠레 지진 2014년 진도 6.4
- 28번째 Baoxing County, Yaan, Sichuan, China 중국 쓰촨성 지진 2013년 진도 6.6
(2010~2015 시도_결과 좋지 않음?!)
- 23번째 2010년 West New Britain Province, Papua New Guinea 진도 7.3
- 24번째 2011년 Kuzawa Terayama, Tanagura, Higashishirakawa District, Fukushima 963-5671, Japan 진도 6.6
- 29번째 2015년 Kishim, Afghanistan 진도 7.5
- vis
- make instance for analysis
moo_global=MooYaHo2(df_global.query("2015 <= Year <= 2020"))
- get distance
moo_global.get_distance()
moo_global.D[moo_global.D>0].mean()
plt.hist(moo_global.D[moo_global.D>0])
- weight matrix
moo_global.get_weightmatrix(theta=(8814.318793468068),kappa=2500)
- fit
moo_global.fit2()
moo_global.df.sort_values("Residual",ascending=False).iloc[:30,:]
바다 아닌 거
- 8번째 2016년 Rotherham, New Zealand 뉴질랜드 카이코우라 지진 진도 7.8
- 9번째 2015년 Langkuru Utara, Pureman, Alor Regency, East Nusa Tenggara, Indonesia 수마트라 진도 6.5
- 15번째 2018년 Hela Province, Papua New Guinea 파푸아뉴기니 진도 7.5
- 20번째 2015년 Kalinchok, Nepal 네팔 진도 7.3
- 26번째 2019년 Coquimbo, Chile 칠레 코킴보주 진도 6.7
- vis
pd.read_html('https://en.wikipedia.org/wiki/Lists_of_21st-century_earthquakes',encoding='utf-8')[0].query('Magnitude<=7')# List of deadliest earthquakes
pd.read_html('https://en.wikipedia.org/wiki/Lists_of_21st-century_earthquakes',encoding='utf-8')[3] # Deadliest earthquakes by year
class eachlocation(MooYaHo3):
def haiti(self,MagThresh=7,ResThresh=1):
fig = px.density_mapbox(self.df,
lat='Latitude',
lon='Longitude',
z='Magnitude',
radius=5,
center=dict(lat=18.4430, lon=-72.5710),
zoom=5,
height=900,
opacity = 0.6,
mapbox_style="stamen-terrain",
range_color=[-7,7])
fig.update_layout(margin={"r":0,"t":0,"l":0,"b":0})
fig.add_scattermapbox(lat = self.df.query('Magnitude > @MagThresh')['Latitude'],
lon = self.df.query('Magnitude > @MagThresh')['Longitude'],
text = self.df.query('Magnitude > @MagThresh')['Magnitude'],
marker_size= 5,
marker_color= 'red',
opacity = 0.1
)
fig.add_scattermapbox(lat = self.df.query('Residual**2 > @ResThresh')['Latitude'],
lon = self.df.query('Residual**2 > @ResThresh')['Longitude'],
text = self.df.query('Magnitude > @ResThresh')['Magnitude'],
marker_size= 20,
marker_color= 'blue',
opacity = 1
)
fig.show()
def lquique(self,MagThresh=7,ResThresh=1):
fig = px.density_mapbox(self.df,
lat='Latitude',
lon='Longitude',
z='Magnitude',
radius=5,
center=dict(lat=-32.6953, lon=-71.4416),
zoom=5,
height=900,
opacity = 0.6,
mapbox_style="stamen-terrain",
range_color=[-7,7])
fig.update_layout(margin={"r":0,"t":0,"l":0,"b":0})
fig.add_scattermapbox(lat = self.df.query('Magnitude > @MagThresh')['Latitude'],
lon = self.df.query('Magnitude > @MagThresh')['Longitude'],
text = self.df.query('Magnitude > @MagThresh')['Magnitude'],
marker_size= 5,
marker_color= 'red',
opacity = 0.1
)
fig.add_scattermapbox(lat = self.df.query('Residual**2 > @ResThresh')['Latitude'],
lon = self.df.query('Residual**2 > @ResThresh')['Longitude'],
text = self.df.query('Magnitude > @ResThresh')['Magnitude'],
marker_size= 20,
marker_color= 'blue',
opacity = 1
)
fig.show()
def sichuan(self,MagThresh=7,ResThresh=1):
fig = px.density_mapbox(self.df,
lat='Latitude',
lon='Longitude',
z='Magnitude',
radius=3,
center=dict(lat=30.3080, lon=102.8880),
zoom=5,
height=900,
opacity = 0.6,
mapbox_style="stamen-terrain",
range_color=[-7,7])
fig.update_layout(margin={"r":0,"t":0,"l":0,"b":0})
fig.add_scattermapbox(lat = self.df.query('Magnitude > @MagThresh')['Latitude'],
lon = self.df.query('Magnitude > @MagThresh')['Longitude'],
text = self.df.query('Magnitude > @MagThresh')['Magnitude'],
marker_size= 5,
marker_color= 'red',
opacity = 0.1
)
fig.add_scattermapbox(lat = self.df.query('Residual**2 > @ResThresh')['Latitude'],
lon = self.df.query('Residual**2 > @ResThresh')['Longitude'],
text = self.df.query('Magnitude > @ResThresh')['Magnitude'],
marker_size= 20,
marker_color= 'blue',
opacity = 1
)
fig.show()
each_location=eachlocation(df_global.query("2010 <= Year < 2015"))
- get distance
each_location.get_distance()
each_location.D[each_location.D>0].mean()
plt.hist(each_location.D[each_location.D>0])
- weight matrix
each_location.get_weightmatrix(theta=(8810.865423093777),kappa=2500)
- fit
each_location.fit2()
each_location.haiti(MagThresh=6.9,ResThresh=0.5)
each_location.lquique(MagThresh=6.4,ResThresh=0.4)
each_location.sichuan(MagThresh=6.5,ResThresh=0.4)