# 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.