# Library yang dibutuhkan
library(gstat)
## Warning: package 'gstat' was built under R version 4.5.3
library(sp)
# 1. Inisialisasi Lokasi dan Parameter
set.seed(123)
n_stasiun <- 5
n_hari <- 20
total_obs <- n_stasiun * n_hari
# Koordinat stasiun dan ketinggian (X)
coords <- data.frame(x = runif(n_stasiun, 0, 10), y = runif(n_stasiun, 0, 10))
elevasi <- rnorm(n_stasiun, 200, 50)
beta0 <- 15
beta1 <- 0.5
b0hat <- NULL; b1hat <- NULL
# 2. Loop Simulasi (100 kali)
for (i in 1:100) {
# Bangkitkan Error Spasial (sederhana)
dist_mat <- as.matrix(dist(coords))
sp_effect <- chol(exp(-dist_mat/5)) %*% matrix(rnorm(n_stasiun * n_hari), nrow = n_stasiun)
# Data Frame Spasial Temporal
df <- data.frame(
stasiun = rep(1:n_stasiun, each = n_hari),
hari = rep(1:n_hari, n_stasiun),
X = rep(elevasi, each = n_hari)
)
# Bangkitkan Y (Curah Hujan)
eps <- rnorm(total_obs, 0, 1)
df$Y <- beta0 + beta1 * df$X + as.vector(t(sp_effect)) + eps
# Estimasi OLS
model <- lm(Y ~ X, data = df)
b0hat <- c(b0hat, coef(model)[1])
b1hat <- c(b1hat, coef(model)[2])
}
# 3. Hasil Estimasi
hasil <- matrix(c(mean(b0hat), sd(b0hat), mean(b1hat), sd(b1hat)), nrow = 2)
rownames(hasil) <- c("mean", "sd")
colnames(hasil) <- c("b0 (Intersep)", "b1 (Elevasi)")
print(hasil)
## b0 (Intersep) b1 (Elevasi)
## mean 15.0165000 0.499828159
## sd 0.4999195 0.002622574
# 1. Load Library
library(ggplot2)
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.5.3
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(tidyr)
# 2. Menyiapkan Data (Menggunakan logika simulasi sebelumnya)
set.seed(123)
n_stasiun <- 5
n_hari <- 10 # Kita ambil 10 hari untuk visualisasi agar tidak terlalu penuh
# Koordinat stasiun (Spasial)
stasiun_coords <- data.frame(
stasiun = paste0("Stasiun ", 1:n_stasiun),
lon = c(106.8, 106.9, 107.0, 106.85, 107.1),
lat = c(-6.1, -6.2, -6.15, -6.3, -6.25)
)
# Simulasi data curah hujan (Temporal)
df_plot <- expand.grid(
stasiun = stasiun_coords$stasiun,
hari = 1:n_hari
) %>%
left_join(stasiun_coords, by = "stasiun") %>%
mutate(
# Logika: Curah hujan dipengaruhi lokasi + tren waktu + noise
curah_hujan = 20 + (lon * 0.5) + (lat * 0.2) + sin(hari) * 5 + rnorm(n(), 0, 2)
)
# 3. Grafik Spasial-Temporal (Facet Grid)
# Menampilkan peta sebaran hujan per hari
ggplot(df_plot, aes(x = lon, y = lat)) +
geom_point(aes(size = curah_hujan, color = curah_hujan), alpha = 0.7) +
scale_color_gradient(low = "lightblue", high = "darkblue") +
facet_wrap(~hari, ncol = 5) + # Membagi tampilan per hari
theme_minimal() +
labs(
title = "Simulasi Sebaran Curah Hujan Spasial-Temporal",
subtitle = "Visualisasi 5 Stasiun Pengamatan selama 10 Hari",
x = "Bujur (Longitude)",
y = "Lintang (Latitude)",
size = "Curah Hujan (mm)",
color = "Intensitas"
) +
theme(legend.position = "bottom")

# 4. Alternatif: Grafik Time-Series per Stasiun
# Untuk melihat fluktuasi temporal secara lebih jelas
ggplot(df_plot, aes(x = hari, y = curah_hujan, color = stasiun)) +
geom_line(size = 1) +
geom_point() +
theme_classic() +
labs(
title = "Tren Temporal Curah Hujan per Stasiun",
x = "Hari ke-",
y = "Curah Hujan (mm)"
)
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
