1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13

#Importing the required libraries
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.cluster import KMeans
import geopandas as gpd
from descartes import PolygonPatch
%matplotlib inline
from matplotlib.pylab import rcParams
import pygmt

1
2
3
4

#import the required dataset
#df_full = pd.read_csv("Data gempa bumi pulau jawa.csv")
df_full = pd.read_csv("ihsan_baru.csv")
df_full.head()

1
2
3
4
5
6
7
8

df_full['Date'] = pd.to_datetime(df_full['Date'],  errors='coerce')
# only keeping the columns we require
df = df_full[['Latitude','Longitude','Magnitudo']]
df_year = df_full[['Latitude','Longitude','Magnitudo','Date']]
# only keeping earthquake categorized as moderate
df = df[df.Magnitudo >= 3]
df_year = df_year[df_year.Magnitudo >= 3]
len(df_year),len(df)

1
2
3
4
5
6
7

#Finding the optimal K using the elbow method
Sum_of_squared_distances = []
K = range(1,10)
for k in K:
    km = KMeans(n_clusters=k)
    km = km.fit(df)    
    Sum_of_squared_distances.append(km.inertia_)

1

plt.plot(K, Sum_of_squared_distances, 'rx-')

1
2

kmeans = KMeans(n_clusters=2).fit(df)
df['Clusters'] = kmeans.labels_

1

map_1 = gpd.read_file('KB_JATIM.SHP')

1
2
3
4
5

#find the centroid location
centroid_1, centroid_2 = kmeans.cluster_centers_
centroid_1 = np.delete(centroid_1,2)
centroid_2 = np.delete(centroid_2,2)

1
2

params={'figure.figsize':(15,12)}
rcParams.update(params)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11

map_1.plot(column = 'NM_PROP')
sns.scatterplot(x=df.Longitude, y=df.Latitude, hue = df.Clusters ,palette='Set1',
                size=df.Magnitudo, sizes=(10,200), 
                    legend = 'brief')
plt.legend(loc=2, bbox_to_anchor=(0.80,1.30))

plt.scatter(centroid_1[1],centroid_1[0],color="yellow")
plt.scatter(centroid_2[1],centroid_2[0],color="yellow")
plt.axis('off')
plt.show()
plt.savefig("earthquake_clusters_1.png", dpi=200)

1
2
3
4
5
6
7
8
9

region = [
    df.Longitude.min() - 1,
    df.Longitude.max() + 1,
    df.Latitude.min() - 1,
    df.Latitude.max() + 1,
]

print(region)
print(data.head())

1
2
3
4
5

fig = pygmt.Figure()
fig.basemap(region=region, projection="M8i", frame=True)
fig.coast(land="white", water="skyblue")
fig.plot(x=df.Longitude, y=df.Latitude, style="c0.2c", color="red", pen="black")
fig.show()

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12

fig = pygmt.Figure()
fig.basemap(region=region, projection="M8i", frame=True)
fig.coast(land="white", water="skyblue")
fig.plot(
    x=df.Longitude,
    y=df.Latitude,
    sizes=0.02 * (1.8 ** df.Magnitudo),
    style="cc",
    color="red",
    pen="black",
)
fig.show()

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11

#map_1.plot(column = 'PROPINSI')
sns.scatterplot(x=df.Longitude, y=df.Latitude, hue = df.Clusters ,palette='Set1',
                size=df.Magnitudo, sizes=(10,200), 
                    legend = 'brief')
plt.legend(loc=2, bbox_to_anchor=(0.80,1.30))

plt.scatter(centroid_1[1],centroid_1[0],color="yellow")
plt.scatter(centroid_2[1],centroid_2[0],color="yellow")
plt.axis('off')
plt.show()
plt.savefig("earthquake_clusters_1.png", dpi=200)

1
2
3
4
5

df_2018 = df_full[(df_full['Date'].dt.year == 2018) & (df_full['Magnitudo']>=3)]

df_2019 = df_full[(df_full['Date'].dt.year == 2019) & (df_full['Magnitudo']>=3)]

df_2020 = df_full[(df_full['Date'].dt.year > 2020) & (df_full['Magnitudo']>=3)]

1
2

#Defining the location for April 8th 2019 earthquake
may_22_location = [110.78,-6.69]

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16

#map_1.plot(column = 'DATA')

plt.scatter(df_2018.Longitude, df_2018.Latitude,20,color='cyan',
            label="2018")
plt.scatter(df_2019.Longitude, df_2019.Latitude,20,label="2019",
            color='orange')
plt.scatter(df_2020.Longitude, df_2020.Latitude,20,label = "2020",
            color='purple')
plt.scatter(may_22_location[0],may_22_location[1],20,color='black',
            alpha=1, label="May 22 epicenter")

plt.legend(loc=2, bbox_to_anchor=(0.80,1.30))

plt.axis('Off')
plt.show()
plt.savefig("earthquake_clusters_2.png", dpi=200)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39

import pandas as pd
import pygmt

# Load sample iris data, and convert 'species' column to categorical dtype
df = pd.read_csv("https://github.com/mwaskom/seaborn-data/raw/master/iris.csv")
df["species"] = df.species.astype(dtype="category")

# Use pygmt.info to get region bounds (xmin, xmax, ymin, ymax, zmin, zmax)
# The below example will return a numpy array like [0., 3., 4., 8., 1., 7.]
region = pygmt.info(
    table=df[["petal_width", "sepal_length", "petal_length"]],  # x, y, z columns
    per_column=True,  # report output as a numpy array
    spacing="1/2/0.5",  # rounds x, y and z intervals by 1, 2 and 0.5 respectively
)

# Make our 3D scatter plot, coloring each of the 3 species differently
fig = pygmt.Figure()
pygmt.makecpt(cmap="cubhelix", color_model="+c", series=(0, 3, 1))
fig.plot3d(
    x=df.petal_width,
    y=df.sepal_length,
    z=df.petal_length,
    sizes=0.1 * df.sepal_width,  # Vary each symbol size according to a data column
    color=df.species.cat.codes.astype(int),  # Points colored by categorical number code
    cmap=True,  # Use colormap created by makecpt
    region=region,  # (xmin, xmax, ymin, ymax, zmin, zmax)
    frame=[
        "WsNeZ3",  # z axis label positioned on 3rd corner
        'xafg+l"Petal Width"',
        'yafg+l"Sepal Length"',
        'zafg+l"Petal Length"',
    ],
    style="uc",  # 3D cUbe, with size in centimeter units
    perspective=[315, 25],  # Azimuth NorthWest (315°), at elevation 25°
    zscale=1.5,  # Vertical exaggeration factor
)
fig.show()

1

!ls *.SHP

1

map_1 = gpd.read_file('KB_JATIM.SHP')

1

map_1.plot(column = 'NM_PROP')

1