knitr::opts_chunk$set(
echo = TRUE
)
# ==============================================================================
# LOAD PACKAGE
# ==============================================================================
library(sf)
library(dplyr)
library(ggplot2)
library(spdep)
library(spatstat)
library(gstat)
library(sp)
library(viridis)
library(readr)
library(stringr)
# ==============================================================================
# PENDAHULUAN
# ==============================================================================
cat("Praktikum Statistika Spasial\n")
## Praktikum Statistika Spasial
cat("Analisis Kemiskinan dan Curah Hujan Jawa Tengah\n\n")
## Analisis Kemiskinan dan Curah Hujan Jawa Tengah
# ==============================================================================
# MENENTUKAN LOKASI FILE
# ==============================================================================
path_shp <- "C:/Users/Lutfi/Downloads/shp_jateng/shp_gadm/gadm41_IDN_2.shp"
path_kemiskinan <- "C:/Users/Lutfi/Downloads/shp_jateng/Kemiskinan Menurut Kabupaten_Kota di Provinsi Jawa Tengah, 2024.csv"
path_curahhujan <- "C:/Users/Lutfi/Downloads/shp_jateng/POWER_Regional_Monthly_2023_2024.csv"
# ==============================================================================
# BAB 1 - PERSIAPAN DATA
# ==============================================================================
cat("\n====================================================\n")
##
## ====================================================
cat("BAB 1 - PERSIAPAN DATA\n")
## BAB 1 - PERSIAPAN DATA
cat("====================================================\n")
## ====================================================
# ------------------------------------------------------------------------------
# 1A. MEMBACA SHAPEFILE
# ------------------------------------------------------------------------------
cat("\nMembaca shapefile Indonesia...\n")
##
## Membaca shapefile Indonesia...
shp_indonesia <- st_read(path_shp)
## Reading layer `gadm41_IDN_2' from data source
## `C:\Users\Lutfi\Downloads\shp_jateng\shp_gadm\gadm41_IDN_2.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 502 features and 13 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 95.00971 ymin: -11.00761 xmax: 141.0194 ymax: 6.076941
## Geodetic CRS: WGS 84
cat("\nNama kolom:\n")
##
## Nama kolom:
print(names(shp_indonesia))
## [1] "GID_2" "GID_0" "COUNTRY" "GID_1" "NAME_1" "NL_NAME_1"
## [7] "NAME_2" "VARNAME_2" "NL_NAME_2" "TYPE_2" "ENGTYPE_2" "CC_2"
## [13] "HASC_2" "geometry"
# Filter Jawa Tengah
shp_jateng <- shp_indonesia %>%
filter(NAME_1 == "Jawa Tengah")
cat("\nJumlah kabupaten/kota Jawa Tengah:\n")
##
## Jumlah kabupaten/kota Jawa Tengah:
print(nrow(shp_jateng))
## [1] 36
print(shp_jateng[, c("NAME_2", "geometry")])
## Simple feature collection with 36 features and 1 field
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 108.5559 ymin: -8.21193 xmax: 111.6914 ymax: -5.725698
## Geodetic CRS: WGS 84
## First 10 features:
## NAME_2 geometry
## 1 Banjarnegara MULTIPOLYGON (((109.42 -7.5...
## 2 Banyumas MULTIPOLYGON (((109.4411 -7...
## 3 Batang MULTIPOLYGON (((109.917 -7....
## 4 Blora MULTIPOLYGON (((111.4427 -7...
## 5 Boyolali MULTIPOLYGON (((110.534 -7....
## 6 Brebes MULTIPOLYGON (((109.0322 -7...
## 7 Cilacap MULTIPOLYGON (((109.3856 -7...
## 8 Demak MULTIPOLYGON (((110.5618 -7...
## 9 Grobogan MULTIPOLYGON (((110.9202 -7...
## 10 Jepara MULTIPOLYGON (((110.7048 -6...
# Peta dasar Jawa Tengah
plot(
st_geometry(shp_jateng),
main = "Peta Jawa Tengah - 35 Kabupaten/Kota"
)
# ------------------------------------------------------------------------------
# 1B. MEMBACA DATA KEMISKINAN
# ------------------------------------------------------------------------------
cat("\nMembaca data kemiskinan...\n")
##
## Membaca data kemiskinan...
kemiskinan_raw <- read_csv(
path_kemiskinan,
skip = 3,
col_names = FALSE
)
colnames(kemiskinan_raw) <- c(
"nama_wilayah",
"garis_kemiskinan",
"jumlah_miskin_ribu",
"persen_miskin"
)
cat("\nData awal kemiskinan:\n")
##
## Data awal kemiskinan:
print(head(kemiskinan_raw))
## # A tibble: 6 × 4
## nama_wilayah garis_kemiskinan jumlah_miskin_ribu persen_miskin
## <chr> <dbl> <dbl> <dbl>
## 1 <NA> 2024 2024 2024
## 2 3300 PROVINSI JAWA TENGAH 507001 3704. 10.5
## 3 3301 Kabupaten Cilacap 441093 186. 10.7
## 4 3302 Kabupaten Banyumas 500861 208. 12.0
## 5 3303 Kabupaten Purbalingga 460870 137. 14.2
## 6 3304 Kabupaten Banjarnegara 398344 138. 14.7
# Membersihkan data
kemiskinan <- kemiskinan_raw %>%
slice(-1) %>%
filter(!is.na(nama_wilayah)) %>%
mutate(
nama_bersih = str_trim(
str_replace(nama_wilayah, "^\\d+ ", "")
),
persen_miskin = as.numeric(persen_miskin),
jumlah_miskin_ribu = as.numeric(jumlah_miskin_ribu),
garis_kemiskinan = as.numeric(garis_kemiskinan)
)
cat("\nData setelah dibersihkan:\n")
##
## Data setelah dibersihkan:
print(kemiskinan[, c("nama_bersih", "persen_miskin")])
## # A tibble: 36 × 2
## nama_bersih persen_miskin
## <chr> <dbl>
## 1 PROVINSI JAWA TENGAH 10.5
## 2 Kabupaten Cilacap 10.7
## 3 Kabupaten Banyumas 12.0
## 4 Kabupaten Purbalingga 14.2
## 5 Kabupaten Banjarnegara 14.7
## 6 Kabupaten Kebumen 15.7
## 7 Kabupaten Purworejo 10.9
## 8 Kabupaten Wonosobo 15.3
## 9 Kabupaten Magelang 10.8
## 10 Kabupaten Boyolali 9.63
## # ℹ 26 more rows
# ------------------------------------------------------------------------------
# 1C. PENYESUAIAN NAMA UNTUK JOIN
# ------------------------------------------------------------------------------
cat("\nNama wilayah SHP:\n")
##
## Nama wilayah SHP:
print(sort(shp_jateng$NAME_2))
## [1] "Banjarnegara" "Banyumas" "Batang" "Blora"
## [5] "Boyolali" "Brebes" "Cilacap" "Demak"
## [9] "Grobogan" "Jepara" "Karanganyar" "Kebumen"
## [13] "Kendal" "Klaten" "Kota Magelang" "Kota Pekalongan"
## [17] "Kota Semarang" "Kota Tegal" "Kudus" "Magelang"
## [21] "Pati" "Pekalongan" "Pemalang" "Purbalingga"
## [25] "Purworejo" "Rembang" "Salatiga" "Semarang"
## [29] "Sragen" "Sukoharjo" "Surakarta" "Tegal"
## [33] "Temanggung" "Waduk Kedungombo" "Wonogiri" "Wonosobo"
cat("\nNama wilayah BPS:\n")
##
## Nama wilayah BPS:
print(sort(kemiskinan$nama_bersih))
## [1] "Kabupaten Banjarnegara" "Kabupaten Banyumas" "Kabupaten Batang"
## [4] "Kabupaten Blora" "Kabupaten Boyolali" "Kabupaten Brebes"
## [7] "Kabupaten Cilacap" "Kabupaten Demak" "Kabupaten Grobogan"
## [10] "Kabupaten Jepara" "Kabupaten Karanganyar" "Kabupaten Kebumen"
## [13] "Kabupaten Kendal" "Kabupaten Klaten" "Kabupaten Kudus"
## [16] "Kabupaten Magelang" "Kabupaten Pati" "Kabupaten Pekalongan"
## [19] "Kabupaten Pemalang" "Kabupaten Purbalingga" "Kabupaten Purworejo"
## [22] "Kabupaten Rembang" "Kabupaten Semarang" "Kabupaten Sragen"
## [25] "Kabupaten Sukoharjo" "Kabupaten Tegal" "Kabupaten Temanggung"
## [28] "Kabupaten Wonogiri" "Kabupaten Wonosobo" "Kota Magelang"
## [31] "Kota Pekalongan" "Kota Salatiga" "Kota Semarang"
## [34] "Kota Surakarta" "Kota Tegal" "PROVINSI JAWA TENGAH"
# Hapus Waduk Kedungombo
shp_jateng <- shp_jateng %>%
filter(NAME_2 != "Waduk Kedungombo")
# Membuat nama kunci
shp_jateng <- shp_jateng %>%
mutate(
nama_kunci = str_to_lower(
str_trim(
str_replace_all(
NAME_2,
c("Kabupaten " = "", "Kota " = "")
)
)
)
) %>%
mutate(
nama_kunci = ifelse(
NAME_2 == "Kota Magelang",
"kota magelang",
ifelse(
NAME_2 == "Magelang",
"kabupaten magelang",
nama_kunci
)
)
)
kemiskinan <- kemiskinan %>%
mutate(
nama_kunci = str_to_lower(
str_trim(
str_replace_all(
nama_bersih,
c("Kabupaten " = "", "Kota " = "")
)
)
)
) %>%
mutate(
nama_kunci = ifelse(
nama_bersih == "Kota Magelang",
"kota magelang",
ifelse(
nama_bersih == "Kabupaten Magelang",
"kabupaten magelang",
nama_kunci
)
)
)
# Cek data yang tidak cocok
tidak_cocok <- anti_join(
data.frame(nama_kunci = shp_jateng$nama_kunci),
data.frame(nama_kunci = kemiskinan$nama_kunci),
by = "nama_kunci"
)
cat("\nData yang tidak cocok:\n")
##
## Data yang tidak cocok:
print(tidak_cocok)
## [1] nama_kunci
## <0 rows> (or 0-length row.names)
# ------------------------------------------------------------------------------
# 1D. MENGGABUNGKAN DATA
# ------------------------------------------------------------------------------
cat("\nMenggabungkan data shapefile dan kemiskinan...\n")
##
## Menggabungkan data shapefile dan kemiskinan...
shp_gabung <- shp_jateng %>%
left_join(
kemiskinan %>%
select(
nama_kunci,
persen_miskin,
jumlah_miskin_ribu,
garis_kemiskinan
),
by = "nama_kunci"
)
cat("\nJumlah missing value:\n")
##
## Jumlah missing value:
print(sum(is.na(shp_gabung$persen_miskin)))
## [1] 0
# Peta kemiskinan
ggplot(shp_gabung) +
geom_sf(aes(fill = persen_miskin)) +
scale_fill_viridis_c(
name = "% Miskin",
option = "plasma"
) +
labs(
title = "Persentase Kemiskinan Jawa Tengah 2024",
subtitle = "Sumber: BPS Jawa Tengah"
) +
theme_minimal()
# ------------------------------------------------------------------------------
# 1E. DATA CURAH HUJAN NASA POWER
# ------------------------------------------------------------------------------
cat("\nMembaca data curah hujan NASA POWER...\n")
##
## Membaca data curah hujan NASA POWER...
curahhujan_raw <- read_csv(
path_curahhujan,
skip = 9
)
print(head(curahhujan_raw))
## # A tibble: 6 × 17
## PARAMETER YEAR LAT LON JAN FEB MAR APR MAY JUN JUL AUG
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 PRECTOTCORR… 2023 -6 109. 281. 294. 149. 200. 58.5 69.0 43.8 7.73
## 2 PRECTOTCORR… 2023 -6 109. 316. 303. 115. 175. 59.9 69.1 74.4 11.7
## 3 PRECTOTCORR… 2023 -6 110 312. 307. 110. 165. 58.9 68.7 81.4 12.4
## 4 PRECTOTCORR… 2023 -6 111. 314. 382. 125. 139. 64.0 65.1 116. 13.9
## 5 PRECTOTCORR… 2023 -6 111. 284. 416. 138. 140. 67.9 62.8 91.8 13.4
## 6 PRECTOTCORR… 2023 -6.5 109. 358. 392. 216. 322. 84.8 63.0 92.8 4.7
## # ℹ 5 more variables: SEP <dbl>, OCT <dbl>, NOV <dbl>, DEC <dbl>, ANN <dbl>
curahhujan <- curahhujan_raw %>%
filter(YEAR == 2023) %>%
select(
LAT,
LON,
curah_hujan_tahunan = ANN
) %>%
filter(curah_hujan_tahunan != -999)
cat("\nJumlah titik curah hujan:\n")
##
## Jumlah titik curah hujan:
print(nrow(curahhujan))
## [1] 30
print(head(curahhujan))
## # A tibble: 6 × 3
## LAT LON curah_hujan_tahunan
## <dbl> <dbl> <dbl>
## 1 -6 109. 1413.
## 2 -6 109. 1481.
## 3 -6 110 1479.
## 4 -6 111. 1595.
## 5 -6 111. 1614.
## 6 -6.5 109. 1909.
# Visualisasi titik curah hujan
ggplot(
curahhujan,
aes(
x = LON,
y = LAT,
color = curah_hujan_tahunan
)
) +
geom_point(size = 3) +
scale_color_viridis_c(
name = "CH Tahunan\n(mm)"
) +
labs(
title = "Titik Curah Hujan NASA POWER 2023",
subtitle = "Wilayah Jawa Tengah",
x = "Longitude",
y = "Latitude"
) +
theme_minimal()
# ==============================================================================
# BAB 2 - ANALISIS POLA SPASIAL
# ==============================================================================
cat("\n====================================================\n")
##
## ====================================================
cat("BAB 2 - ANALISIS POLA SPASIAL\n")
## BAB 2 - ANALISIS POLA SPASIAL
cat("====================================================\n")
## ====================================================
# Konversi ke format spatstat
centroid <- st_centroid(shp_gabung)
coords <- st_coordinates(centroid)
bbox_jateng <- st_bbox(shp_gabung)
window_jateng <- owin(
xrange = c(
bbox_jateng["xmin"],
bbox_jateng["xmax"]
),
yrange = c(
bbox_jateng["ymin"],
bbox_jateng["ymax"]
)
)
pola_titik <- ppp(
x = coords[,1],
y = coords[,2],
window = window_jateng,
marks = shp_gabung$persen_miskin
)
# Plot sebaran titik
plot(
pola_titik,
main = "Sebaran Centroid Kabupaten/Kota Jawa Tengah"
)
# Quadrat test
qt <- quadrat.test(
pola_titik,
nx = 5,
ny = 5
)
cat("\nHasil Quadrat Test:\n")
##
## Hasil Quadrat Test:
print(qt)
##
## Chi-squared test of CSR using quadrat counts
##
## data: pola_titik
## X2 = 85.22, df = 24, p-value = 1.756e-08
## alternative hypothesis: two.sided
##
## Quadrats: 5 by 5 grid of tiles
plot(
qt,
main = "Quadrat Test - Pola Spasial Jawa Tengah"
)
# Fungsi G
G_obs <- Gest(pola_titik)
plot(
G_obs,
main = "Fungsi G - KNN Analysis"
)
# Kernel Density
kde <- density(
pola_titik,
sigma = 0.5
)
plot(
kde,
main = "Kernel Density"
)
contour(kde, add = TRUE)
# Fungsi L
L_obs <- Lest(pola_titik)
plot(
L_obs,
main = "Fungsi L (Ripley's K)"
)
# ==============================================================================
# BAB 3 - MATRIKS PEMBOBOT SPASIAL
# ==============================================================================
cat("\n====================================================\n")
##
## ====================================================
cat("BAB 3 - MATRIKS PEMBOBOT SPASIAL\n")
## BAB 3 - MATRIKS PEMBOBOT SPASIAL
cat("====================================================\n")
## ====================================================
# ------------------------------------------------------------------------------
# 3A. MATRIKS CONTIGUITY - QUEEN
# ------------------------------------------------------------------------------
cat("\nMembuat matriks Queen Contiguity...\n")
##
## Membuat matriks Queen Contiguity...
nb_queen <- poly2nb(
shp_gabung,
queen = TRUE
)
summary(nb_queen)
## Neighbour list object:
## Number of regions: 41
## Number of nonzero links: 214
## Percentage nonzero weights: 12.73052
## Average number of links: 5.219512
## Link number distribution:
##
## 1 2 3 4 5 6 7 8 9 10
## 1 3 7 5 8 4 7 3 1 2
## 1 least connected region:
## 15 with 1 link
## 2 most connected regions:
## 32 33 with 10 links
W_queen <- nb2listw(
nb_queen,
style = "W",
zero.policy = TRUE
)
plot(
st_geometry(shp_gabung),
border = "grey",
main = "Matriks Queen Contiguity"
)
plot(
nb_queen,
coords,
add = TRUE,
col = "red",
lwd = 0.5
)
# ------------------------------------------------------------------------------
# 3B. MATRIKS CONTIGUITY - ROOK
# ------------------------------------------------------------------------------
cat("\nMembuat matriks Rook Contiguity...\n")
##
## Membuat matriks Rook Contiguity...
nb_rook <- poly2nb(
shp_gabung,
queen = FALSE
)
summary(nb_rook)
## Neighbour list object:
## Number of regions: 41
## Number of nonzero links: 212
## Percentage nonzero weights: 12.61154
## Average number of links: 5.170732
## Link number distribution:
##
## 1 2 3 4 5 6 7 8 9 10
## 1 4 6 5 8 4 8 2 1 2
## 1 least connected region:
## 15 with 1 link
## 2 most connected regions:
## 32 33 with 10 links
W_rook <- nb2listw(
nb_rook,
style = "W",
zero.policy = TRUE
)
plot(
st_geometry(shp_gabung),
border = "grey",
main = "Matriks Rook Contiguity"
)
plot(
nb_rook,
coords,
add = TRUE,
col = "blue",
lwd = 0.5
)
# ------------------------------------------------------------------------------
# 3C. MATRIKS BERBASIS JARAK (DISTANCE-BASED)
# ------------------------------------------------------------------------------
cat("\nMembuat matriks Distance-Based...\n")
##
## Membuat matriks Distance-Based...
# Menentukan jarak kritis
jarak_kritis <- max(
unlist(
nbdists(nb_queen, coords)
)
)
cat("\nJarak kritis:\n")
##
## Jarak kritis:
print(jarak_kritis)
## [1] 0.74543
# Membuat tetangga berdasarkan jarak
nb_jarak <- dnearneigh(
coords,
d1 = 0,
d2 = jarak_kritis
)
summary(nb_jarak)
## Neighbour list object:
## Number of regions: 41
## Number of nonzero links: 632
## Percentage nonzero weights: 37.59667
## Average number of links: 15.41463
## Link number distribution:
##
## 5 6 7 9 10 11 12 13 14 16 17 18 19 20 21
## 1 1 2 2 1 1 1 5 2 7 2 2 4 5 5
## 1 least connected region:
## 30 with 5 links
## 5 most connected regions:
## 3 18 19 39 41 with 21 links
W_jarak <- nb2listw(
nb_jarak,
style = "W",
zero.policy = TRUE
)
plot(
st_geometry(shp_gabung),
border = "grey",
main = "Matriks Distance-Based"
)
plot(
nb_jarak,
coords,
add = TRUE,
col = "darkgreen",
lwd = 0.5
)
# ------------------------------------------------------------------------------
# 3D. MATRIKS K-NEAREST NEIGHBOR (KNN)
# ------------------------------------------------------------------------------
cat("\nMembuat matriks K-Nearest Neighbor...\n")
##
## Membuat matriks K-Nearest Neighbor...
# Menggunakan 5 tetangga terdekat
nb_knn <- knn2nb(
knearneigh(coords, k = 5)
)
summary(nb_knn)
## Neighbour list object:
## Number of regions: 41
## Number of nonzero links: 205
## Percentage nonzero weights: 12.19512
## Average number of links: 5
## Non-symmetric neighbours list
## Link number distribution:
##
## 5
## 41
## 41 least connected regions:
## 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 40 41 with 5 links
## 41 most connected regions:
## 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 40 41 with 5 links
W_knn <- nb2listw(
nb_knn,
style = "W",
zero.policy = TRUE
)
plot(
st_geometry(shp_gabung),
border = "grey",
main = "Matriks KNN (k = 5)"
)
plot(
nb_knn,
coords,
add = TRUE,
col = "purple",
lwd = 0.5
)
cat("\n====================================================\n")
##
## ====================================================
cat("SEMUA MATRIKS PEMBOBOT BERHASIL DIBUAT\n")
## SEMUA MATRIKS PEMBOBOT BERHASIL DIBUAT
cat("====================================================\n")
## ====================================================
# ==============================================================================
# BAB 4 - AUTOKORELASI SPASIAL
# ==============================================================================
cat("\n====================================================\n")
##
## ====================================================
cat("BAB 4 - AUTOKORELASI SPASIAL\n")
## BAB 4 - AUTOKORELASI SPASIAL
cat("====================================================\n")
## ====================================================
y <- shp_gabung$persen_miskin
# Moran
moran_queen <- moran.test(
y,
W_queen,
zero.policy = TRUE
)
moran_rook <- moran.test(
y,
W_rook,
zero.policy = TRUE
)
moran_jarak <- moran.test(
y,
W_jarak,
zero.policy = TRUE
)
moran_knn <- moran.test(
y,
W_knn,
zero.policy = TRUE
)
cat("\n===== MORAN GLOBAL =====\n")
##
## ===== MORAN GLOBAL =====
print(moran_queen)
##
## Moran I test under randomisation
##
## data: y
## weights: W_queen
##
## Moran I statistic standard deviate = 2.5932, p-value = 0.004754
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.226084817 -0.025000000 0.009374687
print(moran_rook)
##
## Moran I test under randomisation
##
## data: y
## weights: W_rook
##
## Moran I statistic standard deviate = 2.514, p-value = 0.005968
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.220591029 -0.025000000 0.009542845
print(moran_jarak)
##
## Moran I test under randomisation
##
## data: y
## weights: W_jarak
##
## Moran I statistic standard deviate = 1.5985, p-value = 0.05497
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.050742843 -0.025000000 0.002245306
print(moran_knn)
##
## Moran I test under randomisation
##
## data: y
## weights: W_knn
##
## Moran I statistic standard deviate = 2.4638, p-value = 0.006874
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.188592660 -0.025000000 0.007515737
# Moran scatterplot
moran.plot(
y,
W_queen,
labels = shp_gabung$NAME_2,
xlab = "Persentase Kemiskinan",
ylab = "Spatial Lag",
main = "Moran Scatterplot"
)
# ------------------------------------------------------------------------------
# 4B. INDEKS GEARY GLOBAL
# ------------------------------------------------------------------------------
# Mirip Moran tetapi lebih sensitif terhadap perbedaan lokal
# Range:
# C < 1 -> clustered (autokorelasi positif)
# C > 1 -> dispersed (autokorelasi negatif)
# C = 1 -> random
cat("\n====================================================\n")
##
## ====================================================
cat("INDEKS GEARY GLOBAL\n")
## INDEKS GEARY GLOBAL
cat("====================================================\n")
## ====================================================
# Queen Contiguity
geary_queen <- geary.test(
y,
W_queen,
zero.policy = TRUE
)
# Rook Contiguity
geary_rook <- geary.test(
y,
W_rook,
zero.policy = TRUE
)
# Distance-Based
geary_jarak <- geary.test(
y,
W_jarak,
zero.policy = TRUE
)
# K-Nearest Neighbor
geary_knn <- geary.test(
y,
W_knn,
zero.policy = TRUE
)
# Menampilkan hasil
cat(
"\nQueen Contiguity :",
round(geary_queen$estimate[1], 4),
"| p-value:",
round(geary_queen$p.value, 4),
"\n"
)
##
## Queen Contiguity : 0.7408 | p-value: 0.0065
cat(
"Rook Contiguity :",
round(geary_rook$estimate[1], 4),
"| p-value:",
round(geary_rook$p.value, 4),
"\n"
)
## Rook Contiguity : 0.7475 | p-value: 0.0082
cat(
"Distance-Based :",
round(geary_jarak$estimate[1], 4),
"| p-value:",
round(geary_jarak$p.value, 4),
"\n"
)
## Distance-Based : 0.936 | p-value: 0.1024
cat(
"KNN (k = 5) :",
round(geary_knn$estimate[1], 4),
"| p-value:",
round(geary_knn$p.value, 4),
"\n"
)
## KNN (k = 5) : 0.8024 | p-value: 0.0149
cat("\n====================================================\n")
##
## ====================================================
cat("INTERPRETASI GEARY\n")
## INTERPRETASI GEARY
cat("====================================================\n")
## ====================================================
cat("C < 1 dan p-value < 0.05 -> Clustered\n")
## C < 1 dan p-value < 0.05 -> Clustered
cat("C > 1 dan p-value < 0.05 -> Dispersed\n")
## C > 1 dan p-value < 0.05 -> Dispersed
cat("C ≈ 1 -> Random\n")
## C ≈ 1 -> Random
# ------------------------------------------------------------------------------
# 4C. GETIS-ORD G GLOBAL
# ------------------------------------------------------------------------------
# Mendeteksi apakah nilai tinggi/rendah mengelompok secara global
cat("\n====================================================\n")
##
## ====================================================
cat("GETIS-ORD G GLOBAL\n")
## GETIS-ORD G GLOBAL
cat("====================================================\n")
## ====================================================
# Catatan:
# Getis-Ord menggunakan matriks non-standardisasi
# sehingga style = "B"
# Queen Binary Weight
W_queen_B <- nb2listw(
nb_queen,
style = "B",
zero.policy = TRUE
)
# Rook Binary Weight
W_rook_B <- nb2listw(
nb_rook,
style = "B",
zero.policy = TRUE
)
# Getis-Ord Queen
go_queen <- globalG.test(
y,
W_queen_B,
zero.policy = TRUE
)
# Getis-Ord Rook
go_rook <- globalG.test(
y,
W_rook_B,
zero.policy = TRUE
)
# Menampilkan hasil
cat(
"\nQueen Contiguity :",
round(go_queen$estimate[1], 4),
"| p-value:",
round(go_queen$p.value, 4),
"\n"
)
##
## Queen Contiguity : 0.1355 | p-value: 0.1995
cat(
"Rook Contiguity :",
round(go_rook$estimate[1], 4),
"| p-value:",
round(go_rook$p.value, 4),
"\n"
)
## Rook Contiguity : 0.1338 | p-value: 0.224
cat("\n====================================================\n")
##
## ====================================================
cat("INTERPRETASI GETIS-ORD G\n")
## INTERPRETASI GETIS-ORD G
cat("====================================================\n")
## ====================================================
cat("G > E(G) dan p-value < 0.05 -> Hot Spot Global\n")
## G > E(G) dan p-value < 0.05 -> Hot Spot Global
cat("G < E(G) dan p-value < 0.05 -> Cold Spot Global\n")
## G < E(G) dan p-value < 0.05 -> Cold Spot Global
# ------------------------------------------------------------------------------
# 4D. INDEKS MORAN LOKAL (LISA)
# ------------------------------------------------------------------------------
# Mendeteksi cluster dan outlier spasial di level wilayah individual
# LISA
lisa <- localmoran(
y,
W_queen,
zero.policy = TRUE
)
shp_gabung$lisa_I <- lisa[,1]
shp_gabung$lisa_pval <- lisa[,5]
mean_y <- mean(y, na.rm = TRUE)
lag_y <- lag.listw(
W_queen,
y,
zero.policy = TRUE
)
shp_gabung$lisa_kategori <- "Not Significant"
shp_gabung$lisa_kategori[
y > mean_y &
lag_y > mean_y &
lisa[,5] < 0.05
] <- "High-High"
shp_gabung$lisa_kategori[
y < mean_y &
lag_y < mean_y &
lisa[,5] < 0.05
] <- "Low-Low"
shp_gabung$lisa_kategori[
y > mean_y &
lag_y < mean_y &
lisa[,5] < 0.05
] <- "High-Low"
shp_gabung$lisa_kategori[
y < mean_y &
lag_y > mean_y &
lisa[,5] < 0.05
] <- "Low-High"
warna_lisa <- c(
"High-High" = "#d7191c",
"Low-Low" = "#2c7bb6",
"High-Low" = "#fdae61",
"Low-High" = "#abd9e9",
"Not Significant" = "#eeeeee"
)
ggplot(shp_gabung) +
geom_sf(
aes(fill = lisa_kategori),
color = "white",
size = 0.3
) +
scale_fill_manual(
values = warna_lisa,
name = "Kategori LISA"
) +
labs(
title = "Peta LISA Jawa Tengah 2024"
) +
theme_minimal()
# ------------------------------------------------------------------------------
# 4E. GETIS-ORD G* LOKAL (HOT SPOT ANALYSIS)
# ------------------------------------------------------------------------------
# Mendeteksi hot spot dan cold spot secara lokal
cat("\n====================================================\n")
##
## ====================================================
cat("GETIS-ORD G* LOKAL\n")
## GETIS-ORD G* LOKAL
cat("====================================================\n")
## ====================================================
# Matriks yang menyertakan dirinya sendiri
nb_queen_self <- include.self(nb_queen)
W_queen_self <- nb2listw(
nb_queen_self,
style = "B",
zero.policy = TRUE
)
# Menghitung Getis-Ord G* lokal
go_lokal <- localG(
y,
W_queen_self,
zero.policy = TRUE
)
# Menyimpan hasil
shp_gabung$Gstar <- as.numeric(go_lokal)
shp_gabung$Gstar_sig <- ifelse(
abs(shp_gabung$Gstar) > 1.96,
"Signifikan",
"Tidak Signifikan"
)
# Visualisasi Hot Spot / Cold Spot
ggplot(shp_gabung) +
geom_sf(
aes(fill = Gstar),
color = "white",
size = 0.3
) +
scale_fill_gradient2(
low = "#2c7bb6",
mid = "#eeeeee",
high = "#d7191c",
midpoint = 0,
name = "Z-score G*"
) +
labs(
title = "Getis-Ord G* Lokal - Hot/Cold Spot Kemiskinan",
subtitle = "Merah = Hot Spot | Biru = Cold Spot"
) +
theme_minimal()
# ------------------------------------------------------------------------------
# 4F. RANGKUMAN PERBANDINGAN SEMUA INDEKS
# ------------------------------------------------------------------------------
cat("\n====================================================\n")
##
## ====================================================
cat("RANGKUMAN PERBANDINGAN INDEKS AUTOKORELASI\n")
## RANGKUMAN PERBANDINGAN INDEKS AUTOKORELASI
cat("====================================================\n")
## ====================================================
cat(
sprintf(
"%-25s %-12s %-12s\n",
"Indeks",
"Nilai",
"p-value"
)
)
## Indeks Nilai p-value
cat(
sprintf(
"%-25s %-12s %-12s\n",
"Moran (Queen)",
round(moran_queen$estimate[1], 4),
round(moran_queen$p.value, 4)
)
)
## Moran (Queen) 0.2261 0.0048
cat(
sprintf(
"%-25s %-12s %-12s\n",
"Moran (Rook)",
round(moran_rook$estimate[1], 4),
round(moran_rook$p.value, 4)
)
)
## Moran (Rook) 0.2206 0.006
cat(
sprintf(
"%-25s %-12s %-12s\n",
"Moran (Jarak)",
round(moran_jarak$estimate[1], 4),
round(moran_jarak$p.value, 4)
)
)
## Moran (Jarak) 0.0507 0.055
cat(
sprintf(
"%-25s %-12s %-12s\n",
"Moran (KNN)",
round(moran_knn$estimate[1], 4),
round(moran_knn$p.value, 4)
)
)
## Moran (KNN) 0.1886 0.0069
cat(
sprintf(
"%-25s %-12s %-12s\n",
"Geary (Queen)",
round(geary_queen$estimate[1], 4),
round(geary_queen$p.value, 4)
)
)
## Geary (Queen) 0.7408 0.0065
cat(
sprintf(
"%-25s %-12s %-12s\n",
"Geary (Rook)",
round(geary_rook$estimate[1], 4),
round(geary_rook$p.value, 4)
)
)
## Geary (Rook) 0.7475 0.0082
cat(
sprintf(
"%-25s %-12s %-12s\n",
"Getis-Ord G (Queen)",
round(go_queen$estimate[1], 4),
round(go_queen$p.value, 4)
)
)
## Getis-Ord G (Queen) 0.1355 0.1995
cat("====================================================\n")
## ====================================================
cat("LISA : lihat peta kategori cluster lokal\n")
## LISA : lihat peta kategori cluster lokal
cat("G* Lokal : lihat peta hot spot dan cold spot\n")
## G* Lokal : lihat peta hot spot dan cold spot
cat("====================================================\n")
## ====================================================
# ==============================================================================
# BAB 5 - KRIGING
# ==============================================================================
cat("\n====================================================\n")
##
## ====================================================
cat("BAB 5 - KRIGING\n")
## BAB 5 - KRIGING
cat("====================================================\n")
## ====================================================
# Konversi spatial
coordinates(curahhujan) <- ~ LON + LAT
proj4string(curahhujan) <- CRS(
"+proj=longlat +datum=WGS84"
)
# Plot titik curah hujan
spplot(
curahhujan,
"curah_hujan_tahunan",
main = "Curah Hujan NASA POWER 2023",
col.regions = heat.colors(100)
)
# Histogram
hist(
curahhujan$curah_hujan_tahunan,
main = "Distribusi Curah Hujan",
xlab = "Curah Hujan (mm)",
col = "skyblue",
breaks = 10
)
# Statistik deskriptif
summary(curahhujan$curah_hujan_tahunan)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1145 1589 1803 1773 1988 2153
# Shapiro test
shapiro.test(
curahhujan$curah_hujan_tahunan
)
##
## Shapiro-Wilk normality test
##
## data: curahhujan$curah_hujan_tahunan
## W = 0.96146, p-value = 0.3374
# Variogram empiris
# Variogram menggambarkan bagaimana kemiripan nilai berubah seiring jarak
# Ini adalah langkah KUNCI dalam Kriging
vgm_empiris <- variogram(
curah_hujan_tahunan ~ 1,
data = curahhujan,
cutoff = 150,
width = 10
)
print(vgm_empiris)
## np dist gamma dir.hor dir.ver id
## 1 25 55.29551 77412.93 0 0 var1
## 2 24 69.01449 11501.81 0 0 var1
## 3 40 88.43538 83376.70 0 0 var1
## 4 20 110.59099 48906.98 0 0 var1
## 5 50 133.12074 39832.18 0 0 var1
## 6 30 148.69657 80599.19 0 0 var1
plot(
vgm_empiris,
main = "Variogram Empiris",
xlab = "Jarak (Km)",
ylab = "Semivariance",
pch = 19,
col = "blue"
)
#lihat outlier
# 1. Cek rentang nilai koordinat asli data Anda
bbox(curahhujan)
## min max
## LON 108.75 111.25
## LAT -8.50 -6.00
# Baca variogram: perhatikan 3 parameter penting
# - Nugget (c0) : nilai semivariance saat jarak = 0 (variabilitas mikro)
# - Sill (c0+c) : nilai semivariance saat sudah konstan (plateau)
# - Range (a) : jarak di mana variogram mencapai sill
# FITTING MODEL VARIOGRAM
# Parameter awal
nugget_awal <- min(vgm_empiris$gamma)
sill_awal <- max(vgm_empiris$gamma)
range_awal <- vgm_empiris$dist[
which.max(vgm_empiris$gamma)
] / 2
# Model spherical
model_sph <- fit.variogram(
vgm_empiris,
model = vgm(
psill = sill_awal,
model = "Sph",
range = range_awal,
nugget = nugget_awal
)
)
# Model exponential
model_exp <- fit.variogram(
vgm_empiris,
model = vgm(
psill = sill_awal,
model = "Exp",
range = range_awal,
nugget = nugget_awal
)
)
# Model gaussian
model_gau <- fit.variogram(
vgm_empiris,
model = vgm(
psill = sill_awal,
model = "Gau",
range = range_awal,
nugget = nugget_awal
)
)
cat("\n===== MODEL VARIOGRAM =====\n")
##
## ===== MODEL VARIOGRAM =====
print(model_sph)
## model psill range
## 1 Nug 11501.81 0.00000
## 2 Sph 83376.70 44.21769
print(model_exp)
## model psill range
## 1 Nug 64980.25 0.00
## 2 Exp 173285.98 1066.24
print(model_gau)
## model psill range
## 1 Nug 80008.67 0.0000
## 2 Gau 197618.44 854.5296
# SSE
sse_sph <- attr(model_sph, "SSErr")
sse_exp <- attr(model_exp, "SSErr")
sse_gau <- attr(model_gau, "SSErr")
cat("\n===== SSE =====\n")
##
## ===== SSE =====
cat("SSE Spherical :", sse_sph, "\n")
## SSE Spherical : 50481168
cat("SSE Exponential :", sse_exp, "\n")
## SSE Exponential : 28777142
cat("SSE Gaussian :", sse_gau, "\n")
## SSE Gaussian : 32310787
cat("Model Terbaik : model dengan SSE terkecil\n")
## Model Terbaik : model dengan SSE terkecil
# Model terbaik
model_terbaik <- model_sph
if (sse_exp < sse_sph &
sse_exp < sse_gau) {
model_terbaik <- model_exp
}
if (sse_gau < sse_sph &
sse_gau < sse_exp) {
model_terbaik <- model_gau
}
# Plot variogram empiris + model terbaik
plot(vgm_empiris, model_terbaik,
main = "Variogram Empiris vs Model Terpilih",
xlab = "Jarak (km)",
ylab = "Semivariance")
library(gstat)
# 1. Mengunci parameter kecocokan secara manual menggunakan fit.method = 1 atau 6
# Kita ganti taksiran awalnya agar pas dengan sebaran titik di grafik Anda
taksiran_baru <- vgm(psill = 60000, model = "Exp", range = 70, nugget = 20000)
# 2. Lakukan fitting ulang dengan metode penyeimbang bobot jarak
model_perbaikan <- fit.variogram(vgm_empiris, model = taksiran_baru, fit.method = 6)
# 3. Cetak parameter baru untuk melihat apakah range sudah realistis (< 150 km)
print(model_perbaikan)
## model psill range
## 1 Nug 45401.12 0.0000
## 2 Exp 198444.96 979.1348
# 4. Tampilkan plot hasil perbaikan
plot(vgm_empiris, model = model_perbaikan,
main = "Variogram Hasil Perbaikan Parameter",
xlab = "Jarak (km)",
ylab = "Semivariance",
col = "blue",
pch = 19)
# Grid prediksi
shp_jateng_sp <- as(
shp_gabung,
"Spatial"
)
grid_jateng <- makegrid(
shp_jateng_sp,
cellsize = 0.05
)
colnames(grid_jateng) <- c(
"LON",
"LAT"
)
coordinates(grid_jateng) <- ~ LON + LAT
proj4string(grid_jateng) <- CRS(
"+proj=longlat +datum=WGS84"
)
grid_jateng <- grid_jateng[
!is.na(
over(grid_jateng, shp_jateng_sp)[,1]
),
]
# Kriging
hasil_kriging <- krige(
formula = curah_hujan_tahunan ~ 1,
locations = curahhujan,
newdata = grid_jateng,
model = model_terbaik
)
## [using ordinary kriging]
summary(hasil_kriging)
## Object of class SpatialPointsDataFrame
## Coordinates:
## min max
## LON 108.600 111.650
## LAT -8.187 -5.837
## Is projected: FALSE
## proj4string : [+proj=longlat +datum=WGS84 +no_defs]
## Number of points: 1125
## Data attributes:
## var1.pred var1.var
## Min. :1705 Min. :76697
## 1st Qu.:1774 1st Qu.:77549
## Median :1826 Median :77883
## Mean :1825 Mean :78498
## 3rd Qu.:1874 3rd Qu.:78904
## Max. :1931 Max. :87463
# Data frame hasil
kriging_df <- as.data.frame(
hasil_kriging
)
colnames(kriging_df) <- c(
"LON",
"LAT",
"prediksi",
"varians"
)
# Peta prediksi
ggplot() +
geom_tile(
data = kriging_df,
aes(
x = LON,
y = LAT,
fill = prediksi
)
) +
geom_sf(
data = shp_gabung,
fill = NA,
color = "white",
size = 0.4
) +
scale_fill_viridis_c(
name = "CH (mm)",
option = "turbo"
) +
labs(
title = "Prediksi Curah Hujan Jawa Tengah 2023",
subtitle = paste(
"Model:",
model_terbaik$model[2]
),
x = "Longitude",
y = "Latitude"
) +
theme_minimal()
# Peta varians
ggplot() +
geom_tile(
data = kriging_df,
aes(
x = LON,
y = LAT,
fill = varians
)
) +
geom_sf(
data = shp_gabung,
fill = NA,
color = "white",
size = 0.4
) +
scale_fill_viridis_c(
name = "Varians",
option = "magma"
) +
labs(
title = "Varians Prediksi Kriging",
subtitle = "Semakin gelap = semakin tidak pasti",
x = "Longitude",
y = "Latitude"
) +
theme_minimal()
# ==============================================================================
# SELESAI
# ==============================================================================
cat("\n====================================================\n")
##
## ====================================================
cat("SEMUA ANALISIS SELESAI\n")
## SEMUA ANALISIS SELESAI
cat("====================================================\n")
## ====================================================
cat("File yang dihasilkan:\n")
## File yang dihasilkan:
cat("- jateng_kemiskinan.shp\n")
## - jateng_kemiskinan.shp
cat("- Berbagai visualisasi dan analisis spasial\n")
## - Berbagai visualisasi dan analisis spasial
cat("====================================================\n")
## ====================================================