Install packages

library(ncdf4)
library(raster)
## Warning: package 'raster' was built under R version 4.2.3
## Loading required package: sp
## Warning: package 'sp' was built under R version 4.2.3
## The legacy packages maptools, rgdal, and rgeos, underpinning the sp package,
## which was just loaded, will retire in October 2023.
## Please refer to R-spatial evolution reports for details, especially
## https://r-spatial.org/r/2023/05/15/evolution4.html.
## It may be desirable to make the sf package available;
## package maintainers should consider adding sf to Suggests:.
## The sp package is now running under evolution status 2
##      (status 2 uses the sf package in place of rgdal)
library(terra)
## Warning: package 'terra' was built under R version 4.2.3
## terra 1.7.46
library(Metrics)
## Warning: package 'Metrics' was built under R version 4.2.3
library(Hmisc)
## 
## Attaching package: 'Hmisc'
## The following objects are masked from 'package:terra':
## 
##     describe, mask, units, zoom
## The following objects are masked from 'package:raster':
## 
##     mask, zoom
## The following objects are masked from 'package:base':
## 
##     format.pval, units
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.2.3
library(viridis)
## Warning: package 'viridis' was built under R version 4.2.3
## Loading required package: viridisLite
## Warning: package 'viridisLite' was built under R version 4.2.3
library(corrplot)
## Warning: package 'corrplot' was built under R version 4.2.3
## corrplot 0.92 loaded
library(gridExtra)
## Warning: package 'gridExtra' was built under R version 4.2.3
  1. Buka file NetCDF # ——————-
print("Membuka file NetCDF...")
## [1] "Membuka file NetCDF..."
obs_nc <- nc_open("FIX-observasi era5 suhu_1981-2010.nc")
model_nc <- nc_open("FIX-model cmip6 suhu_1981-2010.nc")
koreksi_nc <- nc_open("suhu_koreksi.nc")
  1. Ekstrak data dan koordinat # —————————–
print("Mengekstrak data suhu dan koordinat...")
## [1] "Mengekstrak data suhu dan koordinat..."
# Ekstrak variabel suhu
obs_temp <- ncvar_get(obs_nc, "mean_t2m")
model_temp <- ncvar_get(model_nc, "mean_t2m")
koreksi_temp <- ncvar_get(koreksi_nc, "tas")

# Ekstrak koordinat
lat_obs <- ncvar_get(obs_nc, "lat")
lon_obs <- ncvar_get(obs_nc, "lon")
lat_model <- ncvar_get(model_nc, "lat")
lon_model <- ncvar_get(model_nc, "lon")
lat_koreksi <- ncvar_get(koreksi_nc, "lat")
lon_koreksi <- ncvar_get(koreksi_nc, "lon")

# Ekstrak waktu untuk analisis temporal
time_obs <- ncvar_get(obs_nc, "time")
time_model <- ncvar_get(model_nc, "time")
time_koreksi <- ncvar_get(koreksi_nc, "time")

print("Dimensi data:")
## [1] "Dimensi data:"
print(paste("Observasi:", paste(dim(obs_temp), collapse = " x ")))
## [1] "Observasi: 185 x 81 x 10957"
print(paste("Model:", paste(dim(model_temp), collapse = " x ")))
## [1] "Model: 185 x 81 x 10957"
print(paste("Koreksi:", paste(dim(koreksi_temp), collapse = " x ")))
## [1] "Koreksi: 185 x 81 x 10957"
  1. Periksa urutan latitude dan standarisasi # ——————————————-
print("Urutan latitude:")
## [1] "Urutan latitude:"
print(paste("Observasi (pertama-terakhir):", lat_obs[1], "-", lat_obs[length(lat_obs)]))
## [1] "Observasi (pertama-terakhir): 8 - -12"
print(paste("Model (pertama-terakhir):", lat_model[1], "-", lat_model[length(lat_model)]))
## [1] "Model (pertama-terakhir): -12 - 8"
print(paste("Koreksi (pertama-terakhir):", lat_koreksi[1], "-", lat_koreksi[length(lat_koreksi)]))
## [1] "Koreksi (pertama-terakhir): -12 - 8"

Standardisasi urutan latitude observasi (balik jika descending), hal ini dilakukan karena urutan latitude pada file observasi beda dengan file model dan hasil koreksi.

if (lat_obs[1] > lat_obs[length(lat_obs)]) {
  print("Membalik urutan latitude observasi dari descending ke ascending...")
  lat_obs <- rev(lat_obs)
  # Balik urutan dimensi latitude dalam data (dimensi ke-2)
  obs_temp <- obs_temp[, rev(seq_len(dim(obs_temp)[2])), ]
}
## [1] "Membalik urutan latitude observasi dari descending ke ascending..."
# Verifikasi setelah pembalikan
print("Setelah standardisasi:")
## [1] "Setelah standardisasi:"
print(paste("Observasi (pertama-terakhir):", lat_obs[1], "-", lat_obs[length(lat_obs)]))
## [1] "Observasi (pertama-terakhir): -12 - 8"
  1. Subset data untuk wilayah Indonesia # ————————————–
print("Melakukan subset data untuk wilayah Indonesia...")
## [1] "Melakukan subset data untuk wilayah Indonesia..."
# Definisi batas wilayah Indonesia
# Latitude: 6° LU sampai 11° LS = 6°N sampai 11°S = 6 sampai -11
# Longitude: 95° BT sampai 141° BT = 95°E sampai 141°E = 95 sampai 141
lat_min <- -11  # 11° LS
lat_max <- 6    # 6° LU
lon_min <- 95   # 95° BT
lon_max <- 141  # 141° BT

print(paste("Batas wilayah Indonesia:"))
## [1] "Batas wilayah Indonesia:"
print(paste("Latitude:", lat_min, "sampai", lat_max))
## [1] "Latitude: -11 sampai 6"
print(paste("Longitude:", lon_min, "sampai", lon_max))
## [1] "Longitude: 95 sampai 141"

Fungsi untuk mencari indeks yang sesuai dengan batas koordinat

find_coord_indices <- function(coords, min_val, max_val) {
  indices <- which(coords >= min_val & coords <= max_val)
  return(indices)
}

# Cari indeks untuk wilayah Indonesia pada setiap dataset
lat_idx_obs <- find_coord_indices(lat_obs, lat_min, lat_max)
lon_idx_obs <- find_coord_indices(lon_obs, lon_min, lon_max)

lat_idx_model <- find_coord_indices(lat_model, lat_min, lat_max)
lon_idx_model <- find_coord_indices(lon_model, lon_min, lon_max)

lat_idx_koreksi <- find_coord_indices(lat_koreksi, lat_min, lat_max)
lon_idx_koreksi <- find_coord_indices(lon_koreksi, lon_min, lon_max)

print("Jumlah grid point yang tercakup:")
## [1] "Jumlah grid point yang tercakup:"
print(paste("Observasi - Lat:", length(lat_idx_obs), "Lon:", length(lon_idx_obs)))
## [1] "Observasi - Lat: 69 Lon: 185"
print(paste("Model - Lat:", length(lat_idx_model), "Lon:", length(lon_idx_model)))
## [1] "Model - Lat: 69 Lon: 185"
print(paste("Koreksi - Lat:", length(lat_idx_koreksi), "Lon:", length(lon_idx_koreksi)))
## [1] "Koreksi - Lat: 69 Lon: 185"
# Subset data untuk wilayah Indonesia
print("Melakukan subset data...")
## [1] "Melakukan subset data..."
obs_temp_indo <- obs_temp[lon_idx_obs, lat_idx_obs, ]
model_temp_indo <- model_temp[lon_idx_model, lat_idx_model, ]
koreksi_temp_indo <- koreksi_temp[lon_idx_koreksi, lat_idx_koreksi, ]

# Subset koordinat
lat_obs_indo <- lat_obs[lat_idx_obs]
lon_obs_indo <- lon_obs[lon_idx_obs]
lat_model_indo <- lat_model[lat_idx_model]
lon_model_indo <- lon_model[lon_idx_model]
lat_koreksi_indo <- lat_koreksi[lat_idx_koreksi]
lon_koreksi_indo <- lon_koreksi[lon_idx_koreksi]

print("Dimensi data setelah subset:")
## [1] "Dimensi data setelah subset:"
print(paste("Observasi:", paste(dim(obs_temp_indo), collapse = " x ")))
## [1] "Observasi: 185 x 69 x 10957"
print(paste("Model:", paste(dim(model_temp_indo), collapse = " x ")))
## [1] "Model: 185 x 69 x 10957"
print(paste("Koreksi:", paste(dim(koreksi_temp_indo), collapse = " x ")))
## [1] "Koreksi: 185 x 69 x 10957"
print("Rentang koordinat setelah subset:")
## [1] "Rentang koordinat setelah subset:"
print(paste("Obs - Lat:", min(lat_obs_indo), "sampai", max(lat_obs_indo)))
## [1] "Obs - Lat: -11 sampai 6"
print(paste("Obs - Lon:", min(lon_obs_indo), "sampai", max(lon_obs_indo)))
## [1] "Obs - Lon: 95 sampai 141"
print(paste("Model - Lat:", min(lat_model_indo), "sampai", max(lat_model_indo)))
## [1] "Model - Lat: -11 sampai 6"
print(paste("Model - Lon:", min(lon_model_indo), "sampai", max(lon_model_indo)))
## [1] "Model - Lon: 95 sampai 141"

Update variabel untuk analisis selanjutnya

obs_temp <- obs_temp_indo
model_temp <- model_temp_indo  
koreksi_temp <- koreksi_temp_indo
lat_obs <- lat_obs_indo
lon_obs <- lon_obs_indo
lat_model <- lat_model_indo
lon_model <- lon_model_indo
lat_koreksi <- lat_koreksi_indo
lon_koreksi <- lon_koreksi_indo
  1. Pastikan dimensi dan grid yang sama (untuk wilayah Indonesia) # ————————————————————— print(“Memeriksa konsistensi grid…”)

Periksa apakah grid sama

grid_match_lat <- length(lat_obs) == length(lat_model) && all(abs(lat_obs - lat_model) < 0.01) grid_match_lon <- length(lon_obs) == length(lon_model) && all(abs(lon_obs - lon_model) < 0.01)

if (!grid_match_lat || !grid_match_lon) { warning(“Grid latitude atau longitude tidak sama antara observasi dan model!”) print(“Mungkin perlu interpolasi lebih lanjut…”) }

  1. Pastikan dimensi dan grid yang sama (untuk wilayah Indonesia) # —————————————————————
print("Memeriksa konsistensi grid untuk wilayah Indonesia...")
## [1] "Memeriksa konsistensi grid untuk wilayah Indonesia..."
# Periksa apakah grid sama setelah subset
grid_match_lat <- length(lat_obs) == length(lat_model) && all(abs(lat_obs - lat_model) < 0.01)
grid_match_lon <- length(lon_obs) == length(lon_model) && all(abs(lon_obs - lon_model) < 0.01)

if (!grid_match_lat || !grid_match_lon) {
  warning("Grid latitude atau longitude tidak sama antara observasi dan model untuk wilayah Indonesia!")
  print("Menggunakan grid observasi sebagai referensi...")
  
  # Jika grid tidak sama, interpolasi model ke grid observasi
  print("Melakukan interpolasi model ke grid observasi...")
  # Di sini bisa ditambahkan kode interpolasi jika diperlukan
  # Untuk sementara, lanjutkan dengan asumsi grid cukup dekat
}

# Periksa dimensi waktu
if (dim(obs_temp)[3] != dim(model_temp)[3] || dim(obs_temp)[3] != dim(koreksi_temp)[3]) {
  warning("Dimensi waktu tidak sama antar dataset!")
  print(paste("Obs:", dim(obs_temp)[3], "Model:", dim(model_temp)[3], "Koreksi:", dim(koreksi_temp)[3]))
  
  # Ambil dimensi waktu terkecil untuk konsistensi
  min_time <- min(dim(obs_temp)[3], dim(model_temp)[3], dim(koreksi_temp)[3])
  obs_temp <- obs_temp[, , 1:min_time]
  model_temp <- model_temp[, , 1:min_time]
  koreksi_temp <- koreksi_temp[, , 1:min_time]
  print(paste("Menggunakan", min_time, "timestep untuk analisis"))
}
  1. Identifikasi area yang valid di Indonesia (bukan NaN) # ——————————————————–
print("Mengidentifikasi area valid...")
## [1] "Mengidentifikasi area valid..."
# Buat mask untuk area yang valid di semua dataset
valid_mask_obs <- !is.na(obs_temp[, , 1])
valid_mask_model <- !is.na(model_temp[, , 1])
valid_mask_koreksi <- !is.na(koreksi_temp[, , 1])

# Gabungkan mask - hanya area yang valid di semua dataset
combined_mask <- valid_mask_obs & valid_mask_model & valid_mask_koreksi

print(paste("Jumlah grid point valid:", sum(combined_mask)))
## [1] "Jumlah grid point valid: 12420"
print(paste("Persentase area valid:", round(sum(combined_mask)/length(combined_mask)*100, 2), "%"))
## [1] "Persentase area valid: 97.3 %"

Hitung metrik evaluasi # ————————

print("Menghitung metrik evaluasi...")
## [1] "Menghitung metrik evaluasi..."
# Fungsi untuk menghitung statistik hanya pada area valid
calculate_metrics_masked <- function(pred, obs, mask) {
  # Extract valid points
  pred_valid <- pred[mask]
  obs_valid <- obs[mask]
  
  # Remove any remaining NAs
  complete_cases <- complete.cases(pred_valid, obs_valid)
  pred_clean <- pred_valid[complete_cases]
  obs_clean <- obs_valid[complete_cases]
  
  if (length(pred_clean) == 0) {
    return(list(bias = NA, rmse = NA, mae = NA, correlation = NA, n = 0))
  }
  
  # Calculate metrics
  bias <- mean(pred_clean - obs_clean)
  rmse <- sqrt(mean((pred_clean - obs_clean)^2))
  mae <- mean(abs(pred_clean - obs_clean))
  correlation <- cor(pred_clean, obs_clean)
  
  return(list(
    bias = bias,
    rmse = rmse, 
    mae = mae,
    correlation = correlation,
    n = length(pred_clean)
  ))
}

# Hitung metrik untuk setiap timestep
n_time <- dim(obs_temp)[3]
metrics_model <- list()
metrics_koreksi <- list()

print("Menghitung metrik untuk setiap timestep...")
## [1] "Menghitung metrik untuk setiap timestep..."
for (t in 1:n_time) {
  if (t %% 500 == 0) {
    print(paste("Progress:", t, "/", n_time))
  }
  
  # Model vs Observasi
  metrics_model[[t]] <- calculate_metrics_masked(
    model_temp[, , t], obs_temp[, , t], combined_mask
  )
  
  # Koreksi vs Observasi  
  metrics_koreksi[[t]] <- calculate_metrics_masked(
    koreksi_temp[, , t], obs_temp[, , t], combined_mask
  )
}
## [1] "Progress: 500 / 10957"
## [1] "Progress: 1000 / 10957"
## [1] "Progress: 1500 / 10957"
## [1] "Progress: 2000 / 10957"
## [1] "Progress: 2500 / 10957"
## [1] "Progress: 3000 / 10957"
## [1] "Progress: 3500 / 10957"
## [1] "Progress: 4000 / 10957"
## [1] "Progress: 4500 / 10957"
## [1] "Progress: 5000 / 10957"
## [1] "Progress: 5500 / 10957"
## [1] "Progress: 6000 / 10957"
## [1] "Progress: 6500 / 10957"
## [1] "Progress: 7000 / 10957"
## [1] "Progress: 7500 / 10957"
## [1] "Progress: 8000 / 10957"
## [1] "Progress: 8500 / 10957"
## [1] "Progress: 9000 / 10957"
## [1] "Progress: 9500 / 10957"
## [1] "Progress: 10000 / 10957"
## [1] "Progress: 10500 / 10957"

Ekstrak metrik rata-rata

bias_model_ts <- sapply(metrics_model, function(x) x$bias)
rmse_model_ts <- sapply(metrics_model, function(x) x$rmse)
cor_model_ts <- sapply(metrics_model, function(x) x$correlation)

bias_koreksi_ts <- sapply(metrics_koreksi, function(x) x$bias)
rmse_koreksi_ts <- sapply(metrics_koreksi, function(x) x$rmse)
cor_koreksi_ts <- sapply(metrics_koreksi, function(x) x$correlation)
  1. Hitung metrik keseluruhan (menggunakan semua timestep) # ———————————————————
print("Menghitung metrik keseluruhan...")
## [1] "Menghitung metrik keseluruhan..."
# Fungsi untuk mengekstrak data dari array 3D menggunakan mask 2D
extract_valid_data <- function(data_3d, mask_2d) {
  # data_3d: array dengan dimensi [lon, lat, time]
  # mask_2d: logical array dengan dimensi [lon, lat]
  
  n_time <- dim(data_3d)[3]
  valid_data <- c()
  
  for (t in 1:n_time) {
    slice <- data_3d[, , t]
    valid_slice <- slice[mask_2d]
    valid_data <- c(valid_data, valid_slice)
  }
  
  return(valid_data)
}

# Ekstrak data yang valid untuk semua timestep
print("Mengekstrak data valid untuk semua timestep...")
## [1] "Mengekstrak data valid untuk semua timestep..."
obs_all <- extract_valid_data(obs_temp, combined_mask)
model_all <- extract_valid_data(model_temp, combined_mask)
koreksi_all <- extract_valid_data(koreksi_temp, combined_mask)

print(paste("Jumlah data points:", length(obs_all)))
## [1] "Jumlah data points: 136085940"
# Remove NAs (jika masih ada)
complete_idx <- complete.cases(obs_all, model_all, koreksi_all)
obs_clean <- obs_all[complete_idx]
model_clean <- model_all[complete_idx]
koreksi_clean <- koreksi_all[complete_idx]

print(paste("Jumlah data points setelah cleaning:", length(obs_clean)))
## [1] "Jumlah data points setelah cleaning: 136085940"

Metrik keseluruhan

overall_metrics <- data.frame(
  Dataset = c("Model", "Bias Corrected"),
  Bias = c(mean(model_clean - obs_clean), mean(koreksi_clean - obs_clean)),
  RMSE = c(sqrt(mean((model_clean - obs_clean)^2)), sqrt(mean((koreksi_clean - obs_clean)^2))),
  MAE = c(mean(abs(model_clean - obs_clean)), mean(abs(koreksi_clean - obs_clean))),
  Correlation = c(cor(model_clean, obs_clean), cor(koreksi_clean, obs_clean)),
  N_points = c(length(model_clean), length(koreksi_clean))
)

print("=== METRIK EVALUASI KESELURUHAN ===")
## [1] "=== METRIK EVALUASI KESELURUHAN ==="
print(overall_metrics)
##          Dataset         Bias      RMSE       MAE Correlation  N_points
## 1          Model 1.041232e+00 1.7400892 1.3508841   0.5836115 136085940
## 2 Bias Corrected 2.287483e-06 0.8758516 0.6871537   0.8483105 136085940
  1. Rata-rata spasial untuk analisis temporal # ——————————————–
print("Menghitung rata-rata spasial...")
## [1] "Menghitung rata-rata spasial..."
# Fungsi untuk rata-rata spasial dengan mask
spatial_mean_masked <- function(data, mask) {
  result <- numeric(dim(data)[3])
  for (t in 1:dim(data)[3]) {
    slice <- data[, , t]
    result[t] <- mean(slice[mask], na.rm = TRUE)
  }
  return(result)
}

# Hitung rata-rata spasial
obs_mean <- spatial_mean_masked(obs_temp, combined_mask)
model_mean <- spatial_mean_masked(model_temp, combined_mask)
koreksi_mean <- spatial_mean_masked(koreksi_temp, combined_mask)

print("Rata-rata spasial berhasil dihitung")
## [1] "Rata-rata spasial berhasil dihitung"
print(paste("Panjang time series:", length(obs_mean)))
## [1] "Panjang time series: 10957"
  1. Metrik time series # ——————–
print("=== METRIK TIME SERIES (rata-rata spasial) ===")
## [1] "=== METRIK TIME SERIES (rata-rata spasial) ==="
ts_metrics <- data.frame(
  Dataset = c("Model", "Bias Corrected"),
  Bias = c(mean(model_mean - obs_mean, na.rm = TRUE), mean(koreksi_mean - obs_mean, na.rm = TRUE)),
  RMSE = c(sqrt(mean((model_mean - obs_mean)^2, na.rm = TRUE)), sqrt(mean((koreksi_mean - obs_mean)^2, na.rm = TRUE))),
  Correlation = c(cor(model_mean, obs_mean, use = "complete.obs"), cor(koreksi_mean, obs_mean, use = "complete.obs"))
)
print(ts_metrics)
##          Dataset         Bias      RMSE Correlation
## 1          Model 1.041232e+00 1.1424152   0.5171547
## 2 Bias Corrected 2.287483e-06 0.4081406   0.5124111
  1. Visualisasi # —————
print("Membuat visualisasi...")
## [1] "Membuat visualisasi..."
# Plot time series rata-rata spasial Indonesia
if (length(time_obs) == length(obs_mean)) {
  # Konversi waktu jika dalam format hari sejak epoch
  if (mean(time_obs, na.rm = TRUE) > 10000) {
    # Asumsi: hari sejak 1900-01-01 atau 1950-01-01
    time_dates <- as.Date(time_obs, origin = "1900-01-01")  # Sesuaikan origin
  } else {
    time_dates <- 1:length(obs_mean)  # Gunakan indeks jika konversi gagal
  }
} else {
  time_dates <- 1:length(obs_mean)
}

# Data frame untuk plotting
ts_df <- data.frame(
  Time = rep(time_dates, 3),
  Temperature = c(obs_mean, model_mean, koreksi_mean),
  Dataset = rep(c("Observation", "Model", "Bias Corrected"), each = length(obs_mean))
)

ls()
##  [1] "bias_koreksi_ts"          "bias_model_ts"           
##  [3] "calculate_metrics_masked" "combined_mask"           
##  [5] "complete_idx"             "cor_koreksi_ts"          
##  [7] "cor_model_ts"             "extract_valid_data"      
##  [9] "find_coord_indices"       "grid_match_lat"          
## [11] "grid_match_lon"           "koreksi_all"             
## [13] "koreksi_clean"            "koreksi_mean"            
## [15] "koreksi_nc"               "koreksi_temp"            
## [17] "koreksi_temp_indo"        "lat_idx_koreksi"         
## [19] "lat_idx_model"            "lat_idx_obs"             
## [21] "lat_koreksi"              "lat_koreksi_indo"        
## [23] "lat_max"                  "lat_min"                 
## [25] "lat_model"                "lat_model_indo"          
## [27] "lat_obs"                  "lat_obs_indo"            
## [29] "lon_idx_koreksi"          "lon_idx_model"           
## [31] "lon_idx_obs"              "lon_koreksi"             
## [33] "lon_koreksi_indo"         "lon_max"                 
## [35] "lon_min"                  "lon_model"               
## [37] "lon_model_indo"           "lon_obs"                 
## [39] "lon_obs_indo"             "metrics_koreksi"         
## [41] "metrics_model"            "model_all"               
## [43] "model_clean"              "model_mean"              
## [45] "model_nc"                 "model_temp"              
## [47] "model_temp_indo"          "n_time"                  
## [49] "obs_all"                  "obs_clean"               
## [51] "obs_mean"                 "obs_nc"                  
## [53] "obs_temp"                 "obs_temp_indo"           
## [55] "overall_metrics"          "rmse_koreksi_ts"         
## [57] "rmse_model_ts"            "spatial_mean_masked"     
## [59] "t"                        "time_dates"              
## [61] "time_koreksi"             "time_model"              
## [63] "time_obs"                 "ts_df"                   
## [65] "ts_metrics"               "valid_mask_koreksi"      
## [67] "valid_mask_model"         "valid_mask_obs"
dim(valid_mask_obs)
## [1] 185  69
# Hitung jumlah grid valid pada hasil koreksi
valid_indonesia_grid <- sum(valid_mask_koreksi, na.rm = TRUE)

Plot time series untuk Indonesia

p1 <- ggplot(ts_df, aes(x = Time, y = Temperature, color = Dataset)) +
  geom_line(alpha = 0.7) +
  scale_color_manual(values = c("Observation" = "black", "Model" = "red", "Bias Corrected" = "blue")) +
  labs(title = "Spatially Averaged Temperature Time Series - Indonesia Region",
       subtitle = paste("Area: 6°N-11°S, 95°E-141°E |", valid_indonesia_grid, "valid grid points"),
       x = "Time", y = "Temperature (°C)") +
  theme_minimal() +
  theme(legend.position = "bottom")

# Plot distribusi bias untuk Indonesia
bias_df <- data.frame(
  Bias = c(model_mean - obs_mean, koreksi_mean - obs_mean),
  Dataset = rep(c("Model", "Bias Corrected"), each = length(obs_mean))
)

p2 <- ggplot(bias_df, aes(x = Bias, fill = Dataset)) +
  geom_histogram(alpha = 0.7, bins = 50, position = "identity") +
  scale_fill_manual(values = c("Model" = "red", "Bias Corrected" = "blue")) +
  labs(title = "Distribution of Bias - Indonesia Region", 
       x = "Bias (°C)", y = "Frequency") +
  theme_minimal() +
  theme(legend.position = "bottom")

# Scatter plot untuk Indonesia
scatter_df <- data.frame(
  Observed = rep(obs_mean, 2),
  Predicted = c(model_mean, koreksi_mean),
  Dataset = rep(c("Model", "Bias Corrected"), each = length(obs_mean))
)

p3 <- ggplot(scatter_df, aes(x = Observed, y = Predicted, color = Dataset)) +
  geom_point(alpha = 0.5) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "gray") +
  scale_color_manual(values = c("Model" = "red", "Bias Corrected" = "blue")) +
  labs(title = "Observed vs Predicted Temperature - Indonesia Region",
       x = "Observed Temperature (°C)", y = "Predicted Temperature (°C)") +
  theme_minimal() +
  theme(legend.position = "bottom")

# Plot spatial map untuk Indonesia (contoh satu timestep)
if (require(maps, quietly = TRUE)) {
  # Buat data untuk spatial plot
  sample_time <- 1  # Ambil timestep pertama sebagai contoh
  
  # Data untuk plotting
  spatial_df <- expand.grid(
    Longitude = lon_obs,
    Latitude = lat_obs
  )
  
  spatial_df$Observation <- as.vector(obs_temp[, , sample_time])
  spatial_df$Model <- as.vector(model_temp[, , sample_time])
  spatial_df$Bias_Corrected <- as.vector(koreksi_temp[, , sample_time])
  spatial_df$Bias_Original <- spatial_df$Model - spatial_df$Observation
  spatial_df$Bias_Corrected_diff <- spatial_df$Bias_Corrected - spatial_df$Observation
  
  # Plot spatial temperature
  p4 <- ggplot(spatial_df, aes(x = Longitude, y = Latitude)) +
    geom_tile(aes(fill = Observation)) +
    scale_fill_viridis_c(name = "Temp (°C)") +
    borders("world", xlim = c(lon_min-5, lon_max+5), ylim = c(lat_min-5, lat_max+5), colour = "white") +
    coord_cartesian(xlim = c(lon_min, lon_max), ylim = c(lat_min, lat_max)) +
    labs(title = "Temperature Distribution - Indonesia Region (Sample Time)",
         subtitle = "ERA5 Observation Data") +
    theme_minimal()
  
  # Plot bias spatial
  p5 <- ggplot(spatial_df, aes(x = Longitude, y = Latitude)) +
    geom_tile(aes(fill = Bias_Original)) +
    scale_fill_gradient2(low = "blue", mid = "white", high = "red", midpoint = 0, name = "Bias (°C)") +
    borders("world", xlim = c(lon_min-5, lon_max+5), ylim = c(lat_min-5, lat_max+5), colour = "white") +
    coord_cartesian(xlim = c(lon_min, lon_max), ylim = c(lat_min, lat_max)) +
    labs(title = "Model Bias - Indonesia Region (Sample Time)",
         subtitle = "Model - Observation") +
    theme_minimal()
  
  print(p4)
  print(p5)
}
## 
## Attaching package: 'maps'
## The following object is masked from 'package:viridis':
## 
##     unemp

Tampilkan plot utama Plot Time Series Indonesia

print(p1)

Plot distribusi bias untuk Indonesia

print(p2)

Scatter plot untuk Indonesia

print(p3)

  1. Ringkasan peningkatan untuk Indonesia # —————————————-
print("=== RINGKASAN PENINGKATAN UNTUK INDONESIA ===")
## [1] "=== RINGKASAN PENINGKATAN UNTUK INDONESIA ==="
bias_improvement <- abs(overall_metrics$Bias[1]) - abs(overall_metrics$Bias[2])
rmse_improvement <- overall_metrics$RMSE[1] - overall_metrics$RMSE[2]
cor_improvement <- overall_metrics$Correlation[2] - overall_metrics$Correlation[1]

print(paste("Wilayah Analisis: Indonesia (6°N-11°S, 95°E-141°E)"))
## [1] "Wilayah Analisis: Indonesia (6°N-11°S, 95°E-141°E)"
print(paste("Total Grid Points:", valid_indonesia_grid))
## [1] "Total Grid Points: 12420"
print(paste("Periode Analisis:", dim(obs_temp)[3], "timestep"))
## [1] "Periode Analisis: 10957 timestep"
print("")
## [1] ""
print("PENINGKATAN METRIK:")
## [1] "PENINGKATAN METRIK:"
print(paste("- Perbaikan Bias:", round(bias_improvement, 4), "°C"))
## [1] "- Perbaikan Bias: 1.0412 °C"
print(paste("- Perbaikan RMSE:", round(rmse_improvement, 4), "°C"))
## [1] "- Perbaikan RMSE: 0.8642 °C"
print(paste("- Perbaikan Korelasi:", round(cor_improvement, 4)))
## [1] "- Perbaikan Korelasi: 0.2647"
print(paste("- Persentase perbaikan RMSE:", round(rmse_improvement/overall_metrics$RMSE[1]*100, 2), "%"))
## [1] "- Persentase perbaikan RMSE: 49.67 %"

Tambahan statistik untuk Indonesia

print("")
## [1] ""
print("STATISTIK TAMBAHAN:")
## [1] "STATISTIK TAMBAHAN:"
temp_range_obs <- range(obs_mean, na.rm = TRUE)
temp_range_model <- range(model_mean, na.rm = TRUE)
temp_range_koreksi <- range(koreksi_mean, na.rm = TRUE)

print(paste("Rentang suhu observasi:", round(temp_range_obs[1], 2), "-", round(temp_range_obs[2], 2), "°C"))
## [1] "Rentang suhu observasi: 25.11 - 28.09 °C"
print(paste("Rentang suhu model:", round(temp_range_model[1], 2), "-", round(temp_range_model[2], 2), "°C"))
## [1] "Rentang suhu model: 25.8 - 29.06 °C"
print(paste("Rentang suhu koreksi:", round(temp_range_koreksi[1], 2), "-", round(temp_range_koreksi[2], 2), "°C"))
## [1] "Rentang suhu koreksi: 25.07 - 28.02 °C"

Hitung variabilitas

var_obs <- var(obs_mean, na.rm = TRUE)
var_model <- var(model_mean, na.rm = TRUE)
var_koreksi <- var(koreksi_mean, na.rm = TRUE)

print("")
## [1] ""
print("VARIABILITAS TEMPORAL:")
## [1] "VARIABILITAS TEMPORAL:"
print(paste("Variabilitas observasi:", round(var_obs, 4), "°C²"))
## [1] "Variabilitas observasi: 0.1617 °C²"
print(paste("Variabilitas model:", round(var_model, 4), "°C²"))
## [1] "Variabilitas model: 0.2789 °C²"
print(paste("Variabilitas koreksi:", round(var_koreksi, 4), "°C²"))
## [1] "Variabilitas koreksi: 0.1795 °C²"
  1. Tutup koneksi NetCDF # ————————
nc_close(obs_nc)
nc_close(model_nc)
nc_close(koreksi_nc)

print("Evaluasi koreksi bias selesai!")
## [1] "Evaluasi koreksi bias selesai!"