(시도) Eearthquake
- imports
- load data and clean it
- define class
- analysis
- 한계점 (?)
- analysis_Korea
- analysis_df_korea
- analysis_Global
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
- 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_korea
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:]
df_global
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_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_korea
df_global = df_global.assign(Year=list(map(lambda x: x.split('-')[0], df_global.time))).iloc[:,1:]
df_global
df.Year.unique()
df.Year = df.Year.astype(np.float64)
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)
# nlst = self.df.groupby('Year').aggregate(len).Latitude.tolist()
# nlst = [0]+np.cumsum(nlst).tolist()
# beta_matrix = np.zeros([self.n,self.n])
# for i in range(len(nlst)-1):
# beta_matrix[nlst[i]:nlst[i+1],nlst[i]:nlst[i+1]] = beta
# self.W = self.W * beta_matrix
# def _get_laplacian(self):
# 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 _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',
animation_frame="Year",animation_group="Year",
radius=15,
center=dict(lat=37, lon=160),
zoom=1.5,
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'],
marker_size= 8,
marker_color= 'red',
opacity = 1
)
fig.add_scattermapbox(lat = self.df.query('Residual**2 > @ResThresh')['Latitude'],
lon = self.df.query('Residual**2 > @ResThresh')['Longitude'],
marker_size= 8,
marker_color= 'blue',
opacity = 1
)
return HTML(fig.to_html(include_mathjax=False, config=dict({'scrollZoom':False})))
- make instance for analysis
moo=MooYaHo(df.query("Year >= 2010"))
- 테스트용으로 2016년 자료만 수집
- get distance
moo.get_distance()
moo.D[moo.D>0].mean()
- 0이 아닌 거리의 평균은 8600정도?
plt.hist(moo.D[moo.D>0])
- 히스토그램결과 -> 2500보다 거리가 작으면 거의 같은 지역이라고 봐도 무방할듯
- weight matrix
moo.get_weightmatrix(theta=(8647.92),kappa=3000)
- 평균적으로 노드는
np.exp(-(dist/8647)**2)=np.exp(-1)=0.36787944117144233정도의 연결강도를 가진다. - 거리가 2500이하이면 weight를 1로 설정한다.
- fit
moo.fit(m=1200)
moo.df.query('Year==2010 & Latitude == 18.443')
- Residual ; m & theta 의 조합 2010년 이후 아이티 기준
- m: 1000, kappa: 2000, Residual: 1.227006
- m: 1000, kappa: 2500, Residual: 1.009784
- m: 800, kappa: 2750, Residual: 1.022774
- m: kappa: Residual:
- m: kappa: Residual:
- vis
- 전체자료는 underlying에 위치
- 지진강도가 MagThresh 이상인 자료는 붉은점으로 시각화함
- 이상치정도(우리가 제안하는 새로운 메져) ResThresh 이상인 자료는 파란점으로 시각화함
- 시간효과를 따로 모형화하지는 않음
pd.read_html('https://ko.wikipedia.org/wiki/%EC%A7%80%EC%A7%84_%EB%AA%A9%EB%A1%9D',encoding='utf-8')[0].iloc[[1,2,3],:] # 가장 피해가 컸던 지진
pd.read_html('https://ko.wikipedia.org/wiki/%EC%A7%80%EC%A7%84_%EB%AA%A9%EB%A1%9D',encoding='utf-8')[1].iloc[[2,3,8],:]# 가장 규모가 컸던 지진
pd.read_html('https://www.usgs.gov/programs/earthquake-hazards/science/20-largest-earthquakes-world',encoding='utf-8')[1].iloc[[2,3,5,7,9,10,17,18],:] # A list of the 20 largest earthquakes in the world.
- 예상: 특정 진도 이상이면서 residual 1보다 큰 곳이 의미있는 지점일 것이다.
- 1976년
- 중국 탕산 지진, 일어난 3곳 ( 2곳은 residual 1초과)
- 2004년
- 인도 수마트라섬 지진, 9 이상인 점(red)과 residual 1 초과인 점(blue)이 겹침
- 때때로 해구에 도달하는 지진에서 섭입판 경계면이 파열되어 1907년, 2004년, 2010년과 같은 대형 쓰나미가 발생한다. (https://ko.wikipedia.org/wiki/2022%EB%85%84_%EC%88%98%EB%A7%88%ED%8A%B8%EB%9D%BC_%EC%A7%80%EC%A7%84)
- 2010년
- 아이티지진(7.0)(지정한 MagThresh 기준은 초과이므로 6.9정도 넣고 돌림)
- 칠레 마줄레주해역 지진(8.8)
- 2011년
- 진도 9 이상 일본 한 건, 진도가 9이면서 residual이 1 초과인 1곳 존재
- 2016년
- 경주 지진 있던 해, 진도가 5.4 이상, 해당 데이터가 존재하지 않아 지도에도 표시되지 않음(한국의 지진 데이터가 없는 것 같다.)
class MooYaHo3(MooYaHo2):
def vis(self,MagThresh=7,ResThresh=1):
fig = px.density_mapbox(self.df,
lat='Latitude',
lon='Longitude',
z='Magnitude',
radius=15,
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})))
- make instance for analysis
moo_korea=MooYaHo3(df_korea.query("Year >= 2010"))
- get distance
moo_korea.get_distance()
moo_korea.D[moo_korea.D>0].mean()
plt.hist(moo_korea.D[moo_korea.D>0])
- weight matrix
moo_korea.get_weightmatrix(theta=(280.08),kappa=200)
- fit
moo_korea.fit2()
moo_korea.df.query('Magnitude>5')
plt.plot((moo_korea.df.Residual)**2)
- vis
- make instance for analysis
moo_global=MooYaHo(df_global.query('Year>=2010'))
- 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=(8268.6226),kappa=4000)
- fit
moo_global.fit(m=10000)
- vis
moo_global.df.query('Year==2010 & Latitude == 18.443')