Sekilas Dasar¶
Penentuan episenter metode lingkaran memanfaatkan minimal 3 stasiun seismik. Data yang dibutuhkan dari setiap stasiun adalah waktu tiba gelombang P dan gelombang S yang kemudian diselisihkan untuk mendapatkan waktu S-P. Selisih waktu S-P ini dipakai untuk memperkirakan jarak yang sudah ditempuh yang pada contoh ini akan mengaplikasikan metode Tau P. Metode ini sering juga disebut circle and chord karena memanfaatkan perpotongan lingkaran dan garis perpotongannya.
Dalam gambar (Havskov et al, 2002) di atas S merupakan stasiun seismik dan d merupakan jarak yang didapatkan dari pencocokan nilai S-P dengan jarak.
Data¶
Data yang akan kita gunakan adalah data waktu tiba gelombang P dan gelombang S dari gempa Mamuju M6,2 yang direkam di stasiun TATO: Taiwan, PMG: Port Moresby, dan COCO: Cocos Island. Data waktu tiba ini dikonversi kedalam format UTCDateTime obspy:
In [5]:
from obspy import UTCDateTime
TATO_P = UTCDateTime("2021-01-14T18:34:08")
TATO_S = UTCDateTime("2021-01-14T18:38:53")
PMG_P = UTCDateTime("2021-01-14T18:34:15")
PMG_S = UTCDateTime("2021-01-14T18:39:04")
COCO_P = UTCDateTime("2021-01-14T18:33:28")
COCO_S = UTCDateTime("2021-01-14T18:37:43")
Selisih waktu tiba S dan P¶
Selisih S-P kemudian dapat kita cari dengan cara yang sederhana yaitu:
In [6]:
TATO_SP = TATO_S-TATO_P
PMG_SP = PMG_S-PMG_P
COCO_SP = COCO_S-COCO_P
print(COCO_SP,"sekon")
255.0 sekon
Mencari jarak dengan rule of thumb¶
Dari nilai SP sebenarnya kita langsung dapat mencari radius lingkaran/jarak dengan rumus Havskov et al (2002):
$$jarak=(SP-2)\times10,$$
dengan satuan kilometer, apabila akan dijadikan derajat maka kita bagi dengan angka 111,26 km:
In [7]:
jarak=[(SP-2)*10/111.26 for SP in [TATO_SP,PMG_SP,COCO_SP]]
jarak
Out[7]:
[25.43591587273054, 25.795434118281502, 22.73952903109833]
Kemudian dapat kita plot:
In [8]:
import matplotlib.pyplot as plt
# koordinat stasiun
TATO_coor=(121.50,24.97)
PMG_coor=(147.16,-9.40)
COCO_coor=(96.83,-12.19)
nama_stasiun = ["TATO","PMG","COCO"]
figure, ax = plt.subplots()
ax.set_aspect(1)
for sta,rad,nsta in zip([TATO_coor, PMG_coor, COCO_coor],\
jarak,nama_stasiun):
ax.scatter(sta[0],sta[1], zorder=100, marker="v", color="blue")
draw_circle =plt.Circle(sta, rad, fill=False)
ax.add_artist(draw_circle)
ax.annotate(nsta, sta)
ax.set_title('Episenter Metode Lingkaran Gempa Mamuju M6,2')
ax.set_xlim(70,180)
ax.set_ylim(-40,60)
ax.scatter(118.8893,-2.9711,color="red", label="Episenter sebenarnya (IRIS)")
ax.legend()
plt.show()
Hasil pendekatan ini cukup tetapi belum akurat, kita bisa menambah akurasi dengan memanfaatkan Tau P.
Mencari jarak berdasarkan S-P¶
Pengenalan modul TauPyModel dari obspy¶
Sebelum mencari jarak berdasarkan S-P kita akan berkenalan dahulu dengan modul TauPyModel dari obspy. Modul ini memungkinkan kita mencari nilai waktu tempuh dari gempa ke stasiun dengan lokasi dan kedalaman tertentu. Modul ini berdasarkan algoritma Tau P dan bisa kita aplikasikan dengan berbagai model kecepatan yang sudah ada, kita akan memakai model IASP91. Selain menghitung waktu tempuh, modul juga dapat mengeplot jejak sinar/raypath:
In [9]:
from obspy.taup import TauPyModel
model = TauPyModel(model="iasp91")
#contoh menghitung waktu tempuh gempa kedalaman 60 km dengan jarak episenter-stasiun adalah 60 derajat
arrivals = model.get_travel_times(source_depth_in_km=60, #kedalaman dlm km
distance_in_degree=60)# jarak dalam derajat
print(arrivals)
26 arrivals P phase arrival at 600.258 seconds pP phase arrival at 616.298 seconds sP phase arrival at 623.358 seconds PcP phase arrival at 645.616 seconds PP phase arrival at 733.120 seconds ScP phase arrival at 881.182 seconds PcS phase arrival at 888.006 seconds PKiKP phase arrival at 1024.399 seconds pPKiKP phase arrival at 1042.077 seconds sPKiKP phase arrival at 1048.763 seconds S phase arrival at 1088.878 seconds SP phase arrival at 1100.115 seconds pS phase arrival at 1107.826 seconds PS phase arrival at 1109.778 seconds sS phase arrival at 1116.576 seconds ScS phase arrival at 1185.125 seconds SKiKP phase arrival at 1231.805 seconds SS phase arrival at 1327.584 seconds PKIKKIKP phase arrival at 1864.119 seconds SKIKKIKP phase arrival at 2071.513 seconds PKIKKIKS phase arrival at 2078.199 seconds SKIKKIKS phase arrival at 2285.426 seconds PKIKPPKIKP phase arrival at 2364.644 seconds PKPPKP phase arrival at 2375.132 seconds PKPPKP phase arrival at 2386.588 seconds SKIKSSKIKS phase arrival at 3216.623 seconds
untuk mengeplot kita dapat menggunakan get_raypaths:
In [10]:
raypaths = model.get_ray_paths(source_depth_in_km=60, #kedalaman dlm km
distance_in_degree=60)# jarak dalam derajat
ax = raypaths.plot_rays()
kita juga dapat membatasi fase apa saja yang akan kita ambil dengan parameter phase_list:
In [11]:
arrivals = model.get_travel_times(source_depth_in_km=60, #kedalaman dlm km
distance_in_degree=60,# jarak dalam derajat
phase_list=["P","S"])#daftar fase yang akan diambil
print(arrivals)
2 arrivals P phase arrival at 600.258 seconds S phase arrival at 1088.878 seconds
In [12]:
raypaths = model.get_ray_paths(source_depth_in_km=60, #kedalaman dlm km
distance_in_degree=60, # jarak dalam derajat
phase_list=["P","S"])#daftar fase yang akan diambil
ax = raypaths.plot_rays()
untuk mendapatkan angka waktu tibanya kita menggunakan indeksing dan atribut .time dan .name, contohnya untuk mendapatkan waktu tiba pertama maka:
In [13]:
first_arrival=arrivals[0]
nama_fase=first_arrival.name
waktu_tiba=first_arrival.time
print("Fase pertama datang {} pada {} sekon".format(nama_fase, waktu_tiba))
Fase pertama datang P pada 600.2579575887129 sekon
Untuk waktu tiba P jelas pertama sehingga kita menggunakan indeks 0 seperti contoh di atas, hanya untuk waktu tiba gelombang S sedikit perlu modifikasi karena bisa saja ada banyak fase P dan S dalam 1 kali kejadian, contoh:
In [14]:
arrivals = model.get_travel_times(source_depth_in_km=10,
distance_in_degree=23.6,
phase_list=["P","S"])
print(arrivals)
8 arrivals P phase arrival at 311.128 seconds P phase arrival at 311.211 seconds P phase arrival at 312.946 seconds S phase arrival at 566.326 seconds S phase arrival at 568.932 seconds S phase arrival at 571.185 seconds S phase arrival at 586.578 seconds S phase arrival at 586.600 seconds
Pada contoh di atas ada 3 kedatangan gelombang P adn 5 kedatangan gelombang S, bagaimana cara mengambil gelombang S yang pertama? Kita bisa sedikit membuat algoritma:
Mengambil nama-nama fasenya dengan loop:
In [15]:
import numpy as np
phase_names=np.array([arr.name for arr in arrivals])
phase_names
Out[15]:
array(['P', 'P', 'P', 'S', 'S', 'S', 'S', 'S'],
dtype='<U1') Mencari dimana letak indeks S:
In [16]:
S_indexes=np.where(phase_names=="S")
S_indexes
Out[16]:
(array([3, 4, 5, 6, 7]),)
Mengambil indeks S paling pertama sebagai indeks waktu tiba S pertama kali:
In [17]:
S_index=S_indexes[0][0]
S_index
Out[17]:
3
Memakai indeks tersebut untuk mengambil S pertama:
In [18]:
Sfirst_arrival=arrivals[S_index]
Snama_fase=Sfirst_arrival.name
Swaktu_tiba=Sfirst_arrival.time
print("Fase S pertama datang {} pada {} sekon".format(Snama_fase, Swaktu_tiba))
Fase S pertama datang S pada 566.3255616737924 sekon
Memanfaatkan TauPyModel untuk mencari jarak¶
Seperti yang sudah kita bahas di atas, TauPyModel berfungsi untuk mencari waktu tempuh gelombang. Dalam algoritma metode lingkaran kita harus melakukan sebaliknya yaitu dari waktu tempuh menjadi jarak. Oleh karena itu kita akan membuat kode singkat untuk membalik algoritmanya, dalam kasus ini kedalaman kita anggap tetap yaitu 10 km.
Langkah selanjutnya adalah kita mencoba masing-masing jarak dengan for loop dan menghitung nilai S-P dari TauPyModel: SP_taup=arrivals[-1].time-arrivals[0].time. Di bawah langkah ini ada kondisional yang menyatakan bila waktu S-P dari TauPyModel lebih kecil dari S-P stasiun yang bersangkutan maka nilai tersebut akan dipakai sebagai jarak, tetapi apabila nilai S-P dari TauPyModel lebih besar maka nilai jarak tidak diupdate.
Sebagai bayangan jika kita akan menebak suatu jarak dengan waktu S-P adalah 250 sekon dan kedalaman gempa adalah 10 km dengan jangkauan tebakan kita adalah 22 sampai 25 dengan kelipatan 1:
In [20]:
distance_list=np.arange(22,25,1) #jangkauan tebakan
SP=250 #SP (s)
depth=10 #km
initial_distance=distance_list[0]
for dist in distance_list:
arrivals = model.get_travel_times(source_depth_in_km=depth,
distance_in_degree=dist,
phase_list=["P","S"])
# waktu tiba P pasti paling pertama
P_arrival = arrivals[0].time
# waktu tiba S harus sedikit hati-hati karena bisa ada banyak P dan S
phase_names=np.array([arr.name for arr in arrivals]) # mengambil nama-nama fase
S_indexes=np.where(phase_names=="S") #mengambil indeks lokasi S
S_index=S_indexes[0][0] # mengambil indeks lokasi S pertama kali
S_arrival = arrivals[S_index].time # waktu tiba S
# menghitung S-P dari TauP IASP91
SP_taup = S_arrival-P_arrival
if SP_taup<SP: #jika SP_taup < SP observasi maka jarak diupdate
print("SP taup < SP: {} < {}".format(SP_taup,SP))
print("Initial distance diupdate: {}".format(initial_distance))
initial_distance = dist
else: # jika tidak maka jarak tidak diupdate
print("SP taup > SP: {} > {}".format(SP_taup,SP))
print("Initial distance tidak diupdate: {}".format(initial_distance))
print("Jarak yang terpilih: {}".format(initial_distance))
SP taup < SP: 243.80686273283436 < 250 Initial distance diupdate: 22 SP taup > SP: 251.6500605586528 > 250 Initial distance tidak diupdate: 22 SP taup > SP: 258.02310930613214 > 250 Initial distance tidak diupdate: 22 Jarak yang terpilih: 22
Contoh di atas memperlihatkan bahwa nilai jarak yang terpilih adalah 22 karena saat dicoba nilai SPnya 243, lebih kecil dari 250 atau SP pengamatan. Saat jarak 23 terlihat bahwa nilai sudah lebih besar dari 250 yaitu 251 sehingga jarak 23 seterusnya tidak dipakai. Melihat algoritma di atas masih memiliki kelemahan yaitu bahwa walaupun pada jarak 23 derajat nilai SP>250 tetapi selisihnya hanya sekitar 1,65 detik, jauh lebih kecil dari selisih pada jarak 22 derajat yaitu sekitar 6,2 detik, oleh karena itu kita memperkenalkan kondisional baru:
In [21]:
distance_list=np.arange(22,25,1)
SP=250
depth=10 #km
initial_distance=distance_list[0]
for dist in distance_list:
arrivals = model.get_travel_times(source_depth_in_km=depth,
distance_in_degree=dist,
phase_list=["P","S"])
P_arrival = arrivals[0].time
# waktu tiba S harus sedikit hati-hari karena bisa ada banyak P dan S
phase_names=np.array([arr.name for arr in arrivals]) # mengambil nama-nama fase
S_indexes=np.where(phase_names=="S") #mengambil indeks lokasi S
S_index=S_indexes[0][0] # mengambil indeks lokasi S pertama kali
S_arrival = arrivals[S_index].time # waktu tiba S
SP_taup = S_arrival-P_arrival
if SP_taup<SP:
initial_distance = dist
print("SP taup < SP: {} < {}".format(SP_taup,SP))
print("Initial distance diupdate: {}".format(initial_distance))
SP_error=abs(SP_taup-SP) #menambah perhitungan SP_error
else:
SP_error_check = abs(SP_taup-SP)
if SP_error_check < SP_error: # jika SP_error sekarang < SP_error sebelumnya maka jarak tetap diupdate
print("SP taup > SP: {} > {}".format(SP_taup,SP))
print("Tetapi SP_error < SP_error_sebelumnya: {} < {}".format(SP_error_check, SP_error))
initial_distance = dist
print("Initial distance diupdate: {}".format(initial_distance))
else:
print("SP taup > SP: {} > {}".format(SP_taup,SP))
print("Initial distance tidak diupdate: {}".format(initial_distance))
print("Jarak yang terpilih: {}".format(initial_distance))
SP taup < SP: 243.80686273283436 < 250 Initial distance diupdate: 22 SP taup > SP: 251.6500605586528 > 250 Tetapi SP_error < SP_error_sebelumnya: 1.6500605586527968 < 6.193137267165639 Initial distance diupdate: 23 SP taup > SP: 258.02310930613214 > 250 Initial distance tidak diupdate: 23 Jarak yang terpilih: 23
Kita sudah mendapatkan jarak yang nilai SPnya lebih mendekati yaitu 23 derajat.
Menentukan range tebakan jarak¶
Kita bisa mencocokkan SP dengan kurva waktu tempuh, tetapi SP kita harus konversi dulu dalam menit:
In [27]:
print(TATO_SP/60, PMG_SP/60, COCO_SP/60)
4.75 4.816666666666666 4.25
dengan menggunakan Tau P kita bisa membuat kurva waktu tempuh:
In [43]:
from obspy.taup import plot_travel_times
import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(6, 6))
ax = plot_travel_times(source_depth=10, phase_list=["P", "S"],
ax=ax, fig=fig, verbose=True,
show=False)
ax.set_title("Kurva Waktu Tempuh P S untuk Kedalaman 10 km")
ax.grid()
There were 23 epicentral distances without an arrival
Karena S-P data kita dari 4,25 sampai 4,8, dengan melihat kurva di atas kita bisa zoom di sekitar 10 sampai 40 derajat:
In [41]:
fig, ax = plt.subplots(figsize=(6, 6))
ax = plot_travel_times(source_depth=10, phase_list=["P", "S"],
ax=ax, fig=fig, verbose=True, min_degrees=10, max_degrees=40, #zoom
show=False)
ax.set_title("Kurva Waktu Tempuh P S untuk Kedalaman 10 km")
ax.grid()
There were 0 epicentral distances without an arrival
dari kurva zoom di atas kita kemudian memilih range tebakan dari 20 sampai 30 derajat.
Implementasi untuk mencari jarak stasiun-gempa yang benar¶
Untuk mencari jarak yang tepat di setiap stasiun kita menambahkan satu for loop lagi dan jangkauan tebakan dari 20 sampai 30 derajat dengan lompatan yang lebih halus yaitu 0,1:
In [22]:
distance_list=np.arange(20,30,0.1)
station_SP_list=[TATO_SP,PMG_SP,COCO_SP]
sta_distance_list=[]
for sta_SP in station_SP_list:
initial_distance=distance_list[0]
for dist in distance_list:
arrivals = model.get_travel_times(source_depth_in_km=depth,
distance_in_degree=dist,
phase_list=["P","S"])
P_arrival = arrivals[0].time
# waktu tiba S harus sedikit hati-hari karena bisa ada banyak P dan S
phase_names=np.array([arr.name for arr in arrivals]) # mengambil nama-nama fase
S_indexes=np.where(phase_names=="S") #mengambil indeks lokasi S
S_index=S_indexes[0][0] # mengambil indeks lokasi S pertama kali
S_arrival = arrivals[S_index].time # waktu tiba S
SP_taup = S_arrival-P_arrival
if SP_taup<sta_SP:
initial_distance = dist
SP_error=abs(SP_taup-sta_SP) #menambah perhitungan SP_error
else:
SP_error_check = abs(SP_taup-sta_SP)
if SP_error_check < SP_error: # jika SP_error sekarang < SP_error sebelumnya maka jarak tetap diupdate
initial_distance = dist
sta_distance_list.append(initial_distance)
sta_distance_list
Out[22]:
[27.900000000000112, 28.500000000000121, 23.50000000000005]
dari algoritma di atas kita mendapatkan nilai jarak masing-masing untuk TATO,PMG, dan COCO adalah 27,9,28,5, dan 23,5.
Mengeplot hasil dengan matplotlib¶
Lingkaran dapat kita plot menggunakan fungsi plt.Circle dengan titik tengah berupa koordinat stasiun dan jari-jari berupa jarak yang sebelumnya kita dapatkan:
In [23]:
import matplotlib.pyplot as plt
# koordinat stasiun
TATO_coor=(121.50,24.97)
PMG_coor=(147.16,-9.40)
COCO_coor=(96.83,-12.19)
#list nama stasiun
nama_stasiun = ["TATO","PMG","COCO"]
figure, ax = plt.subplots()
ax.set_aspect(1) #agar proporsi x dan y sama
# plot lingkaran dan titik
for sta,rad,nsta in zip([TATO_coor, PMG_coor, COCO_coor],\
sta_distance_list,nama_stasiun):
ax.scatter(sta[0],sta[1], zorder=100, marker="v", color="blue") #titik stasiun
draw_circle =plt.Circle(sta, rad, fill=False) #lingkaran
ax.add_artist(draw_circle) #lingkatan
ax.annotate(nsta, sta) # nama stasiun
ax.set_title('Episenter Metode Lingkaran Gempa Mamuju M6,2')
ax.set_xlim(70,180)
ax.set_ylim(-40,60)
ax.scatter(118.8893,-2.9711,color="red", label="Episenter sebenarnya (IRIS)")
ax.legend()
plt.show()
Mengeplot hasil dengan cartopy¶
In [64]:
import numpy as np
import cartopy.crs as ccrs
from cartopy.io import srtm
import matplotlib.pyplot as plt
import cartopy.io.img_tiles as cimgt
fig = plt.figure()
ax = plt.axes(projection=ccrs.PlateCarree())
ax.set_extent((70,180,-40,60))
# Add stamen background
stamen_terrain = cimgt.Stamen('terrain-background') ##maps.stamen.com
ax.add_image(stamen_terrain,2)
ax.set_xlabel("Longitude [deg]")
ax.set_ylabel("Latitude [deg]")
gl = ax.gridlines(draw_labels=True, dms=False, x_inline=False, y_inline=False, zorder=11)
gl.xlabels_top = False
gl.ylabels_right = False
for sta,rad,nsta in zip([TATO_coor, PMG_coor, COCO_coor],\
sta_distance_list,nama_stasiun):
ax.scatter(sta[0],sta[1], zorder=100, marker="v", color="blue")
draw_circle =plt.Circle(sta, rad, fill=False)
ax.add_artist(draw_circle)
ax.annotate(nsta, sta)
ax.scatter(118.8893,-2.9711,color="red", label="Episenter sebenarnya (IRIS)")
plt.savefig("output/demo_epi.png",dpi=300)
/home/anangsahroni/Environments/envgeneral/lib/python3.9/site-packages/cartopy/mpl/gridliner.py:324: UserWarning: The .xlabels_top attribute is deprecated. Please use .top_labels to toggle visibility instead.
warnings.warn('The .xlabels_top attribute is deprecated. Please '
/home/anangsahroni/Environments/envgeneral/lib/python3.9/site-packages/cartopy/mpl/gridliner.py:360: UserWarning: The .ylabels_right attribute is deprecated. Please use .right_labels to toggle visibility instead.
warnings.warn('The .ylabels_right attribute is deprecated. Please '
Rekap kode¶
Jika kode-kode di atas direkap dalam 1 bagian maka menjadi:
In [5]:
import numpy as np
from obspy import UTCDateTime
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import cartopy.io.img_tiles as cimgt
from obspy.taup import TauPyModel
# MODEL KECEPATAN
model = TauPyModel(model="iasp91")
# KEDALAMAN
depth=10 #km
# DATA
TATO_P = UTCDateTime("2021-01-14T18:34:08")
TATO_S = UTCDateTime("2021-01-14T18:38:53")
PMG_P = UTCDateTime("2021-01-14T18:34:15")
PMG_S = UTCDateTime("2021-01-14T18:39:04")
COCO_P = UTCDateTime("2021-01-14T18:33:28")
COCO_S = UTCDateTime("2021-01-14T18:37:43")
# KOORDINAT dan NAMA
TATO_coor=(121.50,24.97)
PMG_coor=(147.16,-9.40)
COCO_coor=(96.83,-12.19)
nama_stasiun = ["TATO","PMG","COCO"]
# SP
TATO_SP = TATO_S-TATO_P
PMG_SP = PMG_S-PMG_P
COCO_SP = COCO_S-COCO_P
# MENCARI JARAK
distance_list=np.arange(20,30,0.1) #Range tebakan
station_SP_list=[TATO_SP,PMG_SP,COCO_SP]
sta_distance_list=[]
for sta_SP in station_SP_list:
initial_distance=distance_list[0]
for dist in distance_list:
arrivals = model.get_travel_times(source_depth_in_km=depth,
distance_in_degree=dist,
phase_list=["P","S"])
P_arrival = arrivals[0].time
# waktu tiba S harus sedikit hati-hari karena bisa ada banyak P dan S
phase_names=np.array([arr.name for arr in arrivals]) # mengambil nama-nama fase
S_indexes=np.where(phase_names=="S") #mengambil indeks lokasi S
S_index=S_indexes[0][0] # mengambil indeks lokasi S pertama kali
S_arrival = arrivals[S_index].time # waktu tiba S
SP_taup = S_arrival-P_arrival
if SP_taup<sta_SP:
initial_distance = dist
SP_error=abs(SP_taup-sta_SP) #menambah perhitungan SP_error
else:
SP_error_check = abs(SP_taup-sta_SP)
if SP_error_check < SP_error: # jika SP_error sekarang < SP_error sebelumnya maka jarak tetap diupdate
initial_distance = dist
sta_distance_list.append(initial_distance)
# PLOT
fig = plt.figure()
ax = plt.axes(projection=ccrs.PlateCarree())
ax.set_extent((70,180,-40,60))
# Add stamen background
stamen_terrain = cimgt.Stamen('terrain-background') ##maps.stamen.com
ax.add_image(stamen_terrain,2) # 2 merupakan level zoom, semakin besar semakin detail tetapi semakin lama proses download, gunakan dengan bijak
ax.set_xlabel("Longitude [deg]")
ax.set_ylabel("Latitude [deg]")
gl = ax.gridlines(draw_labels=True, dms=False, x_inline=False, y_inline=False, zorder=11)
gl.xlabels_top = False
gl.ylabels_right = False
for sta,rad,nsta in zip([TATO_coor, PMG_coor, COCO_coor],\
sta_distance_list,nama_stasiun):
ax.scatter(sta[0],sta[1], zorder=100, marker="v", color="blue")
draw_circle =plt.Circle(sta, rad, fill=False)
ax.add_artist(draw_circle)
ax.annotate(nsta, sta)
ax.scatter(118.8893,-2.9711,color="red", label="Episenter sebenarnya (IRIS)")
plt.savefig("output/demo_epi.png",dpi=300)
/home/anangsahroni/Environments/envgeneral/lib/python3.9/site-packages/cartopy/mpl/gridliner.py:324: UserWarning: The .xlabels_top attribute is deprecated. Please use .top_labels to toggle visibility instead.
warnings.warn('The .xlabels_top attribute is deprecated. Please '
/home/anangsahroni/Environments/envgeneral/lib/python3.9/site-packages/cartopy/mpl/gridliner.py:360: UserWarning: The .ylabels_right attribute is deprecated. Please use .right_labels to toggle visibility instead.
warnings.warn('The .ylabels_right attribute is deprecated. Please '
Tambahan: Mendapatkan titik episenter¶
Dengan menyederhanakan menjadi bidang datar kita dapat memanfaatkan pengetahuan geometri untuk mendapatkan dua titik perpotongan dua lingkaran:
In [24]:
# referensi: https://stackoverflow.com/questions/64800285/intersection-coordinates-lat-lon-of-two-circles-given-the-coordinates-of-the
import math
def circle_intersection(x0, y0, r0, x1, y1, r1):
d = math.sqrt((x1 - x0) ** 2 + (y1 - y0) ** 2)
if d > r0 + r1: # non intersecting
return None
if d < abs(r0 - r1): # one circle within other
return None
if d == 0 and r0 == r1: # coincident circles
return None
a = (r0 ** 2 - r1 ** 2 + d ** 2) / (2 * d)
h = math.sqrt(r0 ** 2 - a ** 2)
x2 = x0 + a * (x1 - x0) / d
y2 = y0 + a * (y1 - y0) / d
x3 = x2 + h * (y1 - y0) / d
y3 = y2 - h * (x1 - x0) / d
x4 = x2 - h * (y1 - y0) / d
y4 = y2 + h * (x1 - x0) / d
return (x3, y3), (x4, y4)
kita gunakan fungsi di atas untuk menghitung perpotongan masing-masing pasangan lingkaran stasiun:
In [54]:
TP_1, TP_2 = circle_intersection(TATO_coor[0], TATO_coor[1], sta_distance_list[0],\
PMG_coor[0], PMG_coor[1], sta_distance_list[1])
TC_1, TC_2 = circle_intersection(TATO_coor[0], TATO_coor[1], sta_distance_list[0],\
COCO_coor[0], COCO_coor[1], sta_distance_list[2])
PC_1, PC_2 = circle_intersection(PMG_coor[0], PMG_coor[1], sta_distance_list[1],\
COCO_coor[0], COCO_coor[1], sta_distance_list[2])
print(TP_1, TP_2)
(119.42228880002661, -2.852528931955689) (148.76572171237126, 19.054729972052563)
kemudian kita ubah bentuk ke dalam list untuk mempermudah:
In [55]:
#TATO PMG
TPx=[TP_1[0], TP_2[0]]
TPy=[TP_1[1], TP_2[1]]
#TATO COCO
TCx=[TC_1[0], TC_2[0]]
TCy=[TC_1[1], TC_2[1]]
# PMG_COCO
PCx=[PC_1[0], PC_2[0]]
PCy=[PC_1[1], PC_2[1]]
dengan masing-masing dua titik perpotongan kita kemudian dapat mengplot garis chord:
In [56]:
import matplotlib.pyplot as plt
# koordinat stasiun
TATO_coor=(121.50,24.97)
PMG_coor=(147.16,-9.40)
COCO_coor=(96.83,-12.19)
nama_stasiun = ["TATO","PMG","COCO"]
figure, ax = plt.subplots()
ax.set_aspect(1)
for sta,rad,nsta in zip([TATO_coor, PMG_coor, COCO_coor],\
sta_distance_list,nama_stasiun):
ax.scatter(sta[0],sta[1], zorder=100, marker="v", color="blue")
draw_circle =plt.Circle(sta, rad, fill=False)
ax.add_artist(draw_circle)
ax.annotate(nsta, sta)
# plot garis chord
ax.plot(TPx, TPy)
ax.plot(TCx, TCy)
ax.plot(PCx, PCy)
ax.set_title('Episenter Metode Lingkaran Gempa Mamuju M6,2')
ax.set_xlim(70,180)
ax.set_ylim(-40,60)
ax.scatter(118.8893,-2.9711,color="red", label="Episenter sebenarnya (IRIS)")
ax.legend()
plt.show()
Memanfaatkan fungsi linregress dari scipy kita dengan mudah dapat mendapatkan nilai slope dan intercept masing-masing garis chord:
In [57]:
from scipy import stats
result = stats.linregress(TPx, TPy)
mTP=result.slope
cTP=result.intercept
result = stats.linregress(TCx, TCy)
mTC=result.slope
cTC=result.intercept
result = stats.linregress(PCx, PCy)
mPC=result.slope
cPC=result.intercept
menggunakan persamaan sederhana untuk mendapatkan perpotongan dua garis:
In [39]:
def two_intersection(m1,c1,m2,c2):
x=(c2-c1)/(m1-m2)
y=m1*x+c1
return x,y
perpotongan setiap pasangan akan dikalkulasi kemudian nantinya direrata:
In [58]:
inter_TPPC=two_intersection(mTP, cTP, mPC, cPC)
inter_TPTC=two_intersection(mTP, cTP, mTC, cTC)
inter_PCTC=two_intersection(mPC, cPC, mTC, cTC)
print(inter_TPPC)
print(inter_TPTC)
print(inter_PCTC)
(118.98966796026603, -3.175515569961405) (118.98966796026605, -3.1755155699613908) (118.98966796026603, -3.175515569961135)
dari titik potong setiap pasang ternyata perbedaannya kecil sampai 12 angka di belakang koma, pererataan dilakukan dengan cara sederhana:
In [62]:
x_av=(inter_TPPC[0]+inter_TPTC[0]+inter_PCTC[0])/3
y_av=(inter_TPPC[1]+inter_TPTC[1]+inter_PCTC[1])/3
print("Episenter hasil kalkulasi",x_av, y_av)
print("Episenter sebenarnya: 118.8893,-2.9711")
Episenter hasil kalkulasi 118.98966796 -3.17551556996 Episenter sebenarnya: 118.8893,-2.9711
Kita dapat plotkan dan akan di zoom agar lebih terlihat:
In [53]:
import matplotlib.pyplot as plt
# koordinat stasiun
TATO_coor=(121.50,24.97)
PMG_coor=(147.16,-9.40)
COCO_coor=(96.83,-12.19)
nama_stasiun = ["TATO","PMG","COCO"]
figure, ax = plt.subplots()
ax.set_aspect(1)
for sta,rad,nsta in zip([TATO_coor, PMG_coor, COCO_coor],\
sta_distance_list,nama_stasiun):
ax.scatter(sta[0],sta[1], zorder=100, marker="v", color="blue")
draw_circle =plt.Circle(sta, rad, fill=False)
ax.add_artist(draw_circle)
ax.annotate(nsta, sta)
ax.plot(TPx, TPy)
ax.plot(TCx, TCy)
ax.plot(PCx, PCy)
ax.set_title('Episenter Metode Lingkaran Gempa Mamuju M6,2')
ax.set_xlim(115,125)
ax.set_ylim(-5,2)
ax.scatter(118.8893,-2.9711,color="red", label="Episenter sebenarnya (IRIS)")
ax.scatter(x_av,y_av,color="green", label="Episenter kalkulasi")
ax.legend()
plt.show()
In [91]:
fig = plt.figure()
ax = plt.axes(projection=ccrs.PlateCarree())
ax.set_extent((115,125,-6,1))
# Add custom background
terrain = cimgt.GoogleTiles(url='https://server.arcgisonline.com/ArcGIS/rest/services/World_Shaded_Relief/MapServer/tile/{z}/{y}/{x}.jpg',\
cache=True)
ax.add_image(terrain,8) # 8 merupakan level zoom dari basemap,
#semakin besar semakin detail tetapi semakin lama proses download, gunakan dengan bijak
ax.set_xlabel("Longitude [deg]")
ax.set_ylabel("Latitude [deg]")
gl = ax.gridlines(draw_labels=True, dms=False, x_inline=False, y_inline=False, zorder=1)
gl.xlabels_top = False
gl.ylabels_right = False
for sta,rad,nsta in zip([TATO_coor, PMG_coor, COCO_coor],\
sta_distance_list,nama_stasiun):
ax.scatter(sta[0],sta[1], zorder=100, marker="v", color="blue")
draw_circle =plt.Circle(sta, rad, fill=False)
ax.add_artist(draw_circle)
ax.annotate(nsta, sta)
ax.plot(TPx, TPy)
ax.plot(TCx, TCy)
ax.plot(PCx, PCy)
ax.scatter(118.8893,-2.9711,color="red", label="Episenter sebenarnya (IRIS)")
ax.scatter(x_av,y_av,color="green", label="Episenter kalkulasi")
ax.legend()
plt.savefig("output/demo_epi_pt.png",dpi=300)
Rekap tambahan¶
Jika kode rekap ditambah algoritma untuk mencari titik:
In [6]:
import math
import numpy as np
from scipy import stats
from obspy import UTCDateTime
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import cartopy.io.img_tiles as cimgt
from obspy.taup import TauPyModel
# MODEL KECEPATAN
model = TauPyModel(model="iasp91")
# KEDALAMAN
depth=10 #km
# DATA
TATO_P = UTCDateTime("2021-01-14T18:34:08")
TATO_S = UTCDateTime("2021-01-14T18:38:53")
PMG_P = UTCDateTime("2021-01-14T18:34:15")
PMG_S = UTCDateTime("2021-01-14T18:39:04")
COCO_P = UTCDateTime("2021-01-14T18:33:28")
COCO_S = UTCDateTime("2021-01-14T18:37:43")
# KOORDINAT dan NAMA
TATO_coor=(121.50,24.97)
PMG_coor=(147.16,-9.40)
COCO_coor=(96.83,-12.19)
nama_stasiun = ["TATO","PMG","COCO"]
# SP
TATO_SP = TATO_S-TATO_P
PMG_SP = PMG_S-PMG_P
COCO_SP = COCO_S-COCO_P
# MENCARI JARAK
distance_list=np.arange(20,30,0.1) #Range tebakan
station_SP_list=[TATO_SP,PMG_SP,COCO_SP]
sta_distance_list=[]
for sta_SP in station_SP_list:
initial_distance=distance_list[0]
for dist in distance_list:
arrivals = model.get_travel_times(source_depth_in_km=depth,
distance_in_degree=dist,
phase_list=["P","S"])
P_arrival = arrivals[0].time
# waktu tiba S harus sedikit hati-hari karena bisa ada banyak P dan S
phase_names=np.array([arr.name for arr in arrivals]) # mengambil nama-nama fase
S_indexes=np.where(phase_names=="S") #mengambil indeks lokasi S
S_index=S_indexes[0][0] # mengambil indeks lokasi S pertama kali
S_arrival = arrivals[S_index].time # waktu tiba S
SP_taup = S_arrival-P_arrival
if SP_taup<sta_SP:
initial_distance = dist
SP_error=abs(SP_taup-sta_SP) #menambah perhitungan SP_error
else:
SP_error_check = abs(SP_taup-sta_SP)
if SP_error_check < SP_error: # jika SP_error sekarang < SP_error sebelumnya maka jarak tetap diupdate
initial_distance = dist
sta_distance_list.append(initial_distance)
# MENCARI TITIK
# fungsi titik potong 2 lingkaran
def circle_intersection(x0, y0, r0, x1, y1, r1):
d = math.sqrt((x1 - x0) ** 2 + (y1 - y0) ** 2)
if d > r0 + r1: # non intersecting
return None
if d < abs(r0 - r1): # one circle within other
return None
if d == 0 and r0 == r1: # coincident circles
return None
a = (r0 ** 2 - r1 ** 2 + d ** 2) / (2 * d)
h = math.sqrt(r0 ** 2 - a ** 2)
x2 = x0 + a * (x1 - x0) / d
y2 = y0 + a * (y1 - y0) / d
x3 = x2 + h * (y1 - y0) / d
y3 = y2 - h * (x1 - x0) / d
x4 = x2 - h * (y1 - y0) / d
y4 = y2 + h * (x1 - x0) / d
return (x3, y3), (x4, y4)
# fungsi titik potong dua garis
def two_intersection(m1,c1,m2,c2):
x=(c2-c1)/(m1-m2)
y=m1*x+c1
return x,y
# menghitung titik perpotongan lingkaran
TP_1, TP_2 = circle_intersection(TATO_coor[0], TATO_coor[1], sta_distance_list[0],\
PMG_coor[0], PMG_coor[1], sta_distance_list[1])
TC_1, TC_2 = circle_intersection(TATO_coor[0], TATO_coor[1], sta_distance_list[0],\
COCO_coor[0], COCO_coor[1], sta_distance_list[2])
PC_1, PC_2 = circle_intersection(PMG_coor[0], PMG_coor[1], sta_distance_list[1],\
COCO_coor[0], COCO_coor[1], sta_distance_list[2])
#TATO PMG
TPx=[TP_1[0], TP_2[0]]
TPy=[TP_1[1], TP_2[1]]
#TATO COCO
TCx=[TC_1[0], TC_2[0]]
TCy=[TC_1[1], TC_2[1]]
# PMG_COCO
PCx=[PC_1[0], PC_2[0]]
PCy=[PC_1[1], PC_2[1]]
# mencari slope dan itersep
result = stats.linregress(TPx, TPy)
mTP=result.slope
cTP=result.intercept
result = stats.linregress(TCx, TCy)
mTC=result.slope
cTC=result.intercept
result = stats.linregress(PCx, PCy)
mPC=result.slope
cPC=result.intercept
#mencari perpotongan garis chord
inter_TPPC=two_intersection(mTP, cTP, mPC, cPC)
inter_TPTC=two_intersection(mTP, cTP, mTC, cTC)
inter_PCTC=two_intersection(mPC, cPC, mTC, cTC)
#mererata perpotongan garis
x_av=(inter_TPPC[0]+inter_TPTC[0]+inter_PCTC[0])/3
y_av=(inter_TPPC[1]+inter_TPTC[1]+inter_PCTC[1])/3
# PLOT
fig = plt.figure()
ax = plt.axes(projection=ccrs.PlateCarree())
ax.set_extent((115,125,-6,1))
# Add custom background
terrain = cimgt.GoogleTiles(url='https://server.arcgisonline.com/ArcGIS/rest/services/World_Shaded_Relief/MapServer/tile/{z}/{y}/{x}.jpg',\
cache=True)
ax.add_image(terrain,8)
ax.set_xlabel("Longitude [deg]")
ax.set_ylabel("Latitude [deg]")
gl = ax.gridlines(draw_labels=True, dms=False, x_inline=False, y_inline=False, zorder=1)
gl.xlabels_top = False
gl.ylabels_right = False
for sta,rad,nsta in zip([TATO_coor, PMG_coor, COCO_coor],\
sta_distance_list,nama_stasiun):
ax.scatter(sta[0],sta[1], zorder=100, marker="v", color="blue")
draw_circle =plt.Circle(sta, rad, fill=False)
ax.add_artist(draw_circle)
ax.annotate(nsta, sta)
ax.plot(TPx, TPy)
ax.plot(TCx, TCy)
ax.plot(PCx, PCy)
ax.scatter(118.8893,-2.9711,color="red", label="Episenter sebenarnya (IRIS)")
ax.scatter(x_av,y_av,color="green", label="Episenter kalkulasi")
ax.legend()
plt.savefig("output/demo_epi_pt.png",dpi=300)
Latihan¶
Stasiun:
NET STA LAT LOT ELEV
II ABPO -19.02° 47.23° 1528 m
P 2021-05-14 06:42:16
S 2021-05-14 06:49:41
IC LSA 29.70° 91.13° 3660 m
P 2021-05-14 06:39:15
S 2021-05-14 06:44:14
II WRAB -19.93° 134.36° 366 m
P 2021-05-14 06:40:58
S 2021-05-14 06:47:17
Dari data di atas carilah perkiraan episenter menggunakan Metode Lingkaran!
Bandingkan dengan episenter sebenarnya dari IRIS:
0.1364° N 96.6442° E 2021-05-14 06:33:07 UTC 11.0 km mww6.7 Off W Coast Of Northern Sumatra
Hitung S-P dahulu untuk memperkirakan range penebakan jarak¶
In [48]:
from obspy import UTCDateTime
# masukkan data
# hitung SP
# konversi dalam menit
SP dalam menit 7.416666666666667 4.983333333333333 6.316666666666666
buat kurva waktu tempuh seperti di bagian penentuan range
In [45]:
from obspy.taup import plot_travel_times
import matplotlib.pyplot as plt
# gunakan plot_travel_times untuk membuat kurva waktu tempuh sebagai perkiraan tebakan jarak
There were 23 epicentral distances without an arrival
Out[45]:
Text(0.5, 1.0, 'Kurva Waktu Tempuh P S untuk Kedalaman 10 km')
Dengan mencocokkan SP dalam menit dengan kurva waktu tempuh di atas kita akan menggunakan tebakan jarak dari ... dengan lompatan ....
Solusi tanpa titik episenter¶
Dengan mengubah sedikit kode yang sudah di Rekap kode di atas, yang diubah adalah DATA, KOODINAT dan NAMA, range tebakan (distance list), ax.set_extent, dan episenter sebenarnya (IRIS):
In [46]:
import math
import numpy as np
from scipy import stats
from obspy import UTCDateTime
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import cartopy.io.img_tiles as cimgt
from obspy.taup import TauPyModel
# gunakan kode pada bagian Rekap kode
# Ubah:
#- DATA
#- NAMA dan KOORDINAT
#- distance_list (tentukan dari berapa sampai berapa dan kelipatannya, dari kurva waktu tempuh di atas)
#- koordinat episenter sebenarnya (IRIS)
Solusi dengan titik hiposenter¶
Modifikasi Rekap tambahan
In [12]:
import math
import numpy as np
from scipy import stats
from obspy import UTCDateTime
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import cartopy.io.img_tiles as cimgt
from obspy.taup import TauPyModel
# gunakan kode pada bagian Rekap tambahan
# Ubah:
#- DATA
#- NAMA dan KOORDINAT
#- distance_list (tentukan dari berapa sampai berapa dan kelipatannya, dari kurva waktu tempuh di atas)
#- koordinat episenter sebenarnya (IRIS)
Referensi¶
Havskov J., Bormann P., Schweitzer J., 2002, Earthquake location, New Manual on Seismological Observatory Practice: GFZ Potsdam