#!/usr/bin/env python
# coding: utf-8

# <div style="background-color: #28a745 !important; color: white; min-height: 50px; padding: 10px; margin: auto;">
#     <p style="text-align: right; font-style: italic">STMKGxHMGI Long Course</p>
#     <h1>Sesi: Penentuan Episenter dengan Metode Lingkaran</h1>
# </div>
# <div style="background-color: whitesmoke; padding: 10px ">
#     <ul>
#         <li>Tujuan: Memperkenalkan cara penentuan episenter yang sederhana yaitu menggunakan metode lingkaran</li>
#         <li>Keluaran: Peserta dapat memahami dasar secara sekilas kemudian mengimplementasikan dalam bentuk kode Python</li>
#         <li>Sesi Materi</li>
#         <li>Waktu/Tempat: Sabtu, 25 September 2021/ Zoom</li>
#     </ul>
# </div>

# ## 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.
# 
# ![circle_chord](figures/circle_and_chord.png)
# 
# 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")


# ### 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


# 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)


# 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)


# 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))


# 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)


# 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


# Mencari dimana letak indeks S:

# In[16]:


S_indexes=np.where(phase_names=="S")
S_indexes


# Mengambil indeks S paling pertama sebagai indeks waktu tiba S pertama kali:

# In[17]:


S_index=S_indexes[0][0]
S_index


# 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))


# ### 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))


# 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))


# 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)


# 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()


# 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()


# 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


# 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)


# ### 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)


# ### 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)


# 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)


# 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")


# 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

# Masukkan data dalam `UTCDateTime` ([contoh](#Data)), hitung S-P ([contoh](#Selisih-waktu-tiba-S-dan-P)) dan konversi masing-masing S-P ke dalam menit ([contoh](#Menentukan-range-tebakan-jarak)):

# In[48]:


from obspy import UTCDateTime

# masukkan data

# hitung SP

# konversi dalam menit


# buat kurva waktu tempuh seperti di bagian **[penentuan range](#Menentukan-range-tebakan-jarak)**

# 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


# 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**](#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**](#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