Tugas Praktikum 2 Mata Kuliah STA2553-Analisis Statistika Spasial
Rosita Ria Rusesta | M0501251016 & Markazul Adabiyah | M0501251035
Konfigurasi Titik Ruang pada Data Titik Panas Jawa Barat
Load Libraries
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.6
## ✔ forcats 1.0.1 ✔ stringr 1.6.0
## ✔ ggplot2 4.0.1 ✔ tibble 3.3.1
## ✔ lubridate 1.9.4 ✔ tidyr 1.3.2
## ✔ purrr 1.2.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(sf)## Linking to GEOS 3.13.0, GDAL 3.8.5, PROJ 9.5.1; sf_use_s2() is TRUElibrary(rnaturalearth)
library(ggspatial)
library(spatstat)## Loading required package: spatstat.data
## Loading required package: spatstat.univar
## spatstat.univar 3.1-6
## Loading required package: spatstat.geom
## spatstat.geom 3.7-0
## Loading required package: spatstat.random
## spatstat.random 3.4-4
## Loading required package: spatstat.explore
## Loading required package: nlme
##
## Attaching package: 'nlme'
##
## The following object is masked from 'package:dplyr':
##
## collapse
##
## spatstat.explore 3.7-0
## Loading required package: spatstat.model
## Loading required package: rpart
## spatstat.model 3.6-1
## Loading required package: spatstat.linnet
## spatstat.linnet 3.4-1
##
## spatstat 3.5-1
## For an introduction to spatstat, type 'beginner'library(viridis)## Loading required package: viridisLitePersiapan Data
Data yang digunakan adalah Titik Panas di Jawa Barat per Agustus 2023 (987 titik)
df_raw <- read.csv('titik panas Jabar Agustus 2023.csv')
# Mengubah dataframe menjadi objek spasial (sf) dan ditransformasikan ke UTM Zona 48S
df_sf <- st_as_sf(df_raw, coords = c("bujur", "lintang"), crs = 4326) %>%
st_transform(crs = 32748)
# Ringkasan Data
head(df_sf)## Simple feature collection with 6 features and 10 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 664561 ymin: 9212472 xmax: 898624.6 ymax: 9275101
## Projected CRS: WGS 84 / UTM zone 48S
## id tanggal..WIB. waktu..WIB. tingkat.kepercayaan satelit
## 1 3798572 8/1/2023 0:38:58 8 snpp
## 2 3798573 8/1/2023 0:38:58 8 snpp
## 3 3798591 8/1/2023 0:38:58 8 snpp
## 4 3798599 8/1/2023 0:38:58 8 snpp
## 5 3798680 8/1/2023 1:31:36 8 noaa20
## 6 3798682 8/1/2023 1:31:36 8 noaa20
## radius.kemungkinan kecamatan kabupaten provinsi tipe
## 1 1125 Cikedung Indramayu Jawa Barat Cluster
## 2 1125 Terisi Indramayu Jawa Barat Cluster
## 3 1125 Ciawi Gerbang Kuningan Jawa Barat Cluster
## 4 1125 Ciemas Sukabumi Jawa Barat Cluster
## 5 1125 Terisi Indramayu Jawa Barat Cluster
## 6 1125 Gempol Cirebon Jawa Barat Cluster
## geometry
## 1 POINT (850685.2 9274580)
## 2 POINT (842265.7 9275101)
## 3 POINT (898624.6 9226889)
## 4 POINT (664561 9212472)
## 5 POINT (842293.7 9274989)
## 6 POINT (876022.7 9257521)class(df_sf)## [1] "sf" "data.frame"summary(df_sf)## id tanggal..WIB. waktu..WIB. tingkat.kepercayaan
## Min. :3798572 Length:987 Length:987 Min. :7.000
## 1st Qu.:3814500 Class :character Class :character 1st Qu.:8.000
## Median :3825170 Mode :character Mode :character Median :8.000
## Mean :3825381 Mean :7.891
## 3rd Qu.:3836538 3rd Qu.:8.000
## Max. :3850065 Max. :9.000
## satelit radius.kemungkinan kecamatan kabupaten
## Length:987 Min. :1125 Length:987 Length:987
## Class :character 1st Qu.:1125 Class :character Class :character
## Mode :character Median :1125 Mode :character Mode :character
## Mean :1241
## 3rd Qu.:1125
## Max. :3414
## provinsi tipe geometry
## Length:987 Length:987 POINT :987
## Class :character Class :character epsg:32748 : 0
## Mode :character Mode :character +proj=utm ...: 0
##
##
## Batas Administratif Provinsi Jawa Barat
# Mengambil dan melakukan transformasi batas administratif Provinsi Jawa Barat
jabar_poly <- ne_states(country = "indonesia", returnclass = "sf") %>%
filter(name == "Jawa Barat") %>%
st_transform(crs = 32748)
# Membuat buffer sejauh 5000 meter (5 km) di sekitar Jawa Barat agar titik di dekat batas wilayah tetap terhitung
jabar_buffered <- st_buffer(jabar_poly, dist = 5000)Pada bagian ini bertujuan untuk memperoleh dan menyiapkan batas administratif Provinsi Jawa Barat sebagai wilayah studi dalam analisis spasial.
Selain itu, dilakukan pula pembuatan buffer atau zona penyangga sejauh 5.000 meter (5 km) di sekitar batas Provinsi Jawa Barat untuk menghindari kehilangan titik pada batas wilayah.
Melakukan Konversi ke Objek Planar Point Pattern (PPP)
# Mengubah polygon Jawa Barat menjadi window pengamatan (owin)
jabar_window <- as.owin(jabar_buffered)
coords_utm <- st_coordinates(df_sf)
titikpanas <- ppp(x = coords_utm[,1], y = coords_utm[,2], window = jabar_window)
# Menghapus titik duplikat dan mengubah satuan dari meter ke km
titikpanas <- unique(titikpanas) %>% rescale(1000, "km")Pembentukan objek bertipe ppp (planar point pattern) bertujuan untuk merepresentasikan data titik panas sebagai suatu realisasi proses titik spasial dalam ruang dua dimensi. Sehingga, dimungkinkan penerapan metode statistik spasial berbasis proses titik, seperti: (1) Kernel Density Estimation; (2) Quadrant Analysis; (3) Nearest Neighbor Analisis; dan (4) Ripley’s K-function.
Penghapusan titik duplikat memastikan bahwa setiap lokasi hanya dipresentasikan satu kali dalam pola titik.
Terkait dengan
Visualisasi Sebaran Titik Panas
ggplot() +
annotation_map_tile(type = "osm", zoomin = 0) +
geom_sf(data = jabar_poly, fill = NA, color = "black", lwd = 1) +
geom_sf(data = df_sf, color = "red", size = 1.2, alpha = 0.6) +
labs(title = "Sebaran Lokasi Titik Panas Jabar (Buffer 5km)") +
theme_minimal()## Zoom: 9
Analisis Intensitas (Heatmap)
# Menghitung bandwidth optimal
sigma_opt <- bw.ppl(titikpanas)
# Menghitung estimasi kepadatan titik
ds <- density(titikpanas, sigma = sigma_opt)
# Menampilkan peta, garis kontur dan titik asli
plot(ds, main = "Peta Intensitas Hotspot (Titik/Km²)", col = magma(100), las = 1)
contour(ds, add = TRUE, col = "#FFFFFF80", drawlabels = FALSE)
plot(titikpanas, add = TRUE, pch = 20, cex = 0.2, col = rgb(1,1,1,0.3))
#Menampilkan Kernel Titik Panas
kernel_titikpanas <- density(titikpanas)
plot(kernel_titikpanas)
#Menampilkan Countour
contour(kernel_titikpanas)
Metode Kuadran
# Membagi area menjadi 5x3 grid dan menghitung jumlah titik per kuadran
q_count <- quadratcount(titikpanas, nx = 5, ny = 3)
plot(q_count, main = "Quadrat Count: Jumlah Titik per Kuadran", col = "red")
plot(jabar_window, add = TRUE, border = "blue", lty = 2)
# Menghitung VMR
counts <- as.vector(q_count)
vmr_val <- var(counts) / mean(counts)
cat("\nNilai VMR:", round(vmr_val, 3), "\n")##
## Nilai VMR: 136.014Keterangan :
VMR = 0 titik menyebar Uniform(sistematik)
VMR = 1 titik menyebar acak
VMR > 1 titik menyebar lebih mengelompok
Nilai VMR =136.014 —> titik menyebar lebih mengelompok
Oleh karena itu, dilakukan uji hipotesis menggunakan statistik uji Chisquare dengan Hipotesis sebagai berikut:
H0: Titik menyebar acak
H1: Titik mempunyai pola mengelompok (clustered)
# Uji Chi-square
quadrat.test(q_count, alternative = "clustered")##
## Chi-squared test of CSR using quadrat counts
##
## data:
## X2 = 1061.7, df = 14, p-value < 2.2e-16
## alternative hypothesis: clustered
##
## Quadrats: 15 tiles (irregular windows)Hasil uji hipotesis menghasilkan nilai p-value<2.2e-16 (<0.05) yang berarti Tolak H0. Maka dai itu, cukup bukti untuk menyatakan bahwa pola titik berkelompok dengan taraf nyata 5%.
Empirical K-Function
# Menghitung dan menampilkan histogram jarak tetangga terdekat
nn_titikpanas <- nndist(titikpanas)
hist(nn_titikpanas)
Interpretasi Histogram :
nn_titikpanas merupakan jarak dari setiap titik panas ke tetangga terdekatnya atau menggambarkan seberapa dekat hotspot dengan hotspot lainnya.
Banyak jarak yang sangat kecil
Banyak hotspot memiliki tetangga sangat dekat
Mengindikasikan hotspot sering muncul berdekatan
Ini dapat mengindikasikan adanya hotspot yang cenderung mengelompok (clustered), bukan tersebar merata.
Distribusi menjulur kanan
Ada sedikit hotspot dengan jarak sangat jauh
Artinya ada hotspot yang relatif terisolasi
# Menghitung dan menampilkan fungsi K Ripley
K<- Kest(titikpanas, correction="Ripley")
plot(K)
Berdasarkan hasil visualisasi, terlihat bahwa nilai fungsi K empiris secara konsisten berada di atas fungsi K teoritis pada seluruh rentang jarak pengamatan. Kondisi ini menunjukkan bahwa jumlah rata-rata titik tetangga dalam radius tertentu lebih besar dibandingkan dengan pola distribusi acak. Dengan demikian, pola persebaran titik panas menunjukkan kecenderungan pengelompokan spasial (clustered pattern).
Fungsi K empiris adalah: Rata-rata kumulatif jumlah titik data yang berada dalam jarak r dari suatu titik data acuan, yang telah dikoreksi terhadap efek tepi, dan dinormalisasi dengan membaginya terhadap intensitas titik.
E<-envelope(titikpanas, Kest, nsim=99)## Generating 99 simulations of CSR ...
## 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, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
## 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
## 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98,
## 99.
##
## Done.plot(E)
Berdasarkan grafik, kurva fungsi K empiris secara konsisten berada di atas kurva fungsi K teoritis CSR dan juga berada di luar batas atas envelope simulasi pada hampir seluruh rentang jarak pengamatan. Kondisi ini menunjukkan bahwa jumlah rata-rata titik tetangga dalam radius tertentu secara signifikan lebih besar dibandingkan dengan yang diharapkan pada pola acak.
Secara statistik, keberadaan kurva empiris di luar envelope CSR mengindikasikan bahwa pola persebaran titik tidak mengikuti Complete Spatial Randomness. Sebaliknya, pola titik menunjukkan kecenderungan pengelompokan spasial (clustered pattern) yang signifikan.
# Uji statistik non-grafis
mad.test(titikpanas, Kest, nsim=99, alternative="two.sided")## Generating 99 simulations of CSR ...
## 1, 2, [3:36 remaining] 3,
## [3:27 remaining] 4, [3:14 remaining] 5, [3:06 remaining] 6,
## [3:01 remaining] 7, [3:01 remaining] 8, [2:58 remaining] 9,
## [2:58 remaining] 10, [2:59 remaining] 11, [2:54 remaining] 12,
## [2:52 remaining] 13, [2:50 remaining] 14, [2:49 remaining] 15,
## [2:47 remaining] 16, [2:44 remaining] 17, [2:42 remaining] 18,
## [2:40 remaining] 19, [2:38 remaining] 20, [2:35 remaining] 21,
## [2:34 remaining] 22, [2:31 remaining] 23, [2:30 remaining] 24,
## [2:29 remaining] 25, [2:26 remaining] 26, [2:24 remaining] 27,
## [2:22 remaining] 28, [2:20 remaining] 29, [2:18 remaining] 30,
## [2:16 remaining] 31, [2:14 remaining] 32, [2:12 remaining] 33,
## [2:10 remaining] 34, [2:08 remaining] 35, [2:06 remaining] 36,
## [2:04 remaining] 37, [2:03 remaining] 38, [2:00 remaining] 39,
## [1:58 remaining] 40, [1:56 remaining] 41, [1:54 remaining] 42,
## [1:52 remaining] 43, [1:50 remaining] 44, [1:48 remaining] 45,
## [1:46 remaining] 46, [1:44 remaining] 47, [1:42 remaining] 48,
## [1:40 remaining] 49, [1:38 remaining] 50, [1:36 remaining] 51,
## [1:34 remaining] 52, [1:32 remaining] 53, [1:30 remaining] 54,
## [1:28 remaining] 55, [1:26 remaining] 56, [1:24 remaining] 57,
## [1:22 remaining] 58, [1:21 remaining] 59, [1:19 remaining] 60,
## [1:17 remaining] 61, [1:15 remaining] 62, [1:13 remaining] 63,
## [1:11 remaining] 64, [1:09 remaining] 65, [1:07 remaining] 66,
## [1:05 remaining] 67, [1:03 remaining] 68, [1:01 remaining] 69,
## [59 sec remaining] 70, [57 sec remaining] 71, [55 sec remaining] 72,
## [53 sec remaining] 73, [51 sec remaining] 74, [49 sec remaining] 75,
## [47 sec remaining] 76, [45 sec remaining] 77, [43 sec remaining] 78,
## [41 sec remaining] 79, [39 sec remaining] 80, [37 sec remaining] 81,
## [35 sec remaining] 82, [34 sec remaining] 83, [32 sec remaining] 84,
## [30 sec remaining] 85, [28 sec remaining] 86, [26 sec remaining] 87,
## [24 sec remaining] 88, [22 sec remaining] 89, [20 sec remaining] 90,
## [18 sec remaining] 91, [16 sec remaining] 92, [14 sec remaining] 93,
## [12 sec remaining] 94, [10 sec remaining] 95, [8 sec remaining] 96,
## [6 sec remaining] 97, [4 sec remaining] 98, [2 sec remaining]
## 99.
##
## Done.##
## Maximum absolute deviation test of CSR
## Monte Carlo test based on 99 simulations
## Summary function: K(r)
## Reference function: theoretical
## Alternative: two.sided
## Interval of distance values: [0, 55.4114151592557] km
## Test statistic: Maximum absolute deviation
## Deviation = observed minus theoretical
##
## data: titikpanas
## mad = 5274.9, rank = 1, p-value = 0.01Interpretasi Hasil Uji CSR
Hasil uji Monte Carlo terhadap fungsi K Ripley menunjukkan bahwa nilai statistik uji Maximum Absolute Deviation (MAD) sebesar 5274,9 dengan nilai p-value sebesar 0,01 berdasarkan 99 simulasi.
Uji ini dilakukan untuk menguji hipotesis:
H₀ : Pola persebaran titik panas mengikuti Complete Spatial Randomness (CSR).
H₁ : Pola persebaran titik panas tidak mengikuti Complete Spatial Randomness.
Berdasarkan hasil pengujian, diperoleh nilai p-value sebesar 0,01 yang lebih kecil dari tingkat signifikansi umum (α = 0,05). Oleh karena itu, Tolak Hipotesis nol. Hal ini menunjukkan bahwa terdapat cukup bukti untuk menyatakan Pola persebaran titik panas tidak mengikuti Complete Spatial Randomness.
Selain itu, nilai MAD yang besar menunjukkan adanya deviasi maksimum antara fungsi K empiris dengan fungsi K teoritis CSR dalam rentang jarak pengamatan antara 0 hingga 55,41 km. Nilai deviasi tersebut mengindikasikan bahwa pola persebaran titik panas memiliki struktur spasial yang tidak acak.