draw_map.py

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

# In[1]:


import pygmt
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings("ignore")
import datetime
from matplotlib.dates import DateFormatter
import matplotlib.dates as mdates
from glob import glob
from math import *


# In[2]:


# focal_res = []
# for a in glob("/home/nghianc/focal_mechanism/PyIsola/MCHAU/output_*/log.txt"):
#     event = a.split("/")[-2].split("_")[-1]
#     data = open(a,"r").readlines()
#     index = [i for i, s in enumerate(data) if '  MT [ Mrr    Mtt    Mpp    Mrt    Mrp    Mtp ]:\n' in s][0]
#     index_o = [i for i, s in enumerate(data) if 'Origin time:' in s][0]
#     index_fp1 = [i for i, s in enumerate(data) if 'Fault plane 1' in s][0]
#     index_fp2 = [i for i, s in enumerate(data) if 'Fault plane 2' in s][0]
    
#     fp1 = data[index_fp1]
#     fp1 = [float(x.replace(",","")) for x in [fp1.split()[x] for x in [5,8,11]]]
#     fp2 = data[index_fp2]
#     fp2 = [float(x.replace(",","")) for x in [fp2.split()[x] for x in [5,8,11]]]
#     index_mw = [i for i, s in enumerate(data) if "Scalar Moment: M0 = " in s][0]
    
#     index_dc = [i for i, s in enumerate(data) if 'DC component:' in s][0]
#     dc = data[index_dc]
#     dc = int(dc.split()[2])
    
#     mw = data[index_mw]
#     mw = float(mw.split()[-1].replace(")",""))
    
#     loc = data[index_o+1].replace("Lat","").replace("Lon","").replace("Depth","").replace("km","").replace("\n","").split()
#     f_tmp = [event,*data[index+1].replace("[","").replace("]","").replace("*","").replace("\n","").split(),*loc,*fp1,*fp2,mw,dc]
#     focal_res.append(f_tmp)


# In[3]:


# focal = pd.DataFrame(focal_res)
# focal.columns = ["id","Mrr","Mtt","Mpp","Mrt","Mrp","Mtp","exp"
#                 ,"lat","lon","dp","str1","dip1","rake1","str2","dip2","rake2","mag"]
# focal.to_csv("../MCHAU_SAC/fm.csv")


# In[4]:


focal_res = pd.read_csv("../MCHAU_SAC/fm.csv")
focal_res.head()


# In[5]:


focalo = pd.DataFrame(focal_res,columns=["event","mrr","mtt","mpp","mrt","mrp","mtp","exp","lat","lon","dep","strike1","dip1","rake1","strike2","dip2","rake2","mw","dc"])
focalo[["mrr","mtt","mpp","mrt","mrp","mtp","exp","lat","lon","dep"]] = focalo[["mrr","mtt","mpp","mrt","mrp","mtp","exp","lat","lon","dep"]].astype('float')
#focal = focal.sort_values(by=["lon","lat"])
#focal.to_csv("../MCHAU_SAC/fm.csv")
focalo


# In[6]:


cities = pd.read_csv("../Map/cities.txt",delim_whitespace=True,names=["Lon","Lat"],index_col=0).reset_index(drop=True)
cities.head()


# In[7]:


sta = pd.read_excel("../Map/station_VN.xlsx",nrows=66).dropna(axis=0,thresh=5).reset_index(drop=True)
sta.head()


# In[8]:


stamsh = ['HBVB', 'MTE2', 'BKVB', 'MLVB', 'SLV', 'LSVB', 'MCVB', 'VTVB', 'TYVB', 'CBVB', 'HGVB', 'VCVB', 'THVB', 'CCVB', 'DBVB']


# In[9]:


# sta = sta.query("Station == {}".format(stamsh))[["Station", "Latitude", "Longtitude"]]


# In[10]:


df = pd.read_csv("../Map/catalog19962020N4M3.txt",sep="\t")
df.head()


# In[11]:


sta_a = list(set(a.split(".")[0] for a in [a.split("_")[2] for a in glob("../MCHAU_SAC/data/202007270512/*HH*")]))


# In[12]:


sta_adf = sta.query("Station == {}".format(stamsh)).reset_index(drop=True)
sta_adf


# In[13]:


# fig = pygmt.Figure()
# pygmt.makecpt(
#     cmap='globe',
#     series='-1000/5000/500',
#     continuous=True
# )
# fig.grdimage(
#     grid=topo_data,
#     region=[minlon, maxlon, minlat, maxlat], 
#     projection='M20c',
#     shading=True,
#     frame="a",
#     )
# fig.coast(
#     region=[minlon, maxlon, minlat, maxlat], 
#     shorelines=True,
#     )
# fig.plot(data="../Map/topo/fault_lv1_edit.txt", pen="1.5,red")
# fig.plot(data="../Map/topo/fault_lv2_edit.txt", pen="1.5")
# fig.show()


# In[16]:


dfc = df.iloc[-17:-16,:].reset_index()
minlon, maxlon = 102, 109
minlat, maxlat = 19, 23.5



# Visualization
fig = pygmt.Figure()
topo_data = '@earth_relief_15s' #30 arc second global relief (SRTM15+V2.1 @ 1.0 km)
pygmt.makecpt(
    cmap='grayC',
    series='-100/10000/200',
    continuous=True
)
fig.grdimage(
    grid=topo_data,
    region=[minlon, maxlon, minlat, maxlat], 
    projection='M20c',
    shading=True,
    frame="a",
    )
fig.coast(
    region=[minlon, maxlon, minlat, maxlat], 
    shorelines=True,
    )

# #Plot earthquakes
dfb = df.iloc[0:-17,:]
fig.plot(
    x=dfb.Lon,
    y=dfb.Lat, 
    size=0.1*1.3** dfb.Ml,
    style='c', 
    color="mediumpurple", 
    pen='1p,black', 
    )

fig.coast(borders=["1/1.0p,black,-"],)
fig.plot(data="../Map/topo/fault_lv1_edit.txt", pen="2.0,red3",label='"Rank 1 fault"')
fig.plot(data="../Map/topo/fault_lv2_edit.txt", pen="1.0",label='"Rank 2 rault"')
fig.legend(position='JBL+jBL+o0.2c',box='+gwhite+p1p')
# fig.plot(
#     x=cities.Lon,
#     y=cities.Lat, 
#     style='d0.3c', 
#     color="white", 
#     pen='red',
#     label="Cities"
#     )
# fig.legend(position='JBL+jBL+o0.2c',box='+gwhite+p1p')



#Plot stations
colorl = ["green","cyan1","deepskyblue"]
# for i,nw in enumerate(sta_adf.Network.unique()):
for i,nw in enumerate(['RM', 'VN', 'TW']):
    sta_nw = sta_adf.query("Network == '{}'".format(nw))
    fig.plot(
    x=sta_nw.Longtitude,
    y=sta_nw.Latitude, 
    style='t0.55c', 
    color=colorl[i], 
    pen='1p,black', 
    label="{}".format(nw),
    )

xc,yc = dfc.Lon.values[0], dfc.Lat.values[0]
l = 0.15

l=0.2
pminlon, pmaxlon = xc-l, xc+l
pminlat, pmaxlat = yc-0.65*l, yc+0.65*l


fig.plot(
    x=[pminlon,pminlon,pmaxlon,pmaxlon,pminlon],
    y=[pminlat,pmaxlat,pmaxlat,pminlat,pminlat],
    pen="2,black",
    )

fig.plot(
    x=dfc.Lon,
    y=dfc.Lat,
    color="red",
    size=0.2*1.3** dfc.Ml,
    style='a', 
    pen='black'
    )

fig.legend()




for m in [3,4,5]:
    mag = 0.1*1.3**m
    fig.plot(
        x=[0],
        y=[0],
        color="mediumpurple",
        style='c{}'.format(mag), 
        pen='black',label="M={}".format(m)
        )

fig.legend(position='JBR+jBR+o0.2c',box='+gwhite+p1p')


fig.text(text="Hanoi", x=106.1, y=21.05, font="10p,Helvetica,black",fill="white")
fig.text(text="LAOS", x=103, y=20, font="20p,Helvetica-Bold,black")
#fig.text(text="VIETNAM", x=105.5, y=21.7, font="20p,Helvetica-Bold,white")
fig.text(text="CHINA", x=107.4, y=23, font="20p,Helvetica-Bold,black")

fig.plot(
    x=[103.12,103.34],
    y=[21.163,21.75],
    color="yellow",
    size=0.1*1.3** np.array([6.7,6.8]),
    style='a', 
    pen='black'
    )

#Plot text and arrow
trans="@10%"
fig.text(text="Ms6.7 (1983)", x=102.74, y=21.75, font="12p,Helvetica-Bold,black{}".format(trans),fill="white{}".format(trans))
fig.text(text="Ms6.8 (1935)", x=102.56, y=21.163, font="12p,Helvetica-Bold,black{}".format(trans),fill="white{}".format(trans))
fig.text(text="Mw5.0 (2020)", x=104.8, y=21.2, font="12p,Helvetica-Bold,black{}".format(trans),fill="white{}".format(trans))
# fig.plot(x=[105.0,dfc.Lon[0]],
#          y=[21.20,dfc.Lat[0]],
#          pen="1",
#         )
fig.basemap(map_scale="x6.0i/0.2i+c2+w100")
# save figure as pdf
fig.savefig("MC_eq.pdf", crop=True, dpi=600)
fig.show()


# In[17]:


dfc = df.iloc[-17:-16,:].reset_index()
l=0.2
minlon, maxlon = xc-l, xc+l
minlat, maxlat = yc-0.65*l, yc+0.65*l


#define etopo data file
topo_data = '@earth_relief_01s' #30 arc second global relief (SRTM15+V2.1 @ 1.0 km)


# Visualization
fig = pygmt.Figure()


# make color pallets
pygmt.makecpt(
    cmap='grayC',
    series='-100/10000/200',
    continuous=True
)

#plot high res topography
fig.grdimage(
    grid=topo_data,
    region=[minlon, maxlon, minlat, maxlat], 
    projection='M20c',
    shading=True,
    frame=True
    )


## Plot faults
#fig.plot(data="../topo/fault_lv1.txt", pen="0.75")
fig.plot(data="../Map/topo/fault_lv2.txt", pen="4,black",label="Fault")
fig.legend(position='JBR+jBR+o0.2c',box='+gwhite+p1p')
dfd = df.iloc[-17:-1,:]
dfd=dfd.sort_values("Lon")
pygmt.makecpt(cmap="jet", series=[dfd.Depth.min(), dfd.Depth.max()])
#plot data points

fig.plot(
    x=dfd.Lon,
    y=dfd.Lat, 
    size=0.1 * 1.4 ** dfd.Ml,
    style='c', 
    color='red', 
#    cmap=True,
    pen='black', 
    )
# fig.colorbar(
#     frame='+l"Earthquake depth (km)"'
#     )


for m in [3,4,5]:
    mag = 0.1 * 1.4 **m
    fig.plot(
        x=[0],
        y=[0],
#        color="blue",
        style='c{}'.format(mag), 
        pen='black',label="M={}".format(m)
        )

    
fig.legend()

colorl = ["purple","green","white","violet"]
for i,nw in enumerate(sta.Network.unique()):
    sta_nw = sta.query("Network == '{}'".format(nw))
    fig.plot(
    x=sta_nw.Longtitude,
    y=sta_nw.Latitude, 
    style='t0.6c', 
    color=colorl[i], 
    pen='black', 
    label="{}".format(nw),
    )

#focal = pd.read_csv("../infor_focal_mocchau.txt",delim_whitespace=True)
# plot_longitude = np.array([104.703, 104.708, 104.686, 104.714, 104.716, 104.683, 104.725,
#                            104.734, 104.716, 104.708, 104.693, 104.749, 104.708])
plot_latitude = np.array([20.981,20.861,20.981,20.861,20.981,20.861,20.981,20.861,20.981,20.861,20.981,20.861, 20.981])

#plot_longitude = np.linspace(104.64,104.79,13)
plot_longitude = [104.56,104.56 , 104.607, 104.613, 104.654, 104.659,104.701, 104.707, 104.748, 104.754 , 104.795, 104.815, 104.842]

#mrr, mtt, mff, mrt, mrf, mtf, exponent
focal = focalo.sort_values("lon").reset_index(drop=True)
def meca_plot():
    for i in range(0,len(plot_longitude)):
        fig.plot(x=[focal.lon.values[i],plot_longitude[i]]
                 ,y=[focal.lat.values[i],plot_latitude[i]]
                 , pen="1p,black"
                )
#     if focal.mw.values[i] == max(focal.mw):
        if plot_latitude[i] > 20.9:
            modlat = 0.017
        else:
            modlat = -0.017
        eventid = str(focal.event.values[i])[4:]
        eid = "{}/{}/{}:{}".format(eventid[0:2],eventid[2:4],eventid[4:6],eventid[6:])
        fig.text(x=plot_longitude[i],y=plot_latitude[i]+modlat,text=eid,fill="white",font="11p,Helvetica-Bold,black")
    for i in range(0,len(focal)):
        fig.meca(
dict(
#         mrr=focal.mrr.values,
#         mtt=focal.mtt.values,
#         mff=focal.mpp.values,
#         mrt=focal.mrt.values,
#         mrf=focal.mrp.values,
#         mtf=focal.mtp.values,
#         exponent=focal.exp.values,
       strike = focal.strike1.values[i],
       dip = focal.dip1.values[i],
       rake = focal.rake1.values[i],
       magnitude=focal.mw.values[i]
)
#             ,"-Feblue -S1c"
            , scale="1.5c"
#             , longitude=focal.lon.values[i]
#             , latitude=focal.lat.values[i]
            , longitude=plot_longitude[i]
            , latitude=plot_latitude[i]
            , depth=focal.dep.values[i]
            , offset = "0.1p,black", 
                )
meca_plot()   
fig.text(x=104.64838,y=20.825,text="MCVB",font="10p,Helvetica-Bold,black",fill="white")
fig.text(x=104.55,y=21.03,text="Da River fault",font="14p,Helvetica-Bold,black",fill="white",angle=-25)
fig.text(x=104.750,y=21.03,text="Muong La fault",font="14p,Helvetica-Bold,black",fill="white",angle=-20)
fig.basemap(map_scale="x6.0i/0.2i+c2+w5")    
# save figure as pdf
fig.savefig("MC_eq_local.pdf", crop=True, dpi=600)
fig.show()





# In[ ]:


dfd['time'] = pd.to_datetime(dfd[["year","month","day","hour","minute","second"]])
fig, ax = plt.subplots(figsize=(8, 4))
ax.scatter(dfd.time,dfd.Ml,s=3*2**dfd.Ml,c="r")
ax.set_ylabel("Magnitude (Ml)")
ax.set_xlabel("Time")
date_form = DateFormatter("%m-%d")
ax.xaxis.set_major_formatter(date_form)
ax.xaxis.set_major_locator(mdates.DayLocator(interval=2))
plt.savefig("eq_freq.jpg",dpi=300,bbox_inches="tight")
plt.show()


# In[15]:


it_p = pd.read_csv("../Map/intensity_interview.txt",delimiter="\t",names=["it","lat","lon"])
it_p


# In[16]:


color_dict = {3:"green", 4:"blue",5:"yellow",6:"orange",7:"red"}
it_dict = {3:"III", 4:"IV",5:"V",6:"VI",7:"VII"}


# In[22]:


dfc = df.iloc[-17:-16,:].reset_index()
minlon, maxlon = 104, 105.5
minlat, maxlat = 20.5, 21.5
topo_data = '@earth_relief_01s' #30 arc second global relief (SRTM15+V2.1 @ 1.0 km)

pga = pd.read_csv("pga.txt")
pga = pd.merge(pga[["sta","pga","mag"]],sta[["Station","Latitude","Longtitude"]],left_on="sta",right_on="Station")
pga = pga.query("mag == 5.2").reset_index(drop=True)

# Visualization
fig = pygmt.Figure()
pygmt.makecpt(
    cmap='grayC',
    series='-100/10000/200',
    continuous=True
)
fig.grdimage(
    grid=topo_data,
    region=[minlon, maxlon, minlat, maxlat], 
    projection='M16c',
    shading=True,
  
    )
fig.coast(
    region=[minlon, maxlon, minlat, maxlat], 
    shorelines=True,frame=True
    )
fig.coast(borders=["1/2.0p,orange"])
#fig.plot(data="../topo/fault_lv1.txt", pen="0.75")
#fig.plot(data="../topo/fault_lv2.txt", pen="0.5")
fig.plot(data="../Map/intensity.txt", pen="0.75",label='"Isoseismal line"')
lonc = dfc.Lon.values[0]
latc = dfc.Lat.values[0]

fig.plot(
    x=lonc,
    y=latc,
    color="black",
    style='a0.5c', 
    pen='black',
    label = "Mainshock"
    )

# fig.plot(
#     x=cities.Lon,
#     y=cities.Lat, 
#     style='d0.3c', 
#     color="white", 
#     pen='red',
#     label="Cities"
#     )
fig.legend(position='JBL+jBL+o0.2c',box='+gwhite+p1p')

#fig.plot(x=0,y=0,label ='"Intensity observation"')

for it in it_p.it.unique():
    it_t = it_p.query("it == '{}'".format(it))
    fig.plot(
    x=it_t.lon,
    y=it_t.lat, 
    style='d0.3c', 
    color="{}".format(color_dict[it]), 
    pen='black',
    label="{}".format(it_dict[it])
    )
fig.legend(position='JBR+jBR+o0.2c',box='+gwhite+p1p')
    
#Plot stations
# for i,nw in enumerate(sta.Network.unique()):
#     sta_nw = sta.query("Network == '{}'".format(nw))
#     fig.plot(
#     x=sta_nw.Longtitude,
#     y=sta_nw.Latitude, 
#     style='t0.8c', 
#     color="white", 
#     pen='black', 
#     )
    
pygmt.makecpt(cmap="jet", series=[pga.pga.min(), pga.pga.max()])

# fig.plot(
#     x=pga.Longtitude,
#     y=pga.Latitude, 
#     style='t0.8c', 
#     color=pga.pga, 
#     cmap=True,
#     pen='black', 
#     )
# fig.colorbar(
#     frame='+l"PGA(cm/s2)"'
#     )







fig.text(text="Hanoi", x=106.01, y=21.05, font="12p,Helvetica,black")
fig.text(text="LAOS", x=104, y=20, font="20p,Helvetica-Bold,white")
#fig.text(text="VIETNAM", x=105.5, y=21.7, font="20p,Helvetica-Bold,white")
fig.text(text="CHINA", x=107.4, y=23, font="20p,Helvetica-Bold,white")

fig.text(text="III", x=104.25, y=21, font="12p,Helvetica-Bold,black",fill="white")
fig.text(text="IV", x=104.42, y=21, font="12p,Helvetica-Bold,black",fill="white")
fig.text(text="V", x=104.60, y=21, font="12p,Helvetica-Bold,black",fill="white")
fig.text(text="VI", x=104.77, y=20.90, font="12p,Helvetica-Bold,black",fill="white")

# save figure as pdf
fig.savefig("MC_gm.pdf", crop=True, dpi=450)
fig.show()


# In[30]:


nm = (20.949987723141568, 104.72027426757619)
tlai = (20.959275500813618, 104.68024488143868)
tlap = (20.943315787784552, 104.6218054424378)
hp = (20.921709729008846, 104.74986429641112)
mc = (20.929,104.708)
loc_infra = [nm,tlai,tlap,hp]

lonc = dfc.Lon.values[0]
latc = dfc.Lat.values[0]

dfc = df.iloc[-17:-16,:].reset_index()
la=0.11
minlon, maxlon = lonc-la,lonc+la
minlat, maxlat = latc-la,latc+la
topo_data = '@earth_relief_01s' #30 arc second global relief (SRTM15+V2.1 @ 1.0 km)

# Visualization
fig = pygmt.Figure()
pygmt.makecpt(
    cmap='grayC',
    series='-100/10000/200',
    continuous=True
)
fig.grdimage(
    grid=topo_data,
    region=[minlon, maxlon, minlat, maxlat], 
    projection='M20c',
    shading=True,
    frame=True,
    )
fig.coast(
    region=[minlon, maxlon, minlat, maxlat], 
    shorelines=True,
    frame=True
    )
fig.coast(borders=["1/2.0p,orange"])
#fig.plot(data="../topo/fault_lv1.txt", pen="0.75")
#fig.plot(data="../topo/fault_lv2.txt", pen="0.5")

fig.plot(
    x=lonc,
    y=latc,
    color="black",
    style='a0.5c', 
    pen='black',
    label = "Mainshock"
    )

for loc in loc_infra:
    fig.plot(
    x=loc[1],
    y=loc[0],
    color="black",
    style='c0.1c', 
    pen='black',
    )


# fig.plot(
#     x=cities.Lon,
#     y=cities.Lat, 
#     style='d0.3c', 
#     color="white", 
#     pen='red',
#     label="Cities"
#     )
fig.legend(position='JBL+jBL+o0.2c',box='+gwhite+p1p')

fig.savefig("MC_local.pdf", crop=True, dpi=450)
fig.show()


# In[35]:


from doublecouple import *
focal = pd.read_csv("../MCHAU_SAC/fm.csv",index_col=0).reset_index(drop=True)
focal_w = focal[["event","mw","lat","lon","dep"
                 ,"strike1","dip1","rake1"
                 ,"strike2","dip2","rake2"]].sort_values("mw",ascending=False).reset_index(drop=True)
focal_w[["strike1","dip1","rake1"]] = focal_w[["strike1","dip1","rake1"]].astype(int).astype(str)
focal_w[["strike2","dip2","rake2"]] = focal_w[["strike2","dip2","rake2"]].astype(int).astype(str)
focal_w["nd1"] = focal_w[["strike1","dip1","rake1"]].agg('/'.join, axis=1)
focal_w["nd2"] = focal_w[["strike2","dip2","rake2"]].agg('/'.join, axis=1)
focal_w[["event","mw","lat","lon","dep"
                 ,"nd1"
                 ,"nd2"]].sort_values("event").reset_index(drop=True)

pazml, tazml, pdipl, tdipl  = [],[],[],[]
for i in range(0,len(focal)):
    dc = DoubleCouple([focal.strike1.values[i],focal.dip1.values[i],focal.rake1.values[i]])
    pazm = dc.axis['P']['azimuth']
    tazm = dc.axis['T']['azimuth']
    pdip = dc.axis['P']['dip']
    tdip = dc.axis['T']['dip']
#    print(pazm,tazm,pdip,tdip)
    pazml.append(pazm)
    tazml.append(tazm)
    pdipl.append(pdip)
    tdipl.append(tdip)
focal["T_AZM"] = tazml
focal["P_AZM"] = pazml
focal["P_PL"] = pdipl
focal["T_PL"] = tdipl


# In[37]:


# Visualization
xc,yc = 104.715,20.925    
l = 0.06

minlon, maxlon = xc-l, xc+l
minlat, maxlat = yc-l, yc+l
fig = pygmt.Figure()


# make color pallets
# pygmt.makecpt(
#     cmap='grayC',
#     series='-1000/5000/500',
#     continuous=True
# )

#plot high res topography
# fig.grdimage(
#     grid=topo_data,
#     region=[minlon, maxlon, minlat, maxlat], 
#     projection='M20c',
#     shading=True,
#     frame=True,
#     )

# plot coastlines
fig.coast(
    region=[minlon, maxlon, minlat, maxlat], 
    shorelines=True,
    frame=True
    )
fig.coast(borders=["1/2.0p,orange"])

## Plot faults
#fig.plot(data="../topo/fault_lv1.txt", pen="0.75")
fig.plot(data="../Map/topo/fault_lv2.txt", pen="2")

pygmt.makecpt(cmap="jet", series=[dfd.Depth.min(), dfd.Depth.max()])
#plot data points

fpx, fpy = 104.717,20.925 
lx = 5/110
az=55
xc1,yc1,xc2,yc2 = fpx-lx*cos(radians(az)), fpy-lx*sin(radians(az)),fpx+lx*cos(radians(az)), fpy+lx*sin(radians(az))
#xc1,yc1,xc2,yc2 = 104.69,20.883,104.74,20.96
fig.plot([xc1,xc2],[yc1,yc2])
fig.plot(xc1,yc1,style='c4p',color="black")
fig.plot(xc2,yc2,style='c4p',color="black")
fig.text(x=xc1+0.005,y=yc1,text="A",font="14p")
fig.text(x=xc2+0.005,y=yc2,text="A'",font="14p")
dist = ((xc1-xc2)**2+(yc1-yc2)**2)**0.5*110
print(dist)


for i in range(0,len(focal)):
    azm = radians(float(focal['T_AZM'].to_list()[i]))
    x0,y0 = focal['lon'].to_list()[i],focal['lat'].to_list()[i]
    l = cos(radians(pd.to_numeric(focal['T_PL'],errors="coerce")[i]))*0.01
    x1,y1 = x0 + sin(azm)*l, y0 + cos(azm)*l
    x2,y2 = x0 - sin(azm)*l, y0 - cos(azm)*l
    dep = focal.dep[i]
    fig.plot(
        x=np.array([x1,x2]),
        y=np.array([y1,y2]), 
        pen="4p,blue"
        )
    
    azm = radians(float(focal['P_AZM'].to_list()[i]))
    x0,y0 = focal['lon'].to_list()[i],focal['lat'].to_list()[i]
    l = cos(radians(pd.to_numeric(focal['P_PL'],errors="coerce")[i]))*0.01
    x1,y1 = x0 + sin(azm)*l, y0 + cos(azm)*l
    x2,y2 = x0 - sin(azm)*l, y0 - cos(azm)*l
    fig.plot(
        x=np.array([x1,x2]),
        y=np.array([y1,y2]), 
        pen="4p,darkgreen"
        )
fig.plot(x=[0,0],y=[0,0],pen="4p,darkgreen", label='"P axis"')
fig.plot(x=[0,0],y=[0,0],pen="4p,blue", label='"T axis"')


fig.legend()    


fig.plot(
    x=focal.lon,
    y=focal.lat, 
    style='c10p', 
    color="white",
    pen='black', 
    )


# fig.plot(fpx,fpy,style="c10p",color="red")
fig.savefig("P_T_axis.pdf", crop=True, dpi=300)
fig.show()


# In[47]:


dfc = df.iloc[-17:-16,:].reset_index()
l=0.11
minlon, maxlon = xc-l, xc+l
minlat, maxlat = yc-0.75*l, yc+0.75*l


#define etopo data file
topo_data = '@earth_relief_01s' #30 arc second global relief (SRTM15+V2.1 @ 1.0 km)


# Visualization
fig = pygmt.Figure()


# make color pallets
pygmt.makecpt(
    cmap='grayC',
    series='-100/10000/200',
    continuous=True
)

#plot high res topography
fig.grdimage(
    grid=topo_data,
    region=[minlon, maxlon, minlat, maxlat], 
    projection='M20c',
    shading=True,
    frame=True
    )


## Plot faults
#fig.plot(data="../topo/fault_lv1.txt", pen="0.75")
fig.plot(data="../Map/topo/fault_lv2.txt", pen="4,black",label="Fault")
fig.legend(position='JBR+jBR+o0.2c',box='+gwhite+p1p')
dfd = df.iloc[-17:-1,:]
dfd=dfd.sort_values("Lon")


for m in [3,4,5]:
    mag = 0.1 * 1.4 **m
    fig.plot(
        x=[0],
        y=[0],
#        color="blue",
        style='c{}'.format(mag), 
        pen='black',label="M={}".format(m)
        )

    
fig.legend()

plot_latitude = np.array([20.961,20.881,20.961,20.881,20.961,
                          20.881,20.961,20.881,20.961,20.881,20.961,20.881, 20.961])

plot_longitude = [104.64,104.64 , 104.665, 104.665, 104.69, 104.69,
                  104.715, 104.715, 104.74, 104.74 , 104.765, 104.765, 104.79,
                    ]

#mrr, mtt, mff, mrt, mrf, mtf, exponent
focal = focalo.sort_values("lon")
def meca_plot():
    for i in range(0,len(plot_longitude)):
        fig.plot([focal.lon.values[i],plot_longitude[i]]
                 ,[focal.lat.values[i],plot_latitude[i]]
                 , pen="1p,black"
                )
#     if focal.mw.values[i] == max(focal.mw):

    for i in range(0,len(focal)):
        fig.meca(
dict(
#         mrr=focal.mrr.values,
#         mtt=focal.mtt.values,
#         mff=focal.mpp.values,
#         mrt=focal.mrt.values,
#         mrf=focal.mrp.values,
#         mtf=focal.mtp.values,
#         exponent=focal.exp.values,
       strike = focal.strike1.values[i],
       dip = focal.dip1.values[i],
       rake = focal.rake1.values[i],
       magnitude=focal.mw.values[i]
)
#             ,"-Feblue -S1c"
            , scale="1.5c"
#             , longitude=focal.lon.values[i]
#             , latitude=focal.lat.values[i]
            , longitude=plot_longitude[i]
            , latitude=plot_latitude[i]
            , depth=focal.dep.values[i]
            , offset = "0.1p,black", 
                )
meca_plot()   
fig.text(x=104.62,y=20.975,text="Da River fault",font="14p,Helvetica-Bold,black",fill="white",angle=-30)
fig.text(x=104.750,y=21.03,text="Muong La fault",font="14p,Helvetica-Bold,black",fill="white",angle=-20)




from doublecouple import *
focal = pd.read_csv("../MCHAU_SAC/fm.csv",index_col=0).reset_index(drop=True)
focal_w = focal[["event","mw","lat","lon","dep"
                 ,"strike1","dip1","rake1"
                 ,"strike2","dip2","rake2"]].sort_values("mw",ascending=False).reset_index(drop=True)
focal_w[["strike1","dip1","rake1"]] = focal_w[["strike1","dip1","rake1"]].astype(int).astype(str)
focal_w[["strike2","dip2","rake2"]] = focal_w[["strike2","dip2","rake2"]].astype(int).astype(str)
focal_w["nd1"] = focal_w[["strike1","dip1","rake1"]].agg('/'.join, axis=1)
focal_w["nd2"] = focal_w[["strike2","dip2","rake2"]].agg('/'.join, axis=1)
focal_w[["event","mw","lat","lon","dep"
                 ,"nd1"
                 ,"nd2"]].sort_values("event").reset_index(drop=True)

pazml, tazml, pdipl, tdipl  = [],[],[],[]
for i in range(0,len(focal)):
    dc = DoubleCouple([focal.strike1.values[i],focal.dip1.values[i],focal.rake1.values[i]])
    pazm = dc.axis['P']['azimuth']
    tazm = dc.axis['T']['azimuth']
    pdip = dc.axis['P']['dip']
    tdip = dc.axis['T']['dip']
#    print(pazm,tazm,pdip,tdip)
    pazml.append(pazm)
    tazml.append(tazm)
    pdipl.append(pdip)
    tdipl.append(tdip)
focal["T_AZM"] = tazml
focal["P_AZM"] = pazml
focal["P_PL"] = pdipl
focal["T_PL"] = tdipl
fpx, fpy = 104.717,20.925 
lx = 5/110
az=55
xc1,yc1,xc2,yc2 = fpx-lx*cos(radians(az)), fpy-lx*sin(radians(az)),fpx+lx*cos(radians(az)), fpy+lx*sin(radians(az))
#xc1,yc1,xc2,yc2 = 104.69,20.883,104.74,20.96
# fig.plot([xc1,xc2],[yc1,yc2])
# fig.plot(xc1,yc1,style='c4p',color="black")
# fig.plot(xc2,yc2,style='c4p',color="black")
# fig.text(x=xc1+0.005,y=yc1,text="A",font="14p")
# fig.text(x=xc2+0.005,y=yc2,text="A'",font="14p")
# dist = ((xc1-xc2)**2+(yc1-yc2)**2)**0.5*110
# print(dist)


for i in range(0,len(focal)):
    azm = radians(float(focal['T_AZM'].to_list()[i]))
    x0,y0 = focal['lon'].to_list()[i],focal['lat'].to_list()[i]
    l = cos(radians(pd.to_numeric(focal['T_PL'],errors="coerce")[i]))*0.01
    x1,y1 = x0 + sin(azm)*l, y0 + cos(azm)*l
    x2,y2 = x0 - sin(azm)*l, y0 - cos(azm)*l
    dep = focal.dep[i]
    fig.plot(
        x=np.array([x1,x2]),
        y=np.array([y1,y2]), 
        pen="2p,blue"
        )
    
    azm = radians(float(focal['P_AZM'].to_list()[i]))
    x0,y0 = focal['lon'].to_list()[i],focal['lat'].to_list()[i]
    l = cos(radians(pd.to_numeric(focal['P_PL'],errors="coerce")[i]))*0.01
    x1,y1 = x0 + sin(azm)*l, y0 + cos(azm)*l
    x2,y2 = x0 - sin(azm)*l, y0 - cos(azm)*l
    fig.plot(
        x=np.array([x1,x2]),
        y=np.array([y1,y2]), 
        pen="2p,darkgreen"
        )
fig.plot(x=[0,0],y=[0,0],pen="2p,darkgreen", label='"P axis"')
fig.plot(x=[0,0],y=[0,0],pen="2p,blue", label='"T axis"')


fig.legend(position='JBL+jBL+o0.2c',box='+gwhite+p1p')    




pygmt.makecpt(cmap="jet", series=[dfd.Depth.min(), dfd.Depth.max()])
fig.plot(
    x=dfd.Lon,
    y=dfd.Lat, 
    size=0.1 * 1.4 ** dfd.Ml,
    style='c', 
    color=dfd.Depth, 
    cmap=True,
    pen='black', 
    )
fig.colorbar(
    frame='+l"Earthquake depth (km)"'
    )

# save figure as pdf
fig.savefig("MC_eq_combine.pdf", crop=True, dpi=300)
fig.show()


# In[ ]: