1 Pendahuluan

Modul ini merupakan pendamping praktikum dari materi kuliah “Analisis Eksplorasi Data Point Pattern”. Kita akan mempraktikkan konsep-konsep berikut menggunakan R, khususnya paket spatstat:

  1. Exploratory Spatial Data Analysis (ESDA)
  2. Data Point Pattern
  3. Intensitas (rata-rata, bervariasi spasial, relatif, terboboti)
  4. Complete Spatial Randomness (CSR): uji kuadrat & fungsi-K
  5. Kernel Density Estimator (KDE)
  6. Hot spots, clusters, & local features

Kita akan menggunakan dua jenis data:

  • Data bawaan spatstat: swedishpines dan bei — untuk mempelajari konsep dasar pada pola titik yang (relatif) homogen dan inhomogen.
  • Studi kasus nyata: chorley — data lokasi kasus kanker laring dan kanker paru-paru di wilayah Chorley-Ribble, Inggris, beserta lokasi insinerator sampah lama. Data ini sering dipakai sebagai contoh klasik analisis hot spot dan rasio intensitas kasus penyakit terhadap populasi kontrol.

1.1 Paket yang Dibutuhkan

# install.packages("spatstat")  # jalankan sekali saja jika belum terpasang
# install.packages("sf")        # dibutuhkan di Bagian 7 untuk mengimpor data asli
library(spatstat)
library(sf)

2 Data Point Pattern

Ingat kembali: data point pattern terbentuk ketika peubah yang dianalisis berhubungan dengan lokasi kejadian, dan lokasi tersebut dianggap acak (bukan lokasi yang sudah ditentukan sebelumnya seperti pada data geostatistik).

Di R, objek pola titik direpresentasikan dengan kelas ppp (planar point pattern) dari paket spatstat, yang menyimpan:

  • koordinat titik (x, y)
  • jendela pengamatan/window (W)
  • (opsional) atribut/label pada tiap titik (marks)

2.1 Contoh 1: Data swedishpines

Data ini berisi lokasi 71 pohon pinus muda di suatu petak hutan di Swedia. contoh klasik pola titik yang relatif homogen karena kepadatan titiknya dianggap kurang lebih sama di seluruh area.

data(swedishpines)
X <- swedishpines
X
## Planar point pattern: 71 points
## window: rectangle = [0, 96] x [0, 100] units (one unit = 0.1 metres)
summary(X)
## Planar point pattern:  71 points
## Average intensity 0.007395833 points per square unit (one unit = 0.1 metres)
## 
## Coordinates are integers
## i.e. rounded to the nearest unit (one unit = 0.1 metres)
## 
## Window: rectangle = [0, 96] x [0, 100] units
## Window area = 9600 square units
## Unit of length: 0.1 metres
plot(X, main = "Lokasi Pohon Pinus Swedia (swedishpines)", pch = 16, cols = "darkgreen")

Interpretasi: Data swedishpines berisi 71 lokasi pohon pinus dalam satu area pengamatan berbentuk persegi panjang. Karena tidak ada marks, data ini disebut unmarked point pattern, artinya yang dianalisis hanya letak pohonnya, bukan karakteristik pohonnya.

2.2 Contoh 2: Studi Kasus Nyata — Data chorley

Data chorley berisi lokasi kasus kanker laring dan kanker paru-paru di wilayah Chorley-Ribble, Inggris, serta lokasi bekas insinerator. Data ini termasuk marked point pattern karena setiap titik punya label/kategori, yaitu jenis kasus kanker.

data(chorley)
chorley
## Marked planar point pattern: 1036 points
## Multitype, with levels = larynx, lung 
## window: polygonal boundary
## enclosing rectangle: [343.45, 366.45] x [410.41, 431.79] km
summary(chorley)
## Marked planar point pattern:  1036 points
## Average intensity 3.287268 points per square km
## 
## *Pattern contains duplicated points*
## 
## Coordinates are given to 1 decimal place
## i.e. rounded to the nearest multiple of 0.1 km
## 
## Multitype:
##        frequency proportion intensity
## larynx        58 0.05598456 0.1840363
## lung         978 0.94401540 3.1032320
## 
## Window: polygonal boundary
## single connected closed polygon with 131 vertices
## enclosing rectangle: [343.45, 366.45] x [410.41, 431.79] km
##                      (23 x 21.38 km)
## Window area = 315.155 square km
## Unit of length: 1 km
## Fraction of frame area: 0.641
plot(chorley, main = "Kasus Kanker Laring & Paru-paru (Chorley-Ribble)", cols = c("red", "grey60"))
points(chorley.extra$incin, pch = 17, col = "blue", cex = 1.5)
legend("topleft", legend = c("Kanker laring (kasus)", "Kanker paru-paru (kontrol)", "Insinerator"),
       col = c("red", "grey60", "blue"), pch = c(1, 1, 17), bty = "n", cex = 0.8)

Interpretasi: Setiap titik menunjukkan lokasi kasus kanker, dengan jenis kasus yang berbeda. Titik merah adalah kanker laring, titik abu-abu adalah kanker paru-paru sebagai pembanding, dan segitiga biru menunjukkan lokasi insinerator. Selanjutnya, kita ingin melihat apakah kasus kanker laring cenderung lebih banyak ditemukan di sekitar insinerator dibandingkan kanker paru-paru.

Kita bisa memisahkan menjadi dua pola titik terpisah bila diperlukan:

cases    <- split(chorley)$larynx
controls <- split(chorley)$lung
plot(cases, main = "Kasus Kanker Laring", pch = 16, cols = "red")

Interpretasi: Setelah data dipisahkan, terlihat sebaran lokasi kasus kanker laring saja. Titik-titik merah menunjukkan lokasi kasus dan tampak tersebar di beberapa bagian wilayah, dengan beberapa titik yang cukup berdekatan.

3 Intensitas

Intensitas pada data point pattern menunjukkan tingkat kepadatan titik dalam suatu wilayah pengamatan. Intensitas dihitung dengan membagi jumlah titik dengan luas wilayah atau window. Semakin besar nilai intensitas, semakin padat titik-titik yang terdapat pada wilayah tersebut. Sebaliknya, semakin kecil nilai intensitas, semakin jarang titik yang tersebar.

Secara sederhana, intensitas dihitung dengan rumus:

\[ \hat{\lambda} = \frac{n}{|W|} \]

dengan:

  • \(n\) = jumlah titik
  • \(|W|\) = luas wilayah pengamatan atau window

3.1 Rata-rata Intensitas

Intensitas \(\hat{\lambda}\) untuk pola titik homogen dihitung sederhana sebagai banyaknya titik dibagi luas window.

intensity(X)  # intensitas swedishpines (titik per satuan luas)
## [1] 0.007395833

Pada data swedishpines, terdapat 71 titik pohon dalam window seluas 9600 satuan luas, sehingga diperoleh intensitas sekitar 0,0074 titik per satuan luas. Nilai ini menggambarkan rata-rata kepadatan pohon pada seluruh wilayah pengamatan.

# Cek manual
n_titik <- npoints(X)
luas    <- area(Window(X))
n_titik / luas
## [1] 0.007395833

3.2 Intensitas dengan Metode Kuadrat (quadratcount)

Metode kuadrat membagi window menjadi \(r \times c\) subwilayah untuk melihat sebaran jumlah titik antarwilayah.

Q <- quadratcount(X, nx = 4, ny = 3)
Q
##              x
## y             [0,24] (24,48] (48,72] (72,96]
##   (66.7,100]       7       3       6       5
##   (33.3,66.7]      5       9       7       7
##   [0,33.3]         4       3       6       9
plot(X, pch = 16, cols = "darkgreen", main = "Quadrat Count: swedishpines")
plot(Q, add = TRUE, cex = 1.5, col = "red")

Interpretasi: Area pengamatan dibagi menjadi beberapa kotak, lalu dihitung jumlah titik pada setiap kotak. Jika pola titik benar-benar acak dan relatif merata, jumlah titik di tiap kotak seharusnya tidak berbeda terlalu jauh. Pada gambar ini, jumlah titik tiap kotak masih bervariasi, misalnya ada kotak berisi 3 titik dan ada yang 9 titik. Perbedaan ini menunjukkan bahwa kepadatan titik belum tentu sama di seluruh area, sehingga perlu diuji lebih lanjut apakah pola titik benar-benar acak atau tidak.

3.3 Intensitas yang Bervariasi Secara Spasial

Data bei berisi lokasi pohon Beilschmiedia di hutan tropis Panama. Data ini merupakan contoh pola titik dengan intensitas inhomogen, artinya kepadatan titik tidak sama di seluruh wilayah pengamatan.

Perbedaan kepadatan tersebut dapat dipengaruhi oleh kondisi lingkungan, seperti elevasi dan kemiringan lahan. Jadi, pada beberapa bagian wilayah pohon dapat tumbuh lebih rapat, sedangkan pada bagian lain lebih jarang.

data(bei)
plot(bei, pch = ".", main = "Lokasi Pohon Beilschmiedia (bei)")

plot(bei.extra$elev, main = "Kovariat Elevasi")
points(bei, pch = ".", col = "white")

Interpretasi: Warna pada peta menunjukkan perbedaan elevasi wilayah. Area dengan warna lebih terang memiliki elevasi yang lebih tinggi, sedangkan area dengan warna lebih gelap memiliki elevasi yang lebih rendah. Peta ini digunakan untuk melihat apakah persebaran dan kepadatan pohon bei berkaitan dengan kondisi elevasi di wilayah tersebut.

Kita bisa menduga intensitas sebagai fungsi dari kovariat elevasi dengan rhohat():

rho <- rhohat(bei, bei.extra$elev)
plot(rho, main = "Intensitas sebagai Fungsi Elevasi")

Interpretasi: Kurva menunjukkan bahwa intensitas pohon berubah mengikuti nilai elevasi. Karena kurva tidak datar, kepadatan pohon tidak sama pada setiap tingkat elevasi. Artinya, pola titik bei bersifat inhomogen dan elevasi diduga memengaruhi kepadatan pohon.

3.4 Intensitas Relatif (Rasio Intensitas): Studi Kasus chorley

Untuk membandingkan distribusi spasial dua jenis titik (kasus vs kontrol), kita hitung rasio intensitas dengan relrisk().

risk_ratio <- relrisk(chorley, sigma = 0.15, casecontrol = TRUE)
plot(risk_ratio, main = "Rasio Intensitas: Kanker Laring vs Kanker Paru-paru")
points(chorley.extra$incin, pch = 17, col = "white", cex = 1.5)

Interpretasi: Warna yang lebih terang menunjukkan wilayah dengan proporsi kasus kanker laring yang lebih tinggi dibandingkan kanker paru-paru sebagai pembanding. Selanjutnya, pola ini dapat dibandingkan dengan lokasi insinerator untuk melihat apakah wilayah dengan proporsi kanker laring yang tinggi cenderung berada di sekitar insinerator.

3.5 Intensitas Terboboti

Pada beberapa data point pattern, setiap titik bisa memiliki bobot yang berbeda. Misalnya, satu lokasi mewakili lebih dari satu kasus, atau setiap titik memiliki nilai tertentu seperti massa, jumlah kejadian, atau ukuran lainnya.

Jika setiap titik memiliki bobot, maka perhitungan kepadatan tidak hanya melihat jumlah titik, tetapi juga mempertimbangkan besar bobot pada masing-masing titik.

4 Complete Spatial Randomness (CSR)

Complete Spatial Randomness (CSR) adalah kondisi ketika titik-titik tersebar secara acak di seluruh wilayah pengamatan.

Pada kondisi CSR:

  • setiap titik muncul secara independen,
  • tidak ada kecenderungan titik untuk saling mengelompok,
  • dan tidak ada kecenderungan titik untuk saling menjauh.

Selanjutnya, kita akan menguji apakah pola titik pada data swedishpines dapat dianggap mengikuti pola CSR atau tidak.

4.1 Uji Quadrat (Chi-Square Test)

Uji Quadrat digunakan untuk mengecek apakah pola titik tersebar secara acak atau tidak.

Hipotesis:

  • \(H_0\): pola titik mengikuti CSR atau tersebar secara acak.
  • \(H_1\): pola titik tidak mengikuti CSR atau tidak tersebar secara acak.
Q <- quadratcount(X, nx = 4, ny = 3)
uji_csr <- quadrat.test(X, nx = 4, ny = 3)
uji_csr
## 
##  Chi-squared test of CSR using quadrat counts
## 
## data:  X
## X2 = 7.5915, df = 11, p-value = 0.5013
## alternative hypothesis: two.sided
## 
## Quadrats: 4 by 3 grid of tiles

Interpretasi: Hasil uji menunjukkan p-value sebesar 0,5013. Karena p-value lebih besar dari 0,05, maka gagal menolak \(H_0\). Artinya, belum ada bukti yang cukup bahwa pola titik swedishpines menyimpang dari pola acak. Jadi, pola titik tersebut masih dapat dianggap mengikuti Complete Spatial Randomness (CSR).

4.2 Pendekatan Fungsi-K (Ripley’s K)

Fungsi-K digunakan untuk melihat pola hubungan antartitik pada berbagai jarak. Metode ini membantu mengetahui apakah titik-titik cenderung acak, mengelompok, atau saling berjauhan.

Berbeda dengan uji kuadrat, Fungsi-K tidak bergantung pada pembagian wilayah menjadi kotak-kotak tertentu, sehingga dapat memberikan gambaran pola titik yang lebih detail pada berbagai skala jarak.

K_obs <- Kest(X)
plot(K_obs, main = "Fungsi-K Ripley: swedishpines")

Interpretasi:

  • Jika kurva hasil observasi mendekati kurva teoritis Poisson, maka pola titik cenderung acak.
  • Jika kurva observasi berada di atas kurva teoritis, maka titik cenderung mengelompok (clustered).
  • Jika kurva observasi berada di bawah kurva teoritis, maka titik cenderung teratur atau saling berjauhan (regular).

Pada grafik swedishpines, kurva observasi terlihat cukup dekat dengan kurva teoritis, sehingga pola titik secara umum masih cenderung mendekati pola acak.

4.3 Uji Formal dengan Envelope (Simulasi Monte Carlo)

Agar kesimpulan dari Fungsi-K tidak hanya berdasarkan pengamatan visual, kita dapat menggunakan envelope.

Envelope adalah batas yang diperoleh dari beberapa kali simulasi pola acak (CSR). Kurva observasi kemudian dibandingkan dengan batas tersebut.

  • Jika kurva observasi masih berada di dalam envelope, pola titik masih dapat dianggap acak (CSR).
  • Jika kurva observasi keluar dari envelope, ada indikasi bahwa pola titik tidak acak, misalnya mengelompok atau teratur.
set.seed(123)
env_csr <- envelope(X, Kest, nsim = 99, verbose = FALSE)
plot(env_csr, main = "Uji CSR dengan Envelope Fungsi-K")

Interpretasi: Garis hitam menunjukkan pola yang diamati, sedangkan area abu-abu merupakan envelope atau batas pola acak hasil simulasi CSR.

Pada grafik ini, garis observasi masih berada di dalam area envelope untuk hampir seluruh jarak \(r\). Artinya, pola titik swedishpines masih konsisten dengan CSR atau pola acak.

Jika garis observasi keluar dari batas atas envelope, ada indikasi titik mengelompok. Jika keluar dari batas bawah, ada indikasi titik lebih teratur atau saling berjauhan.

5 Kernel Density Estimator (KDE)

Kernel Density Estimator (KDE) digunakan untuk melihat bagaimana kepadatan titik berubah di seluruh wilayah pengamatan secara lebih halus.

Berbeda dengan metode kuadrat yang membagi wilayah menjadi beberapa kotak, KDE menghasilkan permukaan kepadatan yang lebih halus sehingga area dengan titik yang lebih padat dan lebih jarang dapat terlihat dengan lebih jelas.

5.1 Pengaruh Bandwidth (\(h\))

Bandwidth mengatur tingkat kehalusan hasil KDE.

  • Bandwidth terlalu kecil → hasil terlihat terlalu detail dan berbintik-bintik (undersmoothed).
  • Bandwidth sedang → pola kepadatan lebih mudah terlihat.
  • Bandwidth terlalu besar → hasil terlalu halus sehingga pola lokal bisa hilang (oversmoothed).
par(mfrow = c(2, 2))
plot(density(X, sigma = 0.05), main = "sigma = 0.05 (undersmoothed)")
plot(density(X, sigma = 0.1),  main = "sigma = 0.1")
plot(density(X),               main = "Bandwidth Default")
plot(density(X, sigma = 1),    main = "sigma = 1 (oversmoothed)")

par(mfrow = c(1, 1))

Interpretasi: Nilai bandwidth menentukan seberapa halus peta kepadatan yang dihasilkan.

  • Bandwidth yang terlalu kecil membuat peta terlalu detail sehingga pola utama sulit terlihat (undersmoothed).
  • Bandwidth yang terlalu besar membuat peta terlalu halus sehingga variasi atau pola lokal bisa hilang (oversmoothed).
  • Karena itu, diperlukan bandwidth yang sesuai agar pola kepadatan titik dapat terlihat dengan jelas.

5.2 Pemilihan Bandwidth Otomatis

Agar bandwidth tidak dipilih secara sembarangan, spatstat menyediakan beberapa metode untuk menentukan nilai \(\sigma\) secara otomatis berdasarkan data.

Beberapa metode yang dapat digunakan antara lain:

  • bw.diggle() → memilih bandwidth dengan metode cross-validation.
  • bw.ppl() → memilih bandwidth dengan metode likelihood cross-validation.

Tujuannya adalah memperoleh bandwidth yang paling sesuai agar peta KDE tidak terlalu kasar maupun terlalu halus.

bw_opt <- bw.diggle(X)
bw_opt
##    sigma 
## 5.706458
plot(density(X, sigma = bw_opt), main = paste("KDE dengan Bandwidth Optimal (sigma =", round(bw_opt, 3), ")"))
plot(X, add = TRUE, pch = 16, cex = 0.5, col = "white")

Interpretasi: Nilai bandwidth optimal yang diperoleh adalah 5,706. Dengan bandwidth ini, KDE menghasilkan peta kepadatan yang lebih seimbang, yaitu tidak terlalu kasar dan tidak terlalu halus.

Area dengan warna lebih terang menunjukkan kepadatan titik yang lebih tinggi, sedangkan area yang lebih gelap menunjukkan kepadatan titik yang lebih rendah. Titik putih menunjukkan lokasi pohon swedishpines.

5.3 KDE untuk Studi Kasus chorley

dens_cases <- density(cases, sigma = bw.diggle(cases))
plot(dens_cases, main = "Kepadatan Kasus Kanker Laring (KDE)")
points(chorley.extra$incin, pch = 17, col = "white", cex = 1.5)

Interpretasi: Warna yang semakin terang menunjukkan kepadatan kasus kanker laring yang semakin tinggi, sedangkan warna yang lebih gelap menunjukkan kepadatan yang lebih rendah.

Area terang dapat dianggap sebagai lokasi dengan konsentrasi kasus yang lebih tinggi. Selanjutnya, posisi area tersebut dapat dibandingkan dengan lokasi insinerator untuk melihat apakah konsentrasi kasus cenderung berada di dekat insinerator.

6 Hot Spots, Clusters, & Local Features

Hot spot adalah wilayah yang memiliki kepadatan atau intensitas titik lebih tinggi dibandingkan wilayah di sekitarnya.

Salah satu cara paling sederhana untuk melihat hot spot adalah menggunakan peta KDE. Pada peta KDE, area dengan warna lebih terang menunjukkan kepadatan titik yang lebih tinggi.

6.1 Identifikasi Hot Spot secara Visual dari KDE

dens_bei <- density(bei, sigma = bw.diggle(bei))
plot(dens_bei, main = "Peta Intensitas (Hot Spot) — Data bei")
contour(dens_bei, add = TRUE, col = "white", nlevels = 5)

Interpretasi: Pada peta data bei, area dengan warna lebih terang menunjukkan lokasi dengan konsentrasi pohon yang lebih tinggi. Area tersebut dapat dianggap sebagai hot spot secara visual. Namun, kesimpulan ini masih berdasarkan peta KDE, sehingga jika ingin lebih formal perlu dilakukan analisis lanjutan.

6.2 Hot Spot Berbasis Rasio Intensitas: chorley

Pada data chorley, hot spot tidak hanya dilihat dari banyaknya kasus kanker laring, tetapi dari risiko relatifnya dibandingkan kanker paru-paru sebagai pembanding.

Dengan cara ini, kita bisa melihat wilayah mana yang memiliki proporsi kasus kanker laring relatif lebih tinggi.

plot(risk_ratio, main = "Hot Spot Risiko Relatif Kanker Laring")
contour(risk_ratio, add = TRUE, nlevels = 4)
points(chorley.extra$incin, pch = 17, col = "white", cex = 1.5)

Interpretasi: Warna yang semakin terang menunjukkan wilayah dengan risiko relatif kanker laring yang lebih tinggi dibandingkan kanker paru-paru sebagai pembanding.

Area yang paling terang dapat dianggap sebagai hot spot risiko relatif. Lokasi hot spot ini kemudian dibandingkan dengan posisi insinerator untuk melihat apakah wilayah dengan risiko relatif tinggi cenderung berada di sekitarnya.

6.3 Uji Signifikansi Local Clustering

Fungsi-L merupakan transformasi dari Fungsi-K yang dibuat agar pola pengelompokan lebih mudah dibaca.

Secara umum:

  • Jika \(L(r)-r = 0\), pola titik cenderung acak.
  • Jika \(L(r)-r > 0\), pola titik cenderung mengelompok (clustered).
  • Jika \(L(r)-r < 0\), pola titik cenderung teratur atau saling berjauhan.
L_bei <- Lest(bei, correction = "Ripley")
plot(L_bei, . - r ~ r, main = "Fungsi-L Terpusatkan: Data bei")

Interpretasi: Pada grafik data bei, kurva \(L(r)-r\) berada di atas garis nol hampir di seluruh jarak \(r\). Hal ini menunjukkan bahwa titik-titik pohon bei cenderung mengelompok, bukan tersebar secara acak.

7 Studi Kasus Tambahan: Impor Data Asli Langsung dari Internet

Selain menggunakan data bawaan spatstat, kita juga bisa menggunakan data point pattern asli yang tersedia di internet dan langsung membacanya melalui R.

Pada contoh ini digunakan data gempa bumi dari USGS Earthquake Catalog. Data tersebut bersifat publik, gratis, dan terus diperbarui. Analisis difokuskan pada wilayah Jawa dan sekitarnya, karena wilayah ini memiliki aktivitas gempa yang cukup tinggi.

Catatan: Karena data gempa terus diperbarui (real-time), jumlah titik dan hasil analisis yang diperoleh bisa berbeda setiap kali kode dijalankan.

7.1 Mengimpor Data dari API

R dapat langsung membaca data dari internet menggunakan read.csv() jika sumber data menyediakan file dalam format CSV.

Pada contoh ini, data gempa diambil dari USGS Earthquake Catalog melalui API. URL yang dibuat berisi beberapa batasan, seperti:

  • periode waktu gempa,
  • batas lintang dan bujur wilayah,
  • magnitudo minimum,
  • serta urutan data berdasarkan waktu.
url_gempa <- paste0(
  "https://earthquake.usgs.gov/fdsnws/event/1/query?",
  "format=csv",
  "&starttime=2018-01-01&endtime=", Sys.Date(),
  "&minlatitude=-9&maxlatitude=-5.5",
  "&minlongitude=108&maxlongitude=114",
  "&minmagnitude=4.0",
  "&orderby=time"
)

gempa_df <- read.csv(url_gempa)
str(gempa_df)
## 'data.frame':    306 obs. of  22 variables:
##  $ time           : chr  "2026-08-07T14:07:39.335Z" "2026-07-28T13:29:48.071Z" "2026-07-28T09:21:24.654Z" "2026-07-24T05:00:13.219Z" ...
##  $ latitude       : num  -5.6 -8.97 -7.48 -7.76 -7.21 ...
##  $ longitude      : num  110 113 108 112 112 ...
##  $ depth          : num  535 62.5 155.5 163.6 220 ...
##  $ mag            : num  4.8 4.4 4.9 4.2 4.3 5.2 4.6 4.9 4.4 4 ...
##  $ magType        : chr  "mb" "mb" "mb" "mb" ...
##  $ nst            : int  76 31 113 23 32 114 75 58 32 23 ...
##  $ gap            : int  47 58 35 103 100 55 62 48 106 75 ...
##  $ dmin           : num  1.461 1.403 0.713 1.646 1.499 ...
##  $ rms            : num  0.77 0.56 0.64 1.06 0.91 0.77 0.65 0.87 0.42 1.04 ...
##  $ net            : chr  "us" "us" "us" "us" ...
##  $ id             : chr  "us6000tixn" "us6000tha9" "us6000tgbx" "us6000tgpx" ...
##  $ updated        : chr  "2026-08-22T20:48:06.040Z" "2026-08-20T20:38:30.040Z" "2026-08-20T08:09:16.537Z" "2026-08-18T13:06:39.040Z" ...
##  $ place          : chr  "112 km NNW of Batang, Indonesia" "83 km S of Dampit, Indonesia" "16 km S of Ciamis, Indonesia" "8 km W of Gampengrejo, Indonesia" ...
##  $ type           : chr  "earthquake" "earthquake" "earthquake" "earthquake" ...
##  $ horizontalError: num  10.71 6.94 5.02 8.51 4.35 ...
##  $ depthError     : num  7.94 7.33 5.32 8.75 7.54 ...
##  $ magError       : num  0.053 0.114 0.057 0.153 0.116 0.086 0.08 0.08 0.107 0.125 ...
##  $ magNst         : int  111 22 97 13 22 13 47 49 25 18 ...
##  $ status         : chr  "reviewed" "reviewed" "reviewed" "reviewed" ...
##  $ locationSource : chr  "us" "us" "us" "us" ...
##  $ magSource      : chr  "us" "us" "us" "us" ...
head(gempa_df[, c("time", "latitude", "longitude", "depth", "mag")])
##                       time latitude longitude   depth mag
## 1 2026-08-07T14:07:39.335Z  -5.5972  110.2073 534.988 4.8
## 2 2026-07-28T13:29:48.071Z  -8.9665  112.8235  62.529 4.4
## 3 2026-07-28T09:21:24.654Z  -7.4773  108.3692 155.549 4.9
## 4 2026-07-24T05:00:13.219Z  -7.7602  111.9400 163.615 4.2
## 5 2026-07-13T02:54:43.466Z  -7.2084  111.9426 219.978 4.3
## 6 2026-06-27T07:47:24.293Z  -8.7198  111.2648  71.683 5.2

Penjelasan parameter URL:

  • minlatitude, maxlatitude, minlongitude, dan maxlongitude digunakan untuk menentukan batas wilayah pengamatan.
  • minmagnitude=4.0 berarti hanya gempa dengan magnitudo 4.0 atau lebih yang diambil, supaya gempa yang sangat kecil tidak terlalu mendominasi data.
  • starttime dan endtime digunakan untuk menentukan periode waktu pengamatan.

Kolom yang paling penting untuk analisis point pattern adalah:

  • longitude = posisi bujur
  • latitude = posisi lintang

Kedua kolom tersebut digunakan sebagai koordinat lokasi kejadian gempa.

Selain itu, depth menunjukkan kedalaman gempa dan mag menunjukkan magnitudo gempa.

7.2 Dari Koordinat Geografis ke Koordinat Datar (Proyeksi)

Data longitude dan latitude masih berupa koordinat geografis dalam satuan derajat. Sementara itu, analisis pada spatstat lebih mudah dilakukan menggunakan koordinat datar dengan satuan jarak, misalnya meter.

Karena itu, koordinat gempa perlu diproyeksikan terlebih dahulu.

Pada wilayah Jawa dengan rentang sekitar 108°–114° BT, digunakan sistem proyeksi UTM Zona 49S (EPSG:32749).

# Ubah data.frame menjadi objek spasial (CRS geografis WGS84)
gempa_sf <- st_as_sf(gempa_df, coords = c("longitude", "latitude"), crs = 4326)

# Proyeksikan ke UTM Zona 49S (satuan meter)
gempa_utm <- st_transform(gempa_sf, crs = 32749)

# Ekstrak koordinat hasil proyeksi
koord <- st_coordinates(gempa_utm)
head(koord)
##             X       Y
## [1,] 412207.3 9381262
## [2,] 700479.0 9008354
## [3,] 209631.5 9172619
## [4,] 603651.4 9142097
## [5,] 604068.8 9203106
## [6,] 529127.4 9036116

Interpretasi: Nilai X dan Y merupakan posisi lokasi gempa yang sudah diubah dari koordinat lintang-bujur menjadi koordinat datar dalam satuan meter. Setiap baris mewakili satu lokasi gempa.

Koordinat ini selanjutnya digunakan untuk menghitung jarak antartitik dan membentuk objek point pattern.

7.3 Membentuk Objek ppp

Koordinat gempa yang sudah diproyeksikan kemudian diubah menjadi objek ppp agar dapat dianalisis menggunakan spatstat.

win_gempa berfungsi sebagai batas wilayah pengamatan, sedangkan koordinat X dan Y menjadi lokasi setiap kejadian gempa.

bbox_utm <- st_bbox(gempa_utm)

win_gempa <- owin(
  xrange = c(bbox_utm["xmin"], bbox_utm["xmax"]),
  yrange = c(bbox_utm["ymin"], bbox_utm["ymax"])
)

gempa_ppp <- ppp(x = koord[, "X"], y = koord[, "Y"], window = win_gempa)
unitname(gempa_ppp) <- c("meter", "meters")

gempa_ppp
## Planar point pattern: 306 points
## window: rectangle = [168508.2, 824884.5] x [9004636, 9386635] meters
plot(gempa_ppp, main = "Pola Titik Gempa Bumi M >= 4.0 di Sekitar Jawa (2018-sekarang)",
     pch = 16, cols = "firebrick", cex = 0.6)

Interpretasi: Setiap titik merah pada grafik menunjukkan lokasi kejadian gempa bumi dengan magnitudo ≥ 4,0 di sekitar Jawa. Setelah dibentuk menjadi objek ppp, data gempa sudah siap digunakan untuk analisis point pattern, seperti intensitas, CSR, fungsi-K, dan KDE.

7.4 Menerapkan Kembali Teknik yang Sudah Dipelajari

7.4.1 Intensitas & Metode Quadrat

Pada tahap ini, teknik yang sudah dipelajari sebelumnya diterapkan kembali pada data gempa bumi.

Fungsi intensity() digunakan untuk menghitung rata-rata kepadatan titik gempa dalam wilayah pengamatan. Selanjutnya, metode quadrat count membagi wilayah menjadi beberapa kotak dan menghitung jumlah gempa pada masing-masing kotak.

intensity(gempa_ppp)
## [1] 1.220411e-09
plot(gempa_ppp, pch = 16, cex = 0.6, cols = "firebrick", main = "Quadrat Count: Gempa Bumi Jawa")
plot(quadratcount(gempa_ppp, nx = 4, ny = 3), add = TRUE, col = "blue", cex = 1.2)

Interpretasi: Jumlah titik gempa pada setiap kotak terlihat cukup berbeda. Beberapa kotak memiliki jumlah gempa yang tinggi, sedangkan kotak lainnya jauh lebih sedikit. Hal ini menunjukkan bahwa persebaran gempa di wilayah Jawa tidak terlihat merata secara visual dan cenderung lebih terkonsentrasi pada wilayah tertentu.

7.4.2 Uji CSR

Uji CSR digunakan untuk mengetahui apakah pola titik gempa tersebar secara acak atau tidak.

Hipotesis:

  • \(H_0\): pola titik gempa mengikuti CSR atau tersebar secara acak.
  • \(H_1\): pola titik gempa tidak mengikuti CSR.
quadrat.test(gempa_ppp, nx = 4, ny = 3)
## 
##  Chi-squared test of CSR using quadrat counts
## 
## data:  gempa_ppp
## X2 = 150.78, df = 11, p-value < 2.2e-16
## alternative hypothesis: two.sided
## 
## Quadrats: 4 by 3 grid of tiles

Interpretasi: Hasil uji menunjukkan p-value < 2,2e-16, yang jauh lebih kecil dari 0,05. Maka \(H_0\) ditolak. Artinya, pola lokasi gempa bumi tidak tersebar secara acak dan terdapat indikasi pola spasial tertentu, seperti pengelompokan.

set.seed(123)
env_gempa <- envelope(gempa_ppp, Kest, nsim = 99, verbose = FALSE)
plot(env_gempa, main = "Uji CSR — Fungsi-K: Gempa Bumi Jawa")

Interpretasi: Garis hitam menunjukkan pola gempa yang diamati, sedangkan area abu-abu merupakan envelope dari pola acak CSR.

Pada grafik, kurva observasi berada jauh di atas envelope hampir di seluruh jarak \(r\). Artinya, pola gempa bumi tidak mengikuti CSR dan cenderung mengelompok (clustered).

Hasil ini sejalan dengan uji quadrat sebelumnya yang juga menunjukkan bahwa lokasi gempa tidak tersebar secara acak.

7.4.3 KDE & Hot Spot

KDE digunakan untuk melihat wilayah yang memiliki kepadatan gempa lebih tinggi dibandingkan wilayah lainnya. Bandwidth dipilih otomatis menggunakan bw.diggle() agar hasil peta kepadatan lebih sesuai dengan pola data.

bw_gempa <- bw.diggle(gempa_ppp)
dens_gempa <- density(gempa_ppp, sigma = bw_gempa)
plot(dens_gempa, main = "Peta Kepadatan (Hot Spot) Gempa Bumi di Sekitar Jawa")
contour(dens_gempa, add = TRUE, col = "white", nlevels = 5)
plot(gempa_ppp, add = TRUE, pch = 16, cex = 0.4, col = "black")

Interpretasi: Area dengan warna yang lebih terang menunjukkan kepadatan kejadian gempa yang lebih tinggi, sedangkan area berwarna gelap menunjukkan kepadatan yang lebih rendah. Area dengan kepadatan tinggi dapat dianggap sebagai hot spot gempa.

Pola ini mendukung hasil analisis sebelumnya bahwa kejadian gempa di sekitar Jawa cenderung mengelompok pada wilayah tertentu, bukan tersebar merata.

7.5 Catatan Penting Saat Menggunakan Data dari Internet

  1. Periksa struktur data terlebih dahulu menggunakan str() dan head(). Nama kolom atau format data dari API dapat berubah sewaktu-waktu.

  2. Perhatikan sistem koordinat. Data longitude dan latitude sebaiknya diproyeksikan ke koordinat datar sebelum digunakan untuk analisis jarak atau point pattern.

  3. Perhatikan batas wilayah pengamatan (window). Pada contoh ini digunakan batas berbentuk kotak (bounding box). Untuk analisis yang lebih baik, sebaiknya gunakan batas wilayah geografis yang sebenarnya.

  4. Data dari internet dapat terus diperbarui. Karena itu, jumlah data dan hasil analisis bisa berbeda ketika kode dijalankan pada waktu yang berbeda.

8 Rangkuman

Konsep Fungsi spatstat / R yang Digunakan
Membuat dan membaca pola titik ppp(), data bawaan seperti swedishpines, bei, chorley
Mengimpor data dari internet read.csv()
Mengubah koordinat geografis ke koordinat datar st_as_sf(), st_transform()
Menghitung intensitas dasar intensity()
Menghitung jumlah titik per kuadrat quadratcount()
Menguji intensitas berdasarkan kovariat rhohat()
Menghitung rasio intensitas kasus dan kontrol relrisk()
Uji CSR dengan metode kuadrat quadrat.test()
Uji CSR dengan Fungsi-K Kest()
Uji CSR dengan envelope envelope()
Menentukan bandwidth KDE bw.diggle(), bw.ppl()
Membuat peta kepadatan / KDE density()
Menampilkan garis kontur contour()
Melihat pola clustering dengan Fungsi-L Lest()

Secara umum, analisis point pattern dilakukan untuk melihat bagaimana titik tersebar dalam suatu wilayah, apakah cenderung acak, mengelompok, atau teratur, serta untuk melihat wilayah dengan kepadatan tinggi (hot spot).

9 Latihan Mandiri

Sebagai pendalaman materi, silakan pelajari kembali modul dan latihan tambahan pada:

https://r-spatial.org/book/11-PointPattern.html

Beberapa hal yang dapat dicoba secara mandiri:

  1. Ulangi uji quadrat CSR pada data bei. Apakah hasilnya berbeda dengan data swedishpines? Jelaskan alasannya.

  2. Ubah nilai nx dan ny pada quadrat.test() untuk data swedishpines. Perhatikan apakah hasil uji berubah ketika ukuran dan jumlah kuadrat diubah.

  3. Pada data chorley, coba gunakan beberapa nilai sigma yang berbeda pada relrisk(). Amati bagaimana perubahan bandwidth memengaruhi bentuk area hot spot yang dihasilkan.

  4. Ubah wilayah pengamatan data gempa ke wilayah lain di Indonesia, misalnya Sumatra, dengan mengubah batas latitude dan longitude pada URL. Setelah itu, sesuaikan sistem proyeksi UTM dan periksa apakah pola gempa cenderung acak atau mengelompok.

  5. Ubah nilai minmagnitude pada data gempa, misalnya menjadi 3.0 atau 5.0. Bandingkan jumlah titik yang diperoleh dan lihat apakah kesimpulan dari analisis CSR juga berubah.

Referensi: Baddeley, A., Rubak, E., Turner, R. (2024). Spatial Point Patterns: Methodology and Application with R; Moraga, P. (2024). Spatial Statistics for Data Science: Theory and Practice with R.