Clustering Berdasarkan Bencana Alam dengan Menggunakan Berbagai Jenis Clustering

Library

Import Data

dt1<- read_xlsx("C:\\Users\\Acer\\Downloads\\Full Bencana Alam 2021-2023.xlsx",sheet="2021-2023")
dt1<-(dt1[1:34,1:9])%>%dplyr::select(-Kekeringan)%>%dplyr::select(-CuacaEkstrem)
OBJECTID<-seq(1:34)
dt1 <- dt1 %>% arrange(Provinsi)
dt.code<-cbind(OBJECTID,dt1)

Eksplorasi

summary(dt1)

##    Provinsi             Gempa        GunungMeletus     TanahLongsor   
##  Length:34          Min.   : 0.000   Min.   :0.0000   Min.   :   2.0  
##  Class :character   1st Qu.: 0.000   1st Qu.:0.0000   1st Qu.:   9.5  
##  Mode  :character   Median : 1.000   Median :0.0000   Median :  20.0  
##                     Mean   : 2.471   Mean   :0.1765   Mean   :  75.0  
##                     3rd Qu.: 4.000   3rd Qu.:0.0000   3rd Qu.:  43.5  
##                     Max.   :20.000   Max.   :2.0000   Max.   :1103.0  
##      Banjir         Kebakaran          Abrasi      
##  Min.   : 19.00   Min.   :  0.00   Min.   : 0.000  
##  1st Qu.: 58.75   1st Qu.:  7.75   1st Qu.: 1.000  
##  Median : 83.50   Median : 33.00   Median : 2.000  
##  Mean   :133.15   Mean   : 83.59   Mean   : 4.353  
##  3rd Qu.:171.50   3rd Qu.:107.00   3rd Qu.: 7.000  
##  Max.   :539.00   Max.   :631.00   Max.   :19.000

dtx <- dt1[,-1]
pairs.panels(dtx, method = "pearson", stars=TRUE)

par(mfrow=c(2,4))

# Looping Box-Plot
type <- names(dtx)[1:ncol(dtx)]
for (i in type){
  boxplot(dtx[,i],
       main = paste(i))
}

Prepocess

Scale

dtx <- scale(dtx)
head(dtx)

##           Gempa GunungMeletus TanahLongsor      Banjir   Kebakaran     Abrasi
## [1,] -0.6306381    -0.3847802   -0.2285273  1.42896544  0.98985398 -0.5176674
## [2,]  0.6456533    -0.3847802    0.1291676 -0.85602548 -0.55103790  0.5823758
## [3,] -0.3753798    -0.3847802   -0.2881431  0.03166883 -0.56734363 -0.9576846
## [4,] -0.3753798    -0.3847802   -0.2881431 -0.50259163 -0.59180223 -0.9576846
## [5,] -0.1201215    -0.3847802   -0.1490395 -0.93821940  0.02781566 -0.9576846
## [6,] -0.6306381    -0.3847802   -0.2533672 -0.79027035 -0.68148377 -0.9576846

PCA

pecea<-prcomp(dt1[,-1],scale=TRUE)
dtp<-pecea$x
head(dtp)

##             PC1        PC2         PC3        PC4          PC5         PC6
## [1,]  0.3531026 -1.4078984 -0.00992587 -0.5740510  1.097688398 -0.49629151
## [2,] -0.2148518  0.9036984 -0.68696473  0.4005886 -0.672338776  0.09964422
## [3,] -0.5628869 -0.3867522  0.30841977  0.8949552  0.444997376 -0.22699046
## [4,] -0.8631585 -0.3019800  0.28499325  0.9840254  0.057232191 -0.05425093
## [5,] -0.8013792 -0.6076745  0.09774839  0.7542980 -0.626454304  0.02108923
## [6,] -1.1517063 -0.2775669  0.32413617  1.0322685  0.008045391  0.20465687

K-Means

Scaled

Jumlah Cluster

set.seed(10)

fviz_nbclust(x = dtx, 
             FUNcluster = kmeans,
             method = 'wss',
             k.max = 10)

fviz_nbclust(x = dtx, 
             FUNcluster = kmeans,
             method = 'silhouette',
             k.max = 10)

fviz_nbclust(x = dtx, 
             FUNcluster = kmeans,
             method = 'gap_stat',
             k.max = 10)

kmeans.dt <-kmeans(dtx, centers = 3,nstart=10)
fviz_cluster(data = dtx, object = kmeans.dt)

table(kmeans.dt$cluster)

## 
##  1  2  3 
## 28  1  5

data.akhir3 <- cbind(dt1, Cluster =  kmeans.dt$cluster) %>% 
  relocate(Cluster, .before = 2)

# Radar Plot
ggRadar(
  data = data.akhir3,
  mapping = aes(colours = Cluster),
) + 
theme_light() +
theme(
  text = element_text(size = 10),  # Mengubah ukuran font global
  title = element_text(size = 12),  # Mengubah ukuran font judul
  axis.text = element_text(size = 10),  # Mengubah ukuran font label sumbu
  legend.text = element_text(size = 8)  # Mengubah ukuran font legenda
)

PCA

Jumlah Cluster

set.seed(10)

fviz_nbclust(x = dtp, 
             FUNcluster = kmeans,
             method = 'wss',
             k.max = 10)

fviz_nbclust(x = dtp, 
             FUNcluster = kmeans,
             method = 'silhouette',
             k.max = 10)

fviz_nbclust(x = dtp, 
             FUNcluster = kmeans,
             method = 'gap_stat',
             k.max = 10)

kmeans.dtp <-kmeans(dtp, centers = 5,nstart=10)
fviz_cluster(data = dtp, object = kmeans.dtp)

table(kmeans.dtp$cluster)

## 
##  1  2  3  4  5 
##  1  3 19  1 10

GMM

Scaled

mod1<-Mclust(dtx)
mod1$BIC

## Bayesian Information Criterion (BIC): 
##         EII       VII       EEI       VEI       EVI       VVI       EEE
## 1 -597.5214 -597.5214 -615.1532 -615.1532 -615.1532 -615.1532 -580.3502
## 2 -559.2172 -485.9441 -561.3111 -422.2024        NA        NA -591.0427
## 3 -566.0764 -457.6315 -571.1291 -443.2729        NA        NA -595.2520
## 4 -555.6605 -431.8038 -566.5335 -398.8138        NA        NA -607.2151
## 5 -500.9326        NA -488.4247        NA        NA        NA -562.2093
## 6 -522.3708        NA -496.5711        NA        NA        NA -589.4288
## 7 -532.6870        NA -518.9977        NA        NA        NA -567.6885
## 8 -543.6748        NA -533.2364        NA        NA        NA -538.5200
## 9 -496.7701        NA -492.3112        NA        NA        NA -514.1277
##         VEE       EVE       VVE       EEV       VEV       EVV       VVV
## 1 -580.3502 -580.3502 -580.3502 -580.3502 -580.3502 -580.3502 -580.3502
## 2 -497.1977        NA        NA        NA -464.4904        NA        NA
## 3        NA        NA        NA        NA -431.8205        NA        NA
## 4        NA        NA        NA        NA -456.7600        NA        NA
## 5        NA        NA        NA        NA        NA        NA        NA
## 6        NA        NA        NA        NA        NA        NA        NA
## 7        NA        NA        NA        NA        NA        NA        NA
## 8        NA        NA        NA        NA        NA        NA        NA
## 9        NA        NA        NA        NA        NA        NA        NA
## 
## Top 3 models based on the BIC criterion: 
##     VEI,4     VEI,2     VII,4 
## -398.8138 -422.2024 -431.8038

plot(mod1, what = 'BIC')

Model GMM

mod1b = Mclust(dtx, G = 3, modelNames = c("VEV"))
summary(mod1b, parameters = TRUE)

## ---------------------------------------------------- 
## Gaussian finite mixture model fitted by EM algorithm 
## ---------------------------------------------------- 
## 
## Mclust VEV (ellipsoidal, equal shape) model with 3 components: 
## 
##  log-likelihood  n df       BIC       ICL
##       -87.19811 34 73 -431.8205 -431.8206
## 
## Clustering table:
##  1  2  3 
## 14  5 15 
## 
## Mixing probabilities:
##         1         2         3 
## 0.4117655 0.1470587 0.4411758 
## 
## Means:
##                    [,1]       [,2]        [,3]
## Gempa         0.4997899  0.2882914 -0.56256935
## GunungMeletus 0.5496842 -0.3847802 -0.38478024
## TanahLongsor  0.3704685 -0.2295208 -0.26926482
## Banjir        0.8013541 -0.6604042 -0.52779832
## Kebakaran     0.3085065 -0.6423500 -0.07382371
## Abrasi        0.3937955  0.8463860 -0.64967250
## 
## Variances:
## [,,1]
##                    Gempa GunungMeletus TanahLongsor     Banjir  Kebakaran
## Gempa          5.2338895     0.6372570    5.9462652  3.5765011 -0.3741846
## GunungMeletus  0.6372570     3.5032271    0.2199069  0.3005589 -1.3091277
## TanahLongsor   5.9462652     0.2199069    7.7174393  4.9740211  0.9975762
## Banjir         3.5765011     0.3005589    4.9740211  3.5027439  0.9908961
## Kebakaran     -0.3741846    -1.3091277    0.9975762  0.9908961  3.7700381
## Abrasi        -1.1700194    -0.5191740   -2.0729339 -1.6755291 -1.2691411
##                  Abrasi
## Gempa         -1.170019
## GunungMeletus -0.519174
## TanahLongsor  -2.072934
## Banjir        -1.675529
## Kebakaran     -1.269141
## Abrasi         2.370898
## [,,2]
##                       Gempa GunungMeletus  TanahLongsor        Banjir
## Gempa          7.505166e-02 -3.909436e-19  1.705983e-02  7.223056e-04
## GunungMeletus -3.909436e-19  2.950565e-04  3.099152e-18  1.929974e-18
## TanahLongsor   1.705983e-02  3.099152e-18  1.903993e-02 -8.478868e-03
## Banjir         7.223056e-04  1.929974e-18 -8.478868e-03  1.711736e-02
## Kebakaran      4.702317e-03  2.069319e-18  4.495385e-03 -1.809503e-03
## Abrasi        -1.966505e-02 -3.885325e-18 -1.408811e-02  1.974907e-03
##                   Kebakaran        Abrasi
## Gempa          4.702317e-03 -1.966505e-02
## GunungMeletus  2.069319e-18 -3.885325e-18
## TanahLongsor   4.495385e-03 -1.408811e-02
## Banjir        -1.809503e-03  1.974907e-03
## Kebakaran      3.139205e-03 -4.849126e-03
## Abrasi        -4.849126e-03  2.755112e-02
## [,,3]
##                       Gempa GunungMeletus  TanahLongsor        Banjir
## Gempa          1.824254e-02 -2.879784e-19  2.598657e-03 -2.426079e-03
## GunungMeletus -2.879784e-19  8.964594e-04 -1.433044e-17 -7.836678e-18
## TanahLongsor   2.598657e-03 -1.433044e-17  1.028641e-02  1.235071e-03
## Banjir        -2.426079e-03 -7.836678e-18  1.235071e-03  6.507615e-02
## Kebakaran     -9.808067e-03 -1.718139e-17  1.499030e-02  9.421723e-03
## Abrasi        -1.531552e-02 -1.864893e-17 -9.175317e-03 -1.656528e-02
##                   Kebakaran        Abrasi
## Gempa         -9.808067e-03 -1.531552e-02
## GunungMeletus -1.718139e-17 -1.864893e-17
## TanahLongsor   1.499030e-02 -9.175317e-03
## Banjir         9.421723e-03 -1.656528e-02
## Kebakaran      2.071361e-01  9.325669e-02
## Abrasi         9.325669e-02  1.303863e-01

table(mod1b[["classification"]])

## 
##  1  2  3 
## 14  5 15

fviz_cluster(mod1b, data = dtx, repel = TRUE, labelsize =8)

data.akhir1 <- cbind(dt1, Cluster =  mod1b[["classification"]]) %>% 
  relocate(Cluster, .before = 2)

# Radar Plot
ggRadar(
  data = data.akhir1,
  mapping = aes(colours = Cluster)
) + 
theme_light() +
theme(
  text = element_text(size = 10),  # Mengubah ukuran font global
  title = element_text(size = 12),  # Mengubah ukuran font judul
  axis.text = element_text(size = 10),  # Mengubah ukuran font label sumbu
  legend.text = element_text(size = 8)  # Mengubah ukuran font legenda
)

PCA

mod2<-Mclust(dtp)
mod2$BIC

## Bayesian Information Criterion (BIC): 
##         EII       VII       EEI       VEI       EVI       VVI       EEE
## 1 -597.5214 -597.5214 -527.4548 -527.4548 -527.4548 -527.4548 -580.3502
## 2 -559.2172 -485.9441 -552.1659 -467.8515 -516.0034 -455.3572 -591.0427
## 3 -566.0785 -457.6316 -538.3366 -449.0898 -517.0487 -441.4360 -595.2507
## 4 -555.6606 -431.8043 -532.7229 -446.8640 -568.6461 -455.4909 -607.2150
## 5 -532.9434 -421.3234 -533.1886 -442.9874 -582.7761 -431.1951 -628.4155
## 6 -539.4300        NA -541.5713        NA        NA        NA -575.3024
## 7 -495.8063        NA -518.6842        NA        NA        NA -472.6681
## 8 -495.6143        NA -532.2559        NA        NA        NA -467.4861
## 9 -489.8217        NA -500.9381        NA        NA        NA        NA
##         VEE       EVE       VVE       EEV       VEV       EVV       VVV
## 1 -580.3502 -580.3502 -580.3502 -580.3502 -580.3502 -580.3502 -580.3502
## 2        NA        NA        NA        NA -464.4904        NA        NA
## 3        NA        NA        NA        NA -451.6115        NA        NA
## 4        NA        NA        NA        NA -474.8109        NA        NA
## 5        NA        NA        NA        NA -478.2068        NA        NA
## 6        NA        NA        NA        NA        NA        NA        NA
## 7        NA        NA        NA        NA        NA        NA        NA
## 8        NA        NA        NA        NA        NA        NA        NA
## 9        NA        NA        NA        NA        NA        NA        NA
## 
## Top 3 models based on the BIC criterion: 
##     VII,5     VVI,5     VII,4 
## -421.3234 -431.1951 -431.8043

Model GMM

mod2b = Mclust(dtp, G = 5, modelNames = c("VVI"))
summary(mod2b, parameters = TRUE)

## ---------------------------------------------------- 
## Gaussian finite mixture model fitted by EM algorithm 
## ---------------------------------------------------- 
## 
## Mclust VVI (diagonal, varying volume and shape) model with 5 components: 
## 
##  log-likelihood  n df       BIC       ICL
##        -102.754 34 64 -431.1951 -431.3918
## 
## Clustering table:
##  1  2  3  4  5 
## 12  5  6  9  2 
## 
## Mixing probabilities:
##          1          2          3          4          5 
## 0.35363445 0.14699430 0.17801236 0.26253536 0.05882353 
## 
## Means:
##            [,1]        [,2]        [,3]         [,4]       [,5]
## PC1  1.27767612 -0.50216847 -0.81126676 -0.925799580  0.1607562
## PC2 -0.07204943  1.04865203 -0.47882175 -0.522310111  1.5928040
## PC3 -0.05099231 -0.64607547  0.26493279 -0.008176555  1.1557854
## PC4 -0.45341629  0.09222485  0.72722865  0.248125463 -0.8127751
## PC5  0.20169009 -0.26654844  0.23098831 -0.195447081 -0.3731579
## PC6 -0.18167030  0.01250550 -0.04271319  0.190284184  0.3409144
## 
## Variances:
## [,,1]
##          PC1      PC2      PC3      PC4       PC5       PC6
## PC1 4.279596 0.000000 0.000000 0.000000 0.0000000 0.0000000
## PC2 0.000000 2.165995 0.000000 0.000000 0.0000000 0.0000000
## PC3 0.000000 0.000000 2.047496 0.000000 0.0000000 0.0000000
## PC4 0.000000 0.000000 0.000000 1.371381 0.0000000 0.0000000
## PC5 0.000000 0.000000 0.000000 0.000000 0.5513876 0.0000000
## PC6 0.000000 0.000000 0.000000 0.000000 0.0000000 0.1798848
## [,,2]
##            PC1        PC2         PC3        PC4        PC5        PC6
## PC1 0.07107456 0.00000000 0.000000000 0.00000000 0.00000000 0.00000000
## PC2 0.00000000 0.02233341 0.000000000 0.00000000 0.00000000 0.00000000
## PC3 0.00000000 0.00000000 0.005781948 0.00000000 0.00000000 0.00000000
## PC4 0.00000000 0.00000000 0.000000000 0.04247175 0.00000000 0.00000000
## PC5 0.00000000 0.00000000 0.000000000 0.00000000 0.06682667 0.00000000
## PC6 0.00000000 0.00000000 0.000000000 0.00000000 0.00000000 0.04499671
## [,,3]
##            PC1        PC2         PC3        PC4       PC5         PC6
## PC1 0.02485197 0.00000000 0.000000000 0.00000000 0.0000000 0.000000000
## PC2 0.00000000 0.02567829 0.000000000 0.00000000 0.0000000 0.000000000
## PC3 0.00000000 0.00000000 0.005391526 0.00000000 0.0000000 0.000000000
## PC4 0.00000000 0.00000000 0.000000000 0.07261159 0.0000000 0.000000000
## PC5 0.00000000 0.00000000 0.000000000 0.00000000 0.0140209 0.000000000
## PC6 0.00000000 0.00000000 0.000000000 0.00000000 0.0000000 0.009211555
## [,,4]
##           PC1       PC2        PC3       PC4        PC5         PC6
## PC1 0.0504685 0.0000000 0.00000000 0.0000000 0.00000000 0.000000000
## PC2 0.0000000 0.1376327 0.00000000 0.0000000 0.00000000 0.000000000
## PC3 0.0000000 0.0000000 0.05760935 0.0000000 0.00000000 0.000000000
## PC4 0.0000000 0.0000000 0.00000000 0.3322461 0.00000000 0.000000000
## PC5 0.0000000 0.0000000 0.00000000 0.0000000 0.04159175 0.000000000
## PC6 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000 0.008026406
## [,,5]
##            PC1          PC2         PC3          PC4          PC5        PC6
## PC1 0.08756675 0.000000e+00 0.000000000 0.0000000000 0.000000e+00 0.00000000
## PC2 0.00000000 1.756317e-05 0.000000000 0.0000000000 0.000000e+00 0.00000000
## PC3 0.00000000 0.000000e+00 0.002990541 0.0000000000 0.000000e+00 0.00000000
## PC4 0.00000000 0.000000e+00 0.000000000 0.0002528181 0.000000e+00 0.00000000
## PC5 0.00000000 0.000000e+00 0.000000000 0.0000000000 1.914202e-05 0.00000000
## PC6 0.00000000 0.000000e+00 0.000000000 0.0000000000 0.000000e+00 0.02584654

fviz_cluster(mod2b, data = dtp, repel = TRUE, labelsize =8)

table(mod2b[["classification"]])

## 
##  1  2  3  4  5 
## 12  5  6  9  2

data.akhir2 <- cbind(dt1, Cluster =  mod2b[["classification"]]) %>% 
  relocate(Cluster, .before = 2)

# Radar Plot
ggRadar(
  data = data.akhir2,
  mapping = aes(colours = Cluster)
) + 
theme_light() +
theme(
  text = element_text(size = 10),  # Mengubah ukuran font global
  title = element_text(size = 12),  # Mengubah ukuran font judul
  axis.text = element_text(size = 10),  # Mengubah ukuran font label sumbu
  legend.text = element_text(size = 8)  # Mengubah ukuran font legenda
)

Hierarchy

Scaled

hc <- hclust(dist(dtx), method = "complete")
plot(hc)

clusters <- cutree(hc, k = 4)

# Visualize the clusters
fviz_cluster(list(data = dtx, cluster = clusters),
             labelsize = 8,
             repel=TRUE,
             ellipse.type = "convex", 
             geom = "point", 
             stand = FALSE, 
             show.clust.cent = FALSE,
             ggtheme = theme_minimal())+
  geom_text(aes(label = 1:34), size = 3, vjust = -0.5)

Fuzzy C

# Apply Fuzzy C-Means clustering
efceem <- fcm(dtx, centers = 3, m = 2)

table(efceem$cluster)

## 
##  1  2  3 
## 18  3 13

# Assign each point to the cluster with the highest membership value
clusters <- apply(efceem$u, 1, which.max)

# Example data, assuming you have your original data in a variable called data
data_df <- as.data.frame(dtx)
data_df$Cluster <- as.factor(clusters)

# Visualize the clusters
fviz_cluster(list(data = dtx, cluster = clusters), 
             geom = "point", 
             ellipse.type = "convex", 
             stand = FALSE, 
             show.clust.cent = TRUE, 
             ggtheme = theme_minimal())+
  geom_text(aes(label = 1:34), size = 3, vjust = -0.5)

# Extract cluster centers
cluster_centers <- as.data.frame(efceem$v)
cluster_centers$Cluster <- 1:nrow(cluster_centers)

# Normalize the data for ggradar
max_val <- apply(cluster_centers[,-ncol(cluster_centers)], 2, max)
min_val <- apply(cluster_centers[,-ncol(cluster_centers)], 2, min)

cluster_centers[,-ncol(cluster_centers)] <- sweep(cluster_centers[,-ncol(cluster_centers)], 2, max_val, "/")
cluster_centers[,-ncol(cluster_centers)] <- sweep(cluster_centers[,-ncol(cluster_centers)], 2, min_val, "-")

# Bind the data with a row of maximum values for scaling
max_values <- data.frame(matrix(1, ncol = ncol(cluster_centers) - 1, nrow = 1))
colnames(max_values) <- colnames(cluster_centers)[-ncol(cluster_centers)]
max_values$Cluster <- "Max"


# Create the radar plot
# Radar Plot
ggRadar(
  data = cluster_centers,
  mapping = aes(colours = Cluster)
) + 
theme_light() +
theme(
  text = element_text(size = 10),  # Mengubah ukuran font global
  title = element_text(size = 12),  # Mengubah ukuran font judul
  axis.text = element_text(size = 10),  # Mengubah ukuran font label sumbu
  legend.text = element_text(size = 8)  # Mengubah ukuran font legenda
)

Hasil Akhir

Hasil Pakai Fuzzy C

Peta

GMM

prov <- st_read("C:\\Users\\Acer\\Documents\\BATAS PROVINSI DESEMBER 2019 DUKCAPIL\\BATAS_PROVINSI_DESEMBER_2019_DUKCAPIL.shp")

## Reading layer `BATAS_PROVINSI_DESEMBER_2019_DUKCAPIL' from data source 
##   `C:\Users\Acer\Documents\BATAS PROVINSI DESEMBER 2019 DUKCAPIL\BATAS_PROVINSI_DESEMBER_2019_DUKCAPIL.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 34 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XYZ
## Bounding box:  xmin: 95.01079 ymin: -11.00762 xmax: 141.0194 ymax: 6.07693
## z_range:       zmin: 2.65e-05 zmax: 2.65e-05
## Geodetic CRS:  WGS 84

data.cluster <- cbind(dt.code, Cluster = mod1b[["classification"]])

data.prov <- prov %>%
  inner_join(data.cluster, by = c("OBJECTID" = "OBJECTID"))

ggplot() +  
  geom_sf(data=data.prov, aes(fill=factor(`Cluster`))) +
  scale_fill_manual(values=c("1" = "indianred", "2" = "darkseagreen", "3" = "dodgerblue", "4"="lightgreen","5"="darkorchid1"), 
                    name = "Keterangan") +
  labs(title = "Cluster Bencana Alam \n pada Provinsi Indonesia 2021-2024",
       x = "Longitude",
       y = "Latitude") +
  theme_minimal() +
  theme(legend.text = element_text(size=10),
        legend.title = element_text(size=10, face="bold"),
        axis.text.x = element_text(size = 10),
        axis.text.y = element_text(size = 10),
        plot.title = element_text(size=12, face="bold", hjust = 0.5)) +
  scale_x_continuous(labels = function(x) paste0(x, "°")) +
  scale_y_continuous(labels = function(y) paste0(y, "°"))

Fuzzy C

data.cluster1<-cbind(dt.code, Cluster= efceem$cluster)
data.prov1 <- prov %>%
  inner_join(data.cluster1, by = c("OBJECTID" = "OBJECTID"))

ggplot() +  
  geom_sf(data=data.prov1, aes(fill=factor(`Cluster`))) +
  scale_fill_manual(values=c("1" = "lightgreen", "2" = "lightyellow", "3" = "red"), 
                    name = "Keterangan") +
  labs(title = "Cluster Bencana Alam \n pada Provinsi Indonesia 2021-2024",
       x = "Longitude",
       y = "Latitude") +
  theme_minimal() +
  theme(legend.text = element_text(size=10),
        legend.title = element_text(size=10, face="bold"),
        axis.text.x = element_text(size = 10),
        axis.text.y = element_text(size = 10),
        plot.title = element_text(size=12, face="bold", hjust = 0.5),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank()) +
  scale_x_continuous(labels = function(x) paste0(x, "°")) +
  scale_y_continuous(labels = function(y) paste0(y, "°"))

ggplot() +  
  geom_sf(data=data.prov1, aes(fill=factor(`Cluster`))) +
  scale_fill_manual(values=c("1" = "#008000", "2" = "#FFDF00", "3" = "#4D0000"), 
                    name = "Keterangan") +
  labs(title = "Cluster Bencana Alam \n pada Provinsi Indonesia 2021-2024",
       x = "Longitude",
       y = "Latitude") +
  theme_minimal() +
  theme(legend.text = element_text(size=10),
        legend.title = element_text(size=10, face="bold"),
        axis.text.x = element_text(size = 10),
        axis.text.y = element_text(size = 10),
        plot.title = element_text(size=12, face="bold", hjust = 0.5),
        panel.grid.major = element_blank()) +
  scale_x_continuous(labels = function(x) paste0(x, "°")) +
  scale_y_continuous(labels = function(y) paste0(y, "°"))

ggplot() +  
  geom_sf(data=data.prov1, aes(fill=factor(`Cluster`))) +
  scale_fill_manual(values=c("1" = "#E6E6FA", "2" = "#9A737D", "3" = "#4D0000"), 
                    name = "Keterangan") +
  labs(title = "Cluster Bencana Alam \n pada Provinsi Indonesia 2021-2024",
       x = "Longitude",
       y = "Latitude") +
  theme_minimal() +
  theme(legend.text = element_text(size=10),
        legend.title = element_text(size=10, face="bold"),
        axis.text.x = element_text(size = 10),
        axis.text.y = element_text(size = 10),
        plot.title = element_text(size=12, face="bold", hjust = 0.5),
        panel.grid.major = element_blank()) +
  scale_x_continuous(labels = function(x) paste0(x, "°")) +
  scale_y_continuous(labels = function(y) paste0(y, "°"))