model sesuai tanggal tanam
Hani Zahra Anabela
2025-12-10
PACKAGES
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(lubridate)##
## Attaching package: 'lubridate'## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, unionlibrary(plotly)## Loading required package: ggplot2##
## Attaching package: 'plotly'## The following object is masked from 'package:ggplot2':
##
## last_plot## The following object is masked from 'package:stats':
##
## filter## The following object is masked from 'package:graphics':
##
## layoutlibrary(ggplot2)
library(tidyr)
library(patchwork)
library(writexl)
library(shiny)
library(shinydashboard)##
## Attaching package: 'shinydashboard'## The following object is masked from 'package:graphics':
##
## boxlibrary(readxl)
library(Matrix)##
## Attaching package: 'Matrix'## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpacklibrary(scales)
library(readr)##
## Attaching package: 'readr'## The following object is masked from 'package:scales':
##
## col_factorlibrary(reshape2) # Untuk melt data (jika diperlukan)##
## Attaching package: 'reshape2'## The following object is masked from 'package:tidyr':
##
## smithssetwd("~/KULIAH/5Semester5/MSP")Proses Data
############## Proses Data ##############
kab_bogor_df <- read_excel("bmkgbogor.xlsx", sheet = 1)
str(kab_bogor_df)## tibble [3,652 × 6] (S3: tbl_df/tbl/data.frame)
## $ Tanggal: POSIXct[1:3652], format: "2013-01-01" "2013-01-02" ...
## $ TA : num [1:3652] 21.3 21.6 21.4 20 20.8 20 20.4 20.3 19.3 21.8 ...
## $ RH : num [1:3652] 89 90 89 97 85 94 91 92 98 83 ...
## $ CH : num [1:3652] 10 47 0 0 26 42 8 9 71 66 ...
## $ S : num [1:3652] 1.4 5.2 2.2 2 0 3.4 0 0 0.5 0 ...
## $ FF : num [1:3652] 2.3 1.83 2.02 2.57 2.54 2.65 2.94 3.28 3.87 4.82 ...#kab_bogor_df <- kab_bogor_df[,c(3,12,15,4,8,10)]
colnames(kab_bogor_df) <- c('tgl','suhu','rh','ch', 'lama', 'angin')
#View(kab_bogor_df)
#kab_bogor_df <- arrange(kab_bogor_df, tgl)
str(kab_bogor_df)## tibble [3,652 × 6] (S3: tbl_df/tbl/data.frame)
## $ tgl : POSIXct[1:3652], format: "2013-01-01" "2013-01-02" ...
## $ suhu : num [1:3652] 21.3 21.6 21.4 20 20.8 20 20.4 20.3 19.3 21.8 ...
## $ rh : num [1:3652] 89 90 89 97 85 94 91 92 98 83 ...
## $ ch : num [1:3652] 10 47 0 0 26 42 8 9 71 66 ...
## $ lama : num [1:3652] 1.4 5.2 2.2 2 0 3.4 0 0 0.5 0 ...
## $ angin: num [1:3652] 2.3 1.83 2.02 2.57 2.54 2.65 2.94 3.28 3.87 4.82 ...#Penyesuaian
kab_bogor_df$tgl <- as.Date(kab_bogor_df$tgl, format="%Y-%m-%d")
kab_bogor_df$lama <- as.numeric(kab_bogor_df$lama)
kab_bogor_df$ch <- as.numeric(kab_bogor_df$ch)
kab_bogor_df$suhu <- as.numeric(kab_bogor_df$suhu)
kab_bogor_df$rh <- as.numeric(kab_bogor_df$rh)
kab_bogor_df$angin <- as.numeric(kab_bogor_df$angin)
#Pengecekan Data Kosong
print(paste('Jumlah Data NA pada CH: ', round(sum(is.na(kab_bogor_df$ch))/length(kab_bogor_df$ch)*100,2), '%', sep = ''))## [1] "Jumlah Data NA pada CH: 14.21%"print(paste('Jumlah Data NA pada Suhu: ', round(sum(is.na(kab_bogor_df$suhu))/length(kab_bogor_df$suhu)*100,2), '%', sep = ''))## [1] "Jumlah Data NA pada Suhu: 6.76%"print(paste('Jumlah Data NA pada RH: ', round(sum(is.na(kab_bogor_df$rh))/length(kab_bogor_df$rh)*100,2), '%', sep = ''))## [1] "Jumlah Data NA pada RH: 7.34%"print(paste('Jumlah Data NA pada Angin: ', round(sum(is.na(kab_bogor_df$angin))/length(kab_bogor_df$angin)*100,2), '%', sep = ''))## [1] "Jumlah Data NA pada Angin: 0%"print(paste('Jumlah Data NA pada Lama Penyinaran: ', round(sum(is.na(kab_bogor_df$lama))/length(kab_bogor_df$lama)*100,2), '%', sep = ''))## [1] "Jumlah Data NA pada Lama Penyinaran: 7.5%"summary(kab_bogor_df)## tgl suhu rh ch
## Min. :2013-01-01 Min. :18.3 Min. : 3.2 Min. : 0.00
## 1st Qu.:2015-07-02 1st Qu.:21.0 1st Qu.: 82.0 1st Qu.: 0.00
## Median :2017-12-31 Median :21.5 Median : 87.0 Median : 2.80
## Mean :2017-12-31 Mean :21.5 Mean : 86.0 Mean : 10.69
## 3rd Qu.:2020-06-30 3rd Qu.:22.0 3rd Qu.: 90.0 3rd Qu.: 14.40
## Max. :2022-12-30 Max. :57.5 Max. :100.0 Max. :192.80
## NA's :247 NA's :268 NA's :519
## lama angin
## Min. : 0.000 Min. :0.450
## 1st Qu.: 1.500 1st Qu.:1.050
## Median : 3.600 Median :1.290
## Mean : 3.744 Mean :1.412
## 3rd Qu.: 5.700 3rd Qu.:1.650
## Max. :10.000 Max. :4.820
## NA's :274############## Modifikasi Parameter Cuaca ##############
#Pengisian Data Kosong
kab_bogor_df$ch[is.na(kab_bogor_df$ch)] <- 0
kab_bogor_df$lama[is.na(kab_bogor_df$lama)] <- mean(kab_bogor_df$lama, na.rm = TRUE)
kab_bogor_df$lama <- kab_bogor_df$lama + 0.25
kab_bogor_df$suhu[is.na(kab_bogor_df$suhu)] <- mean(kab_bogor_df$suhu, na.rm = TRUE)
##kab_bogor_df$suhu <- kab_bogor_df$suhu + 2
kab_bogor_df$rh[is.na(kab_bogor_df$rh)] <- mean(kab_bogor_df$rh, na.rm = TRUE)
##kab_bogor_df$rh <- kab_bogor_df$rh - 20
##kab_bogor_df$rh <- pmin(kab_bogor_df$rh + 10, 100)
kab_bogor_df$angin[is.na(kab_bogor_df$angin)] <- mean(kab_bogor_df$angin, na.rm = TRUE)
summary(kab_bogor_df)## tgl suhu rh ch
## Min. :2013-01-01 Min. :18.3 Min. : 3.2 Min. : 0.000
## 1st Qu.:2015-07-02 1st Qu.:21.0 1st Qu.: 83.0 1st Qu.: 0.000
## Median :2017-12-31 Median :21.5 Median : 86.0 Median : 1.250
## Mean :2017-12-31 Mean :21.5 Mean : 86.0 Mean : 9.167
## 3rd Qu.:2020-06-30 3rd Qu.:22.0 3rd Qu.: 90.0 3rd Qu.: 11.400
## Max. :2022-12-30 Max. :57.5 Max. :100.0 Max. :192.800
## lama angin
## Min. : 0.250 Min. :0.450
## 1st Qu.: 1.950 1st Qu.:1.050
## Median : 3.994 Median :1.290
## Mean : 3.994 Mean :1.412
## 3rd Qu.: 5.775 3rd Qu.:1.650
## Max. :10.250 Max. :4.820sum(is.na(kab_bogor_df))## [1] 0############## CODING BARU ##############
##Masukkan parameter tanaman dan pendefinisian
k = as.numeric(0.4)
sla = as.numeric(0.07)
rue = as.numeric(1.3) #10.13031/trans.12144
div = as.numeric(1000000)
tb = as.numeric(15.5)
tgl_tanam = as.Date("2018-01-01", format="%Y-%m-%d")
print(tgl_tanam)## [1] "2018-01-01"# TU
suhutu1<- kab_bogor_df$suhu[1:5]
suhutu2<- kab_bogor_df$suhu[6:47]
suhutu3<- kab_bogor_df$suhu[48:82]
tu1 = abs(10*(tb-mean(suhutu1)))
tu2 = abs(42*(tb-mean(suhutu2)))
tu3 = abs(35*(tb-mean(suhutu3)))
tu = tu1+tu2+tu3
waktu_tanam = 87#Durasi tanam tanaman yang dipilih
#Julian date
kab_bogor_df$jd1970 <- julian(kab_bogor_df$tgl)RUNNING MODEL PER MUSIM TANAM
# ==============================================================================
# 1. PERSIAPAN DATA INPUT (Pastikan kab_bogor_df sudah bersih)
# ==============================================================================
# (Asumsi: kab_bogor_df sudah ada di environment kamu dari script sebelumnya)
# Pastikan kolom Julian Date ada
kab_bogor_df$jd1970 <- julian(kab_bogor_df$tgl)
# Tentukan Tanggal Tanam: Mulai 1 Jan 2018, interval 87 hari
tgl_mulai <- as.Date("2018-01-01")
tgl_akhir <- as.Date("2018-12-31")
list_tgl_tanam <- seq(from = tgl_mulai, to = tgl_akhir, by = 87)
# --- BAGIAN PENTING: FILTERING ---
# Kita hanya ambil tanggal tanam jika (Tanggal Tanam + 87 hari) MASIH di 2018.
# Jika panennya lompat ke 2016, tanggal tanam itu kita buang.
list_tgl_tanam <- list_tgl_tanam[ (list_tgl_tanam + 87) <= tgl_akhir]
print("Jadwal Tanam Final (Panen selesai di 2018):")## [1] "Jadwal Tanam Final (Panen selesai di 2018):"print(list_tgl_tanam)## [1] "2018-01-01" "2018-03-29" "2018-06-24" "2018-09-19"# ==============================================================================
# 2. UPDATE FUNGSI SIMULASI (Output Kolom Lengkap)
# ==============================================================================
run_simulation_full <- function(start_date, weather_db) {
# A. Filter Data Cuaca
data_slice <- weather_db %>%
filter(tgl >= start_date) %>%
slice(1:120)
if(nrow(data_slice) < 90) return(NULL)
# B. Parameter Tanaman
k = 0.4; sla = 0.07; rue = 1.3; div = 1000000; tb = 15.5
# C. Hitung Thermal Unit (TU)
suhu_slice <- data_slice$suhu
suhutu1 <- suhu_slice[1:5]
suhutu2 <- suhu_slice[6:47]
suhutu3 <- suhu_slice[48:82]
tu1 = abs(10*(tb-mean(suhutu1, na.rm=TRUE)))
tu2 = abs(42*(tb-mean(suhutu2, na.rm=TRUE)))
tu3 = abs(35*(tb-mean(suhutu3, na.rm=TRUE)))
tu = tu1+tu2+tu3
# D. Durasi Tanam
fase = rep(NA, nrow(data_slice))
s = rep(0, nrow(data_slice))
ds = rep(0, nrow(data_slice))
waktu_tanam_curr = 87
for (i in 2:nrow(data_slice)) {
if (s[i-1] >= 0.57) {
ds[i] = 0.43 * (suhu_slice[i] - tb) / tu3; fase[i] = 3
} else if (s[i-1] >= 0.06) {
ds[i] = 0.51 * (suhu_slice[i] - tb) / tu2; fase[i] = 2
} else {
ds[i] = 0.06 * (suhu_slice[i] - tb) / tu1; fase[i] = 1
}
s[i] = s[i-1] + ds[i]
if (s[i] >= 1) { waktu_tanam_curr = i; break }
}
# Potong Data
n_rows <- waktu_tanam_curr
data_crop <- data_slice[1:n_rows, ]
# E. Inisialisasi Variabel Lengkap
tggl <- data_crop$tgl; ch_vec <- data_crop$ch
lama_vec <- data_crop$lama; suhu_vec <- data_crop$suhu
rh_vec <- data_crop$rh; angin_vec <- data_crop$angin
lai <- numeric(n_rows); lai[1] <- 1.14
dW <- numeric(n_rows); Qint <- numeric(n_rows); Q10 <- numeric(n_rows)
wdaun <- numeric(n_rows); wbatang <- numeric(n_rows)
wakar <- numeric(n_rows); wbuah <- numeric(n_rows); wtot <- numeric(n_rows)
# Biomassa flow vars
dWdaun <- numeric(n_rows); dWbatang <- numeric(n_rows)
dWakar <- numeric(n_rows); dWbuah <- numeric(n_rows)
rdaun <- numeric(n_rows); rbatang <- numeric(n_rows)
rakar <- numeric(n_rows); rbuah <- numeric(n_rows)
# Neraca Air vars
KAT <- numeric(n_rows); KAT[1] <- 0.35
ic <- numeric(n_rows); ch_ef <- numeric(n_rows); ro <- numeric(n_rows)
ETo <- numeric(n_rows); ETc <- numeric(n_rows); wdf <- numeric(n_rows)
Rs <- numeric(n_rows); ta <- numeric(n_rows); tm <- numeric(n_rows)
em <- numeric(n_rows)
# Constants
KL <- 0.35; TLP <- 0.2832; miring <- 0.8; atm <- 100.299
Lat <- -6.479679; Alt <- 920; KATawal_val <- 0.35
# F. Loop Perhitungan Harian
for (i in 2:n_rows) {
# --- Radiasi & Penman ---
dr = 1 + 0.003 * cos((2 * pi * Lat / 365))
delta = 0.409 * sin((2 * pi * Lat / 365) - 1.39)
phi = Lat * (pi / 180)
omega_s = acos(-tan(phi) * tan((delta)))
Ra = ((cos(phi) * cos(delta) * sin((omega_s)) + (omega_s * sin(phi) * sin(delta))) * ((24 * 60 / pi) * 0.082 * dr))
N = (24 / pi) * omega_s
Rs[i] = (0.25 + 0.5 * (lama_vec[i] / N)) * Ra
Rs0 = (0.75 + 2 * 10 ^ -5 * Alt) * Ra
Rns = (1 - 0.23) * Rs[i]
es = 0.6108 * exp(17.27 * suhu_vec[i] / (suhu_vec[i] + 237.3))
ea_val = rh_vec[i] * es / 100
vpd = es - ea_val
delta_vap = 4098 * (0.6108 * exp((17.27 * suhu_vec[i]) / (suhu_vec[i] + 237.3)) / (suhu_vec[i] + 237.3) ^ 2)
gamma = 0.665 * 10 ^ (-3) * atm
Rn = Rns - (4.903e-9 * (suhu_vec[i] + 273.16)^4 * (0.34 - 0.14 * sqrt(ea_val)) * (1.35 * (Rs[i] / Rs0) - 0.35))
u2 = angin_vec[i]
ETo[i] = (0.408 * delta_vap * (Rn) + gamma * (900 / (suhu_vec[i] + 273)) * u2 * vpd) / (delta_vap + gamma * (1 + 0.34 * u2))
# --- Biomassa ---
Qint[i] = Rs[i] * (1 - exp(-k * lai[i-1]))
dW[i] = rue * Qint[i]
Q10[i] = 2 ^ ((suhu_vec[i] - 20) / 10)
km = 0.03
rdaun[i] = km * wdaun[i-1] * Q10[i]; rbatang[i] = km * wbatang[i-1] * Q10[i]
rakar[i] = km * wakar[i-1] * Q10[i]; rbuah[i] = km * wbuah[i-1] * Q10[i]
b_total = 905.78 + 422.35 + 148.52
if (s[i] < 0.33) {
dWdaun[i] = dW[i] * (905.78/b_total) - rdaun[i]
dWbatang[i] = dW[i] * (422.35/b_total) - rbatang[i]
dWakar[i] = dW[i] * (148.52/b_total) - rakar[i]
wdaun[i] = wdaun[i-1] + dWdaun[i]; wbatang[i] = wbatang[i-1] + dWbatang[i]
wakar[i] = wakar[i-1] + dWakar[i]; wbuah[i] = 0
} else {
dWbuah[i] = dW[i] - rbuah[i]
wbuah[i] = wbuah[i-1] + dWbuah[i]
wdaun[i] = wdaun[i-1] - rdaun[i]; wbatang[i] = wbatang[i-1] - rbatang[i]; wakar[i] = wakar[i-1] - rakar[i]
}
lai[i] = sla * wdaun[i]
wtot[i] = wdaun[i] + wakar[i] + wbatang[i] + wbuah[i]
# --- Neraca Air ---
kc_val <- if(i > 60) 1.15 else (0.5 + (i/60)*(1.15-0.5))
ETc[i] = ETo[i] * kc_val
ic[i] = if(lai[i] <= 3) 0.35 * lai[i] else 1.27
if(ch_vec[i] < ic[i]) ic[i] = ch_vec[i]
ch_ef[i] = ch_vec[i] - ic[i]
ro[i] = (sin(pi * miring / 180) * ch_ef[i] * (KAT[i-1] / KL)) / (lai[i] + 1)
em[i] = ETo[i] * exp(-k * lai[i]); if(em[i] < 0) em[i] = 0
tm[i] = ETo[i] - em[i]
if (KAT[i-1] < KL) {
ta[i] = if(KAT[i-1] > TLP) tm[i] * (KAT[i-1] - TLP) / (KL - TLP) else 0
} else { ta[i] = tm[i] }
delta_kat = (ch_ef[i] - ro[i] - ta[i]) / 1000
prev_kat = if(i==2) KATawal_val else KAT[i-1]
KAT[i] = prev_kat + delta_kat
if (KAT[i] > KL) KAT[i] = KL
if (KAT[i] < 0) KAT[i] = 0
wdf[i] = ta[i] / tm[i]
if (is.na(wdf[i]) || wdf[i] < 0.35) wdf[i] = 0.35
}
# G. PACKING HASIL (STRUKTUR DIPERBAIKI DISINI)
# Menyesuaikan nama kolom agar mirip hh_2013_2022
res_df <- data.frame(
Tanggal = tggl,
HST = 1:n_rows,
# --- DATA CUACA & INPUT ---
Suhu = suhu_vec,
CH = ch_vec,
Lama_Penyinaran = lama_vec, # Menggunakan underscore biar aman
Radiasi = Rs,
RH = rh_vec,
Angin = angin_vec,
# --- DATA BIOMASSA ---
LAI = lai,
Daun = wdaun,
Batang = wbatang,
Akar = wakar,
Buah = wbuah,
Total = wtot,
Fase = as.factor(fase[1:n_rows]),
# --- DATA NERACA AIR ---
Intersepsi = ic,
Run_off = ro,
ETo = ETo,
ETc = ETc,
em = em, # Evaporasi
tm = tm, # Transpirasi
wdf = wdf,
KAT = KAT,
Curah_hujan_efektif = ch_ef,
# --- PRODUKTIVITAS & IDENTITAS ---
Prod1HAinTon = (wbuah * 50000) / div,
Tanggal_Tanam_Siklus = rep(start_date, n_rows),
Siklus_Ke = NA
)
return(res_df)
}
# ==============================================================================
# 3. GENERATE ULANG DATA FRAME df_simulasi_2018
# ==============================================================================
all_results <- list()
for (k in seq_along(list_tgl_tanam)) {
tgl <- list_tgl_tanam[k]
hasil <- run_simulation_full(tgl, kab_bogor_df)
if (!is.null(hasil)) {
hasil$Siklus_Ke <- k
all_results[[k]] <- hasil
}
}
df_simulasi_2018 <- bind_rows(all_results)
# ==============================================================================
# 4. CEK HASIL & VISUALISASI
# ==============================================================================
print("Struktur Data Frame Baru:")## [1] "Struktur Data Frame Baru:"str(df_simulasi_2018)## 'data.frame': 354 obs. of 27 variables:
## $ Tanggal : Date, format: "2018-01-01" "2018-01-02" ...
## $ HST : int 1 2 3 4 5 6 7 8 9 10 ...
## $ Suhu : num 21.5 21.8 21 22.2 22 21.7 22.1 22 21 21.7 ...
## $ CH : num 0 0.6 0.5 20 0.4 21.2 0 5.3 22.2 7.2 ...
## $ Lama_Penyinaran : num 7.45 7.05 0.25 2.25 3.15 0.25 4.95 0.25 5.65 2.25 ...
## $ Radiasi : num 0 19.83 9.64 12.64 13.99 ...
## $ RH : num 80 81 88 87 80 82 79 83 88 84 ...
## $ Angin : num 1.5 1.79 1.25 1.95 2.45 2.35 2.09 1.84 1.27 1.95 ...
## $ LAI : num 1.14 0.405 0.473 0.578 0.719 ...
## $ Daun : num 0 5.79 6.76 8.26 10.27 ...
## $ Batang : num 0 2.7 3.15 3.85 4.79 ...
## $ Akar : num 0 0.949 1.108 1.354 1.685 ...
## $ Buah : num 0 0 0 0 0 0 0 0 0 0 ...
## $ Total : num 0 9.44 11.01 13.46 16.75 ...
## $ Fase : Factor w/ 3 levels "1","2","3": NA 1 1 1 1 1 1 1 1 1 ...
## $ Intersepsi : num 0 0.142 0.165 0.202 0.252 ...
## $ Run_off : num 0 0.00455 0.00317 0.17512 0.0012 ...
## $ ETo : num 0 3.67 2.04 2.58 3.06 ...
## $ ETc : num 0 1.92 1.09 1.4 1.7 ...
## $ em : num 0 3.12 1.69 2.05 2.3 ...
## $ tm : num 0 0.55 0.352 0.533 0.765 ...
## $ wdf : num 0 1 0.999 0.998 1 ...
## $ KAT : num 0.35 0.35 0.35 0.35 0.349 ...
## $ Curah_hujan_efektif : num 0 0.458 0.335 19.798 0.148 ...
## $ Prod1HAinTon : num 0 0 0 0 0 0 0 0 0 0 ...
## $ Tanggal_Tanam_Siklus: Date, format: "2018-01-01" "2018-01-01" ...
## $ Siklus_Ke : int 1 1 1 1 1 1 1 1 1 1 ...VIZ PPRODUKTIVITAS
# Visualisasi Perbandingan Produktivitas Antar Siklus Tanam
plot_multi_cycle <- ggplot(df_simulasi_2018, aes(x = HST, y = Prod1HAinTon, color = factor(Tanggal_Tanam_Siklus))) +
geom_line(size = 1) +
labs(
title = "Simulasi Pertumbuhan Timun: 4 Siklus Tanam di 2018",
subtitle = "Interval Tanam setiap ~87 Hari",
x = "Hari Setelah Tanam (HST)",
y = "Produktivitas (Ton/Ha)",
color = "Tanggal Tanam"
) +
theme_classic()## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.print(plot_multi_cycle)
ggplotly(plot_multi_cycle)#ggplot ver
ggplot(df_simulasi_2018, aes(x = HST, y = Prod1HAinTon, color = factor(Tanggal_Tanam_Siklus))) +
geom_line(size = 1) +
labs(
title = "Simulasi Pertumbuhan Timun: 4 Siklus Tanam di 2018",
x = "Hari Setelah Tanam (HST)",
y = "Produktivitas (Ton/Ha)",
color = "Tanggal Tanam"
) +
scale_y_continuous(limits = c(0, NA)) +
scale_x_continuous(limits = c(0, 90)) +
theme_classic()+
theme_classic(base_size = 14)+
theme(text = element_text(family = "serif"))+
theme(plot.title = element_text(face = "bold"))
MUSIM KERING VS. BASAH
# ==============================================================================
# 1. MENENTUKAN STATUS MUSIM PER SIKLUS (DATA DRIVEN)
# ==============================================================================
# Hitung ringkasan per siklus tanam
summary_siklus <- df_simulasi_2018 %>%
group_by(Siklus_Ke, Tanggal_Tanam_Siklus) %>%
summarise(
Total_CH_Siklus = sum(CH, na.rm = TRUE), # Total air yang diterima tanaman
Rata_Suhu_Siklus = mean(Suhu, na.rm = TRUE), # Rata-rata suhu
Panen_Akhir = max(Prod1HAinTon, na.rm = TRUE), # Hasil panen final (akumulasi maks)
.groups = "drop"
)
# Tentukan Threshold (Ambang Batas) berdasarkan rata-rata tahun itu
threshold_lokal <- mean(summary_siklus$Total_CH_Siklus)
print(paste("Rata-rata Curah Hujan per Siklus Tanam (87 hari):", round(threshold_lokal, 2), "mm"))## [1] "Rata-rata Curah Hujan per Siklus Tanam (87 hari): 731.65 mm"# Beri Label Musim
summary_siklus <- summary_siklus %>%
mutate(Status_Musim = ifelse(Total_CH_Siklus >= threshold_lokal,
"Musim Basah",
"Musim Kering"))
print("Klasifikasi Musim per Siklus:")## [1] "Klasifikasi Musim per Siklus:"print(summary_siklus %>% select(Siklus_Ke, Tanggal_Tanam_Siklus, Total_CH_Siklus, Status_Musim))## # A tibble: 4 × 4
## Siklus_Ke Tanggal_Tanam_Siklus Total_CH_Siklus Status_Musim
## <int> <date> <dbl> <chr>
## 1 1 2018-01-01 1435. Musim Basah
## 2 2 2018-03-29 526. Musim Kering
## 3 3 2018-06-24 233. Musim Kering
## 4 4 2018-09-19 732. Musim Kering# ==============================================================================
# 2. VISUALISASI PERBANDINGAN PRODUKTIVITAS
# ==============================================================================
# Gabungkan status musim kembali ke data utama (Opsional, jika ingin mewarnai grafik harian)
# Tapi untuk membandingkan produksi, kita cukup pakai data summary_siklus
# A. Bar Chart: Curah Hujan per Siklus (Untuk Validasi Musim)
p_ch <- ggplot(summary_siklus, aes(x = factor(Tanggal_Tanam_Siklus), y = Total_CH_Siklus, fill = Status_Musim)) +
geom_col() +
geom_hline(yintercept = threshold_lokal, linetype = "dashed", color = "red") +
annotate("text", x = 1, y = threshold_lokal + 50, label = "Threshold Rata-rata", color = "red", hjust=0) +
scale_fill_manual(values = c("Musim Basah" = "#0072B2", "Musim Kering" = "#E69F00")) +
labs(title = "Validasi: Total Hujan per Siklus Tanam",
y = "Total Curah Hujan (mm/87 hari)", x = "Tanggal Mulai Tanam") +
theme_classic()
# B. Bar Chart: Perbandingan Produktivitas (Basah vs Kering)
p_prod <- ggplot(summary_siklus, aes(x = Status_Musim, y = Panen_Akhir, fill = Status_Musim)) +
stat_summary(fun = mean, geom = "bar", width = 0.6, color = "black") +
geom_jitter(width = 0.1, size = 3, color = "white", stroke = 1) + # Menampilkan titik data asli
geom_text(stat = 'summary', fun = mean, aes(label = round(..y.., 2)), vjust = -0.5, fontface="bold") +
scale_fill_manual(values = c("Musim Basah" = "#0072B2", "Musim Kering" = "#E69F00")) +
labs(
title = "Perbandingan Produktivitas Timun",
subtitle = "Musim Basah vs Musim Kering (Berdasarkan CH Aktual)",
y = "Rata-rata Produktivitas (Ton/Ha)",
x = NULL
) +
theme_classic(base_size = 14) +
theme(legend.position = "none")
# Tampilkan
print(p_ch)
print(p_prod)## Warning: The dot-dot notation (`..y..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(y)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
MUSIM KERING VS. BASAH VER 2
# --- 1. Persiapan Data Agregat (Jaga-jaga jika belum run yang tadi) ---
# Kita butuh rata-rata produktivitas per Kategori Musim untuk Bar Chart
stat_musim <- summary_siklus %>%
group_by(Status_Musim) %>%
summarise(
Mean_Prod = mean(Panen_Akhir, na.rm = TRUE),
SD_Prod = sd(Panen_Akhir, na.rm = TRUE), # Opsional: untuk error bar
Count = n()
)
# --- 2. Membuat Plot A: Bar Chart (Rata-rata) ---
p_bar <- ggplot(stat_musim, aes(x = Status_Musim, y = Mean_Prod, fill = Status_Musim)) +
geom_col(width = 0.6, color = "black", alpha = 0.9) +
# Menambah Label Angka di atas batang
geom_text(aes(label = round(Mean_Prod, 2)),
vjust = -0.5, size = 5, family = "serif", fontface = "bold") +
# Skala & Warna Konsisten
scale_fill_manual(values = c("Musim Basah" = "#0072B2", "Musim Kering" = "#E69F00")) +
scale_y_continuous(expand = expansion(mult = c(0, 0.2))) + # Memberi ruang atas
# Labeling
labs(
title = "Rata-rata Produktivitas Timun",
subtitle = "Komparasi Rata-rata per Musim",
x = NULL,
y = "Produktivitas (Ton/Ha)"
) +
theme_classic(base_size = 14) +
theme(
text = element_text(family = "serif"),
legend.position = "none",
plot.title = element_text(face = "bold", hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5)
)
# --- 3. Membuat Plot B: Boxplot (Distribusi) ---
p_box <- ggplot(summary_siklus, aes(x = Status_Musim, y = Panen_Akhir, fill = Status_Musim)) +
# Jitter: Titik data asli (penting jika datanya sedikit)
geom_jitter(width = 0.1, size = 3, color = "grey40", alpha = 0.6) +
# Boxplot transparan di atasnya
geom_boxplot(width = 0.4, alpha = 0.7, outlier.shape = NA) +
# Tanda Mean (Wajik Putih)
stat_summary(fun = mean, geom = "point", shape = 23, size = 4,
fill = "white", color = "black", stroke = 1.5) +
# Skala & Warna
scale_fill_manual(values = c("Musim Basah" = "#0072B2", "Musim Kering" = "#E69F00")) +
# Labeling
labs(
title = "Distribusi Hasil Panen",
subtitle = "Sebaran Data per Siklus Tanam",
x = NULL,
y = NULL # Y axis dihapus biar tidak redundan
) +
scale_y_continuous()+
theme_classic(base_size = 14) +
theme(
text = element_text(family = "serif"),
legend.position = "none",
axis.text.y = element_blank(), # Hilangkan angka sumbu Y kanan
axis.ticks.y = element_blank(),
axis.line.y = element_blank(),
plot.title = element_text(face = "bold", hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5)
)
# --- 4. Menggabungkan dengan Patchwork ---
final_plot <- p_bar + p_box +
plot_layout(widths = c(1, 1)) + # Proporsi lebar 50:50
# Judul Utama Gabungan
plot_annotation(
title = "Analisis Dampak Iklim Terhadap Produktivitas Timun",
subtitle = "Simulasi Model Pertumbuhan Tanaman (Data Driven Season)",
theme = theme(
plot.title = element_text(size = 18, face = "bold", family = "serif", hjust = 0.5),
plot.subtitle = element_text(size = 13, family = "serif", hjust = 0.5, margin = margin(b=20))
)
)
# Tampilkan
print(final_plot)
## VIZ LAI SEASONAL
####VIZ####
str(df_simulasi_2018)## 'data.frame': 354 obs. of 27 variables:
## $ Tanggal : Date, format: "2018-01-01" "2018-01-02" ...
## $ HST : int 1 2 3 4 5 6 7 8 9 10 ...
## $ Suhu : num 21.5 21.8 21 22.2 22 21.7 22.1 22 21 21.7 ...
## $ CH : num 0 0.6 0.5 20 0.4 21.2 0 5.3 22.2 7.2 ...
## $ Lama_Penyinaran : num 7.45 7.05 0.25 2.25 3.15 0.25 4.95 0.25 5.65 2.25 ...
## $ Radiasi : num 0 19.83 9.64 12.64 13.99 ...
## $ RH : num 80 81 88 87 80 82 79 83 88 84 ...
## $ Angin : num 1.5 1.79 1.25 1.95 2.45 2.35 2.09 1.84 1.27 1.95 ...
## $ LAI : num 1.14 0.405 0.473 0.578 0.719 ...
## $ Daun : num 0 5.79 6.76 8.26 10.27 ...
## $ Batang : num 0 2.7 3.15 3.85 4.79 ...
## $ Akar : num 0 0.949 1.108 1.354 1.685 ...
## $ Buah : num 0 0 0 0 0 0 0 0 0 0 ...
## $ Total : num 0 9.44 11.01 13.46 16.75 ...
## $ Fase : Factor w/ 3 levels "1","2","3": NA 1 1 1 1 1 1 1 1 1 ...
## $ Intersepsi : num 0 0.142 0.165 0.202 0.252 ...
## $ Run_off : num 0 0.00455 0.00317 0.17512 0.0012 ...
## $ ETo : num 0 3.67 2.04 2.58 3.06 ...
## $ ETc : num 0 1.92 1.09 1.4 1.7 ...
## $ em : num 0 3.12 1.69 2.05 2.3 ...
## $ tm : num 0 0.55 0.352 0.533 0.765 ...
## $ wdf : num 0 1 0.999 0.998 1 ...
## $ KAT : num 0.35 0.35 0.35 0.35 0.349 ...
## $ Curah_hujan_efektif : num 0 0.458 0.335 19.798 0.148 ...
## $ Prod1HAinTon : num 0 0 0 0 0 0 0 0 0 0 ...
## $ Tanggal_Tanam_Siklus: Date, format: "2018-01-01" "2018-01-01" ...
## $ Siklus_Ke : int 1 1 1 1 1 1 1 1 1 1 ...# --- Definisi Tema Konsisten (Serif + Clean) ---
my_theme <- theme_classic(base_size = 14) +
theme(
text = element_text(family = "serif"),
plot.title = element_text(face = "bold", hjust = 0.5, size = 16),
plot.subtitle = element_text(hjust = 0.5, color = "gray40"),
legend.position = "bottom",
strip.background = element_rect(fill = "#f0f0f0", color = NA), # Latar judul facet
strip.text = element_text(face = "bold")
)
# --- Plot 1: Dinamika LAI per Siklus Tanam ---
p_lai <- ggplot(df_simulasi_2018, aes(x = HST, y = LAI, color = factor(Tanggal_Tanam_Siklus))) +
geom_line(linewidth = 1.2, alpha = 0.8) +
scale_color_brewer(palette = "Dark2") + # Warna pembeda siklus
labs(
title = "Dinamika Leaf Area Index (LAI)",
subtitle = "Perbandingan Antar Siklus Tanam (2018)",
x = "Hari Setelah Tanam (HST)",
y = "LAI",
color = "Tanggal Tanam"
) +
my_theme
print(p_lai)
## VIZ BIOMASS
# Pivot data agar formatnya long (untuk stacked chart)
biomassa_long <- df_simulasi_2018 %>%
select(HST, Tanggal_Tanam_Siklus, Daun, Batang, Akar, Buah) %>%
pivot_longer(cols = c(Daun, Batang, Akar, Buah),
names_to = "Organ", values_to = "Berat_Kering")
p_biom <- ggplot(biomassa_long, aes(x = HST, y = Berat_Kering, fill = Organ)) +
geom_area(alpha = 0.85, color = "white", size = 0.1) +
facet_wrap(~Tanggal_Tanam_Siklus, ncol = 2) + # Pisah per siklus
labs(
title = "Akumulasi Biomassa Timun",
subtitle = "Sebaran Organ per Siklus Tanam",
x = "HST", y = "Berat Kering (g/tanaman)"
) +
theme_classic(base_size = 14) +
theme(
text = element_text(family = "serif"),
plot.title = element_text(face = "bold", hjust = 0.5),
legend.position = "bottom"
)
print(p_biom)
KAT VS. CH SEASONAL
# --- 1. Ambil Tanggal Batas Musim Tanam ---
# Kita ambil tanggal unik dari kolom 'Tanggal_Tanam_Siklus' sebagai batas
batas_musim <- unique(as.Date(df_simulasi_2018$Tanggal_Tanam_Siklus))
# --- 2. Plot KAT (Atas) ---
p_kat <- ggplot(df_simulasi_2018, aes(x = Tanggal, y = KAT)) +
# Garis KAT
geom_line(color = "#0072B2", linewidth = 1) +
# Garis Vertikal Pemisah Musim
geom_vline(xintercept = batas_musim, linetype = "dashed", color = "gray9", alpha = 0.8) +
# Label Siklus di Atas
annotate("text", x = batas_musim + 5, y = max(df_simulasi_2018$KAT) * 1.02,
label = paste("Musim", seq_along(batas_musim)),
hjust = 0, vjust = 0, size = 3.5, family = "serif", color = "black", fontface="bold") +
# Format Axis X (Bulan)
scale_x_date(date_breaks = "1 month", date_labels = "%b") +
labs(title = "Dinamika Kadar Air Tanah (KAT)",
y = "KAT (Fraksi)",
x = NULL) + # X label dihapus agar tidak duplikat dengan grafik bawah
theme_classic(base_size = 14) +
theme(
text = element_text(family = "serif"),
plot.title = element_text(face = "bold"),
axis.text.x = element_blank(), # Hilangkan teks bulan di grafik atas agar rapi
axis.ticks.x = element_blank()
)
# --- 3. Plot Curah Hujan (Bawah) ---
p_ch <- ggplot(df_simulasi_2018, aes(x = Tanggal, y = CH)) +
# Bar Curah Hujan
geom_col(fill = "#E69F00", alpha = 0.8, width = 1) +
# Garis Vertikal Pemisah Musim (Agar sejajar dengan atas)
geom_vline(xintercept = batas_musim, linetype = "dashed", color = "gray9", alpha = 0.8) +
# Format Axis X (Bulan)
scale_x_date(date_breaks = "1 month", date_labels = "%b") +
labs(title = "Curah Hujan Harian",
y = "Curah Hujan (mm)",
x = "Bulan (Tahun 2018)") +
theme_classic(base_size = 14) +
theme(
text = element_text(family = "serif"),
plot.title = element_text(face = "bold")
)
# --- 4. Gabungkan dengan Patchwork ---
combined_plot <- p_kat / p_ch +
plot_layout(heights = c(1, 1)) + # Proporsi tinggi sama
plot_annotation(
title = "Hubungan Curah Hujan dan Kadar Air Tanah",
subtitle = "Simulasi Sepanjang Tahun 2018 dengan Pemisah Musim Tanam",
theme = theme(
plot.title = element_text(size = 18, face = "bold", family = "serif", hjust = 0.5),
plot.subtitle = element_text(size = 12, family = "serif", hjust = 0.5)
)
)
# Tampilkan
print(combined_plot)
#THEME
# ==============================================================================
# 1. DEFINISI TEMA KUSTOM (Agar tampilan konsisten dan rapi)
# ==============================================================================
# Kita membuat tema dasar dari theme_classic lalu memodifikasinya
theme_publication <- theme_classic(base_size = 14, base_family = "serif") +
theme(
# Judul Plot
plot.title = element_text(face = "bold", hjust = 0.5, size = 18, color = "black"),
plot.subtitle = element_text(hjust = 0.5, size = 12, color = "gray30"),
# Sumbu dan Grid
axis.line = element_line(color = "black", linewidth = 0.8), # Garis sumbu lebih tegas
axis.ticks = element_line(color = "black"),
axis.text = element_text(color = "black"),
panel.grid.major.y = element_line(color = "gray90", linewidth = 0.5), # Tambah grid horizontal tipis
# Legenda
legend.position = "bottom",
legend.title = element_blank(), # Hilangkan judul legenda
legend.text = element_text(size = 12),
legend.key = element_rect(fill = "transparent", color = NA)
)EVAPOTRANSPIRATION
# ==============================================================================
# 2. PERSIAPAN DATA
# ==============================================================================
# Ambil tanggal batas musim tanam untuk garis vertikal
batas_musim <- unique(as.Date(df_simulasi_2018$Tanggal_Tanam_Siklus))
# ==============================================================================
# 3. PLOT ATAS: EVAPOTRANSPIRASI (ET)
# ==============================================================================
p_et <- ggplot(df_simulasi_2018, aes(x = Tanggal)) +
# Garis ETo (Potensial) - Oranye/Coklat
geom_line(aes(y = ETo, color = "ETo (Potensial)"), linetype = "solid", linewidth = 1) +
# Garis ETc (Tanaman) - Biru
geom_line(aes(y = ETc, color = "ETc (Aktual)"), linewidth = 1) +
# Garis Vertikal Pemisah Musim
geom_vline(xintercept = batas_musim, linetype = "dashed", color = "gray40", alpha = 0.7) +
# Label Siklus di Atas
annotate("text", x = batas_musim + 5, y = max(df_simulasi_2018$ETo) * 1.02,
label = paste("Musim Tanam", seq_along(batas_musim)),
hjust = 0, vjust = 0, size = 3.5, family = "serif", color = "black", fontface="bold") +
# Pengaturan Warna (Menyesuaikan contoh gambar: Biru dan Oranye)
scale_color_manual(values = c("ETc (Aktual)" = "#0072B2", "ETo (Potensial)" = "#D55E00")) +
# Axis X
scale_x_date(date_breaks = "1 month", date_labels = "%b") +
# Label dan Tema
labs(title = "Evapotranspirasi — Bogor (2018)",
y = "mm/hari",
x = NULL) + # X label dihapus agar tidak duplikat
theme_publication +
theme(
axis.ticks.x = element_blank()
)
print(p_et)
SUN RADIATION
# ==============================================================================
# 4. PLOT BAWAH: RADIASI & LAMA PENYINARAN
# ==============================================================================
# Faktor Skala untuk sumbu ganda
scale_rad <- 2.5
p_rad <- ggplot(df_simulasi_2018, aes(x = Tanggal)) +
# Area Radiasi (Background) - Pink/Ungu Transparan
geom_area(aes(y = Radiasi, fill = "Radiasi (Rs)"), alpha = 0.4) +
# Garis Lama Penyinaran (Foreground) - Hijau Teal
geom_line(aes(y = Lama_Penyinaran * scale_rad, color = "Lama Penyinaran"), linewidth = 1) +
# Garis Vertikal Pemisah Musim
geom_vline(xintercept = batas_musim, linetype = "dashed", color = "gray9", alpha = 0.7) +
# Label Siklus di Atas
annotate("text", x = batas_musim + 5, y = max(df_simulasi_2018$Radiasi) * 1.02,
label = paste("Musim Tanam", seq_along(batas_musim)),
hjust = 0, vjust = 0, size = 3.5, family = "serif", color = "black", fontface="bold") +
# Dual Axis Y
scale_y_continuous(
name = "Radiasi (MJ/m²/hari)",
sec.axis = sec_axis(~ . / scale_rad, name = "Lama Penyinaran (Jam)")
) +
# Pengaturan Warna & Fill (Menyesuaikan contoh gambar: Pink dan Hijau Teal)
scale_fill_manual(values = c("Radiasi (Rs)" = "#CC79A7")) +
scale_color_manual(values = c("Lama Penyinaran" = "#009E73")) +
# Axis X
scale_x_date(date_breaks = "1 month", date_labels = "%b") +
# Label dan Tema
labs(title = "Lama Penyinaran & Radiasi",
x = "Tanggal (Bulan)") +
theme_publication +
theme(
axis.title.y.right = element_text(vjust = 1, angle = 90) # Merapikan label sumbu kanan
)
print(p_rad)
#KLIMOGRAF 10 TAHUN
# KLIMOGRAF
# Pastikan kolom 'tgl' sudah berformat Date
monthly_clim <- kab_bogor_df %>%
# 1. Buat kolom baru yang membulatkan tanggal ke awal bulan
mutate(bulan = floor_date(tgl, unit = "month")) %>%
# 2. Group berdasarkan bulan tersebut
group_by(bulan) %>%
# 3. Lakukan summarise
# Pastikan nama kolom (CH/Suhu) sesuai besar-kecil hurufnya dengan data aslimu
summarise(CH = sum(ch, na.rm=TRUE), # Curah hujan biasanya di-TOTAL per bulan
Suhu = mean(suhu, na.rm=TRUE)) # Suhu biasanya di-RATA-RATA per bulan
#vis bulanan
# 1. Agregasi data berdasarkan Nama Bulan
# Kita menghitung: "Berapa rata-rata curah hujan di bulan Januari selama 10 tahun terakhir?"
df_seasonal <- monthly_clim %>%
mutate(Bulan_Nama = month(bulan, label = TRUE, abbr = FALSE, locale = "ID")) %>% # Mengambil nama bulan (Januari, dst)
group_by(Bulan_Nama) %>%
summarise(
# Rata-rata CH: Rata-rata dari total hujan per bulan
Rata_CH = mean(CH, na.rm = TRUE),
# Rata-rata Suhu
Rata_Suhu = mean(Suhu, na.rm = TRUE)
)
# Cek hasil (harus ada 12 baris)
print(df_seasonal)## # A tibble: 12 × 3
## Bulan_Nama Rata_CH Rata_Suhu
## <ord> <dbl> <dbl>
## 1 Januari 446. 20.9
## 2 Februari 520. 20.9
## 3 Maret 343. 21.5
## 4 April 370. 21.8
## 5 Mei 231. 22.1
## 6 Juni 158. 21.7
## 7 Juli 128. 21.3
## 8 Agustus 100. 21.3
## 9 September 119. 21.6
## 10 Oktober 248. 21.7
## 11 November 331. 21.7
## 12 Desember 353. 21.5# --- Hitung Faktor Skala Baru ---
# Kita hitung ulang karena nilai rata-rata mungkin lebih rendah/tinggi dari nilai ekstrem bulanan
max_ch_avg <- max(df_seasonal$Rata_CH, na.rm = TRUE)
max_suhu_avg <- max(df_seasonal$Rata_Suhu, na.rm = TRUE)
#vis biar lebih keliatan
# --- 1. Tentukan Batas dan Skala ---
ylim_ch <- max(df_seasonal$Rata_CH, na.rm = TRUE) * 1.1 # Kasih sedikit ruang di atas bar
suhu_min <- 15 # Batas bawah yang kamu minta
suhu_max <- 25 # Batas atas yang kamu minta
# Rumus Skala:
# Kita harus memetakan rentang suhu (10 derajat) ke rentang hujan (0 - ylim_ch)
scale_factor <- ylim_ch / (suhu_max - suhu_min)
ggplot(df_seasonal, aes(x = Bulan_Nama)) +
# --- 2. Grafik Batang (Curah Hujan) - Tetap Normal ---
geom_col(aes(y = Rata_CH), fill = "#80c1ff", alpha = 0.7, width = 0.7) +
# --- 3. Grafik Garis (Suhu) - DITRANSFORMASI ---
# Rumusnya: (Suhu - 20) * scale_factor
# Artinya: Kalau suhunya 20, maka (20-20)*scale = 0 (titik terbawah)
geom_line(aes(y = (Rata_Suhu - suhu_min) * scale_factor, group = 1),
color = "#E31A1C", linewidth = 1.2) +
geom_point(aes(y = (Rata_Suhu - suhu_min) * scale_factor),
color = "#E31A1C", size = 3) +
# Label nilai suhu agar lebih jelas (opsional)
# Label untuk Bar (Curah Hujan)
geom_text(aes(y = Rata_CH, label = round(Rata_CH, 0)),
vjust = -0.5, color = "#2c5aa0", size = 3, fontface = "bold") +
# Label untuk Line (Suhu)
geom_text(aes(y = (Rata_Suhu - suhu_min) * scale_factor,
label = round(Rata_Suhu, 1)),
vjust = -1.5, color = "#E31A1C", size = 3, fontface = "bold") +
# --- 4. Mengatur Sumbu Ganda ---
scale_y_continuous(
name = "Rata-rata Curah Hujan (mm)",
limits = c(0, ylim_ch), # Kunci batas sumbu kiri
# Transformasi Balik untuk Label Kanan:
# (Nilai / scale_factor) + 20
sec.axis = sec_axis(~ (. / scale_factor) + suhu_min,
name = "Rata-rata Suhu (°C)",
breaks = seq(suhu_min, suhu_max, 1)) # Label per 1 derajat
) +
labs(
title = "Pola Iklim Bulanan (2013–2022)",
x = NULL
) +
theme_classic(base_size = 14) +
theme(
text = element_text(family = "serif"),
axis.title.y.left = element_text(color = "#2c5aa0", face = "bold"),
axis.title.y.right = element_text(color = "#E31A1C", face = "bold", angle = 90),
plot.title = element_text(face = "bold", hjust = 0.5),
axis.text.x = element_text(angle = 45, hjust = 1)
)
KLIMOGRAF 2018
# ==============================================================================
# 1. DATA PREPARATION: FILTER KHUSUS 2018
# ==============================================================================
# Kita ambil dari monthly_clim yang sudah dibuat sebelumnya
df_2018 <- monthly_clim %>%
# Filter hanya tahun 2018
filter(year(bulan) == 2018) %>%
# Buat kolom Nama Bulan untuk Sumbu X
mutate(Bulan_Nama = month(bulan, label = TRUE, abbr = FALSE, locale = "ID"))
# Cek data (Pastikan ada 12 baris untuk Jan-Des 2018)
print(df_2018)## # A tibble: 12 × 4
## bulan CH Suhu Bulan_Nama
## <date> <dbl> <dbl> <ord>
## 1 2018-01-01 333. 20.8 Januari
## 2 2018-02-01 672. 20.7 Februari
## 3 2018-03-01 432 21.4 Maret
## 4 2018-04-01 291. 21.9 April
## 5 2018-05-01 108 22.1 Mei
## 6 2018-06-01 152. 21.7 Juni
## 7 2018-07-01 9 21.0 Juli
## 8 2018-08-01 20.5 21.1 Agustus
## 9 2018-09-01 162. 21.3 September
## 10 2018-10-01 130. 21.9 Oktober
## 11 2018-11-01 382. 21.6 November
## 12 2018-12-01 196. 21.6 Desember# ==============================================================================
# 2. HITUNG FAKTOR SKALA (KHUSUS DATA 2018)
# ==============================================================================
# Kita hitung batas atas/bawah berdasarkan data aktual 2018
max_ch_2018 <- max(df_2018$CH, na.rm = TRUE)
# --- Setting Batas Skala ---
ylim_ch <- max_ch_2018 * 1.1 # Memberi ruang kosong 10% di atas bar chart
suhu_min <- 17 # Batas bawah sumbu suhu (sesuaikan jika suhu 2018 ekstrem)
suhu_max <- 25 # Batas atas sumbu suhu (saya naikkan ke 30 untuk jaga-jaga El Nino 2018)
# Rumus Skala Konversi
scale_factor <- ylim_ch / (suhu_max - suhu_min)
# ==============================================================================
# 3. VISUALISASI KLIMOGRAF 2018
# ==============================================================================
clim_2018 <- ggplot(df_2018, aes(x = Bulan_Nama)) +
# --- Grafik Batang (Curah Hujan) ---
# Perhatikan: variable y menggunakan 'CH' (bukan Rata_CH)
geom_col(aes(y = CH), fill = "#80c1ff", alpha = 0.7, width = 0.7) +
# --- Grafik Garis (Suhu) ---
# Perhatikan: variable y menggunakan 'Suhu' (bukan Rata_Suhu)
geom_line(aes(y = (Suhu - suhu_min) * scale_factor, group = 1),
color = "#E31A1C", linewidth = 1.2) +
geom_point(aes(y = (Suhu - suhu_min) * scale_factor),
color = "#E31A1C", size = 3) +
# --- Label Angka ---
# Label CH
geom_text(aes(y = CH, label = round(CH, 0)),
vjust = -0.5, color = "#2c5aa0", size = 3, fontface = "bold") +
# Label Suhu
geom_text(aes(y = (Suhu - suhu_min) * scale_factor,
label = round(Suhu, 1)),
vjust = -1.5, color = "#E31A1C", size = 3, fontface = "bold") +
# --- Sumbu Ganda (Dual Axis) ---
scale_y_continuous(
name = "Curah Hujan (mm)",
limits = c(0, ylim_ch),
sec.axis = sec_axis(~ (. / scale_factor) + suhu_min,
name = "Suhu (°C)",
breaks = seq(suhu_min, suhu_max, 1))
) +
# --- Judul & Tema ---
labs(
title = "Klimograf Bogor Tahun 2018",
subtitle = "Total Curah Hujan & Rata-rata Suhu Bulanan",
x = NULL
) +
theme_classic(base_size = 14) +
theme(
text = element_text(family = "serif"),
axis.title.y.left = element_text(color = "#2c5aa0", face = "bold"),
axis.title.y.right = element_text(color = "#E31A1C", face = "bold", angle = 90),
plot.title = element_text(face = "bold", hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5, color = "gray40"),
axis.text.x = element_text(angle = 45, hjust = 1)
)
print(clim_2018)
PENGARUH SUHU
# ==============================================================================
# 1. PERSIAPAN DATA (Filter Tahun 2018)
# ==============================================================================
# Menggunakan df_simulasi_start2018 karena struktur datanya paling lengkap
# Kita filter hanya tahun 2018 dan pastikan tidak ada nilai NA
mod_df <- df_simulasi_2018
# ==============================================================================
# 2. PEMODELAN & VISUALISASI
# ==============================================================================
if(nrow(mod_df) >= 30) {
# --- A. Membuat Model Regresi Linier ---
# Model: Produktivitas diprediksi oleh Suhu harian
lm_mod <- lm(Prod1HAinTon ~ Suhu, data = mod_df)
# Tampilkan ringkasan model di console (opsional, untuk cek signifikansi)
print(summary(lm_mod))
# --- B. Membuat Skenario Prediksi ---
base_mean_temp <- mean(mod_df$Suhu, na.rm=TRUE)
scenarios <- data.frame(
Suhu = c(base_mean_temp, base_mean_temp+1, base_mean_temp+2, base_mean_temp+3),
scenario = c("Base", "+1°C", "+2°C", "+3°C")
)
# Hitung prediksi berdasarkan model
scenarios$predicted_prod <- predict(lm_mod, newdata = scenarios)
# --- C. Mengatur Urutan Factor (Agar grafik rapi) ---
scenarios$scenario <- factor(scenarios$scenario,
levels = c("Base", "+1°C", "+2°C", "+3°C"))
# --- D. Visualisasi ---
p10 <- ggplot(scenarios, aes(x = scenario, y = predicted_prod, fill = scenario)) +
# Membuat Bar Chart
geom_col(width = 0.6, alpha = 0.9, color = "black") +
# Menambahkan Label Angka di Atas Batang
geom_text(aes(label = round(predicted_prod, 2)),
vjust = -0.5, fontface = "bold", family = "serif", size = 4) +
# Pewarnaan: Kuning (Base) ke Merah Tua (+3°C)
scale_fill_manual(values = c("Base" = "#F0E442",
"+1°C" = "#E69F00",
"+2°C" = "#D55E00",
"+3°C" = "#800000")) +
# Skala Y: Memberi ruang di atas bar
scale_y_continuous(expand = expansion(mult = c(0, 0.2))) +
# Label dan Judul
labs(
title = "Proyeksi Produktivitas Timun vs Kenaikan Suhu",
subtitle = paste("Model Regresi Linier Data Simulasi Tahun 2018 (n =", nrow(mod_df), ")"),
x = "Skenario Kenaikan Suhu",
y = "Prediksi Rata-rata Produktivitas Harian (Ton/Ha)",
fill = "Skenario"
) +
# Tema Konsisten
theme_classic(base_size = 14) +
theme(
text = element_text(family = "serif"),
plot.title = element_text(face = "bold", hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5, color = "gray40"),
axis.text.y = element_text(color = "black"),
legend.position = "none"
)
print(p10)
} else {
message("Data tahun 2015 tidak mencukupi untuk analisis regresi (< 30 observasi).")
}##
## Call:
## lm(formula = Prod1HAinTon ~ Suhu, data = mod_df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.7608 -9.7884 0.3107 8.2014 18.0419
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -13.1998 12.6374 -1.045 0.2970
## Suhu 1.0908 0.5895 1.850 0.0651 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.365 on 352 degrees of freedom
## Multiple R-squared: 0.009632, Adjusted R-squared: 0.006818
## F-statistic: 3.423 on 1 and 352 DF, p-value: 0.06512
VALIDASI MODEL
# ==============================================================================
# 1. INPUT DATA AKTUAL (BPS KAB. BOGOR)
# ==============================================================================
# Data diketik manual berdasarkan gambar Excel yang Anda lampirkan
# Satuan awal: Kuintal/Hektar -> Dikonversi ke Ton/Ha (dibagi 10)
df_simulasi_start2018 <- read_excel("df_simulasi_start2018.xlsx")
df_bps <- data.frame(
Tahun = c(2018, 2019, 2020, 2021, 2022),
Prod_Kuintal = c(157.04, 152.24, 151.57, 201.61, 243.83)
) %>%
mutate(
Prod_Aktual = Prod_Kuintal / 10, # Konversi ke Ton/Ha
Sumber = "Observasi (BPS)"
)
print("Data Aktual BPS (Converted to Ton/Ha):")## [1] "Data Aktual BPS (Converted to Ton/Ha):"print(df_bps)## Tahun Prod_Kuintal Prod_Aktual Sumber
## 1 2018 157.04 15.704 Observasi (BPS)
## 2 2019 152.24 15.224 Observasi (BPS)
## 3 2020 151.57 15.157 Observasi (BPS)
## 4 2021 201.61 20.161 Observasi (BPS)
## 5 2022 243.83 24.383 Observasi (BPS)# ==============================================================================
# 2. PERSIAPAN DATA MODEL (AGREGASI TAHUNAN)
# ==============================================================================
# Karena data BPS adalah data tahunan, kita harus merata-ratakan hasil
# dari 4 musim tanam dalam satu tahun di model kita.
df_model_validasi <- df_simulasi_start2018 %>%
# Filter tahun yang sama dengan data BPS
filter(Tahun >= 2018 & Tahun <= 2022) %>%
# Langkah 1: Ambil Panen Maksimum per Musim Tanam
group_by(Tahun, Musim_Ke) %>%
summarise(Panen_Per_Musim = max(Prod1HAinTon, na.rm=TRUE), .groups = "drop") %>%
# Langkah 2: Rata-ratakan 4 musim tersebut menjadi 1 nilai per tahun
group_by(Tahun) %>%
summarise(
Prod_Model_Mean = mean(Panen_Per_Musim, na.rm=TRUE),
Prod_Model_Min = min(Panen_Per_Musim, na.rm=TRUE), # Untuk Error Bar
Prod_Model_Max = max(Panen_Per_Musim, na.rm=TRUE) # Untuk Error Bar
) %>%
mutate(Sumber = "Simulasi Model")# ==============================================================================
# 3. PENGGABUNGAN & PERHITUNGAN ERROR (RMSE & MAPE)
# ==============================================================================
df_validasi <- left_join(df_bps, df_model_validasi, by = "Tahun")
# Hitung Statistik Error
rmse <- sqrt(mean((df_validasi$Prod_Model_Mean - df_validasi$Prod_Aktual)^2))
mape <- mean(abs((df_validasi$Prod_Model_Mean - df_validasi$Prod_Aktual) / df_validasi$Prod_Aktual)) * 100
korelasi <- cor(df_validasi$Prod_Model_Mean, df_validasi$Prod_Aktual)
print(paste("RMSE (Root Mean Square Error):", round(rmse, 2), "Ton/Ha"))## [1] "RMSE (Root Mean Square Error): 4.91 Ton/Ha"print(paste("MAPE (Mean Absolute Percentage Error):", round(mape, 2), "%"))## [1] "MAPE (Mean Absolute Percentage Error): 26.43 %"print(paste("Korelasi (r):", round(korelasi, 2)))## [1] "Korelasi (r): -0.72"# ==============================================================================
# 4. VISUALISASI VALIDASI MODEL (GGPLOT)
# ==============================================================================
# Menyiapkan label performa untuk ditempel di grafik
label_stats <- paste0(
"Evaluasi Model (2018-2022):\n",
"RMSE: ", round(rmse, 2), " Ton/Ha\n",
"MAPE: ", round(mape, 2), " %"
)
p_val <- ggplot() +
# --- A. AREA RANGE MODEL (Min-Max 4 Musim) ---
# Ini menunjukkan bahwa meskipun rata-rata mungkin beda, data BPS
# mungkin masih masuk dalam "Range" hasil simulasi musim basah/kering
geom_ribbon(data = df_model_validasi,
aes(x = Tahun, ymin = Prod_Model_Min, ymax = Prod_Model_Max, fill = "Rentang Model (4 Musim)"),
alpha = 0.2) +
# --- B. GARIS MODEL (Rata-rata) ---
geom_line(data = df_model_validasi,
aes(x = Tahun, y = Prod_Model_Mean, color = "Simulasi Model (Rata-rata)"),
linewidth = 1.2, linetype = "dashed") +
geom_point(data = df_model_validasi,
aes(x = Tahun, y = Prod_Model_Mean, color = "Simulasi Model (Rata-rata)"),
size = 3) +
# --- C. GARIS OBSERVASI (BPS) ---
geom_line(data = df_validasi,
aes(x = Tahun, y = Prod_Aktual, color = "Observasi Aktual (BPS)"),
linewidth = 1.2) +
geom_point(data = df_validasi,
aes(x = Tahun, y = Prod_Aktual, color = "Observasi Aktual (BPS)"),
size = 3, shape = 17) + # Shape segitiga untuk pembeda
# --- D. LABEL ANGKA ---
geom_text(data = df_validasi, aes(x = Tahun, y = Prod_Aktual, label = round(Prod_Aktual, 1)),
vjust = -1, color = "#E31A1C", fontface = "bold", size = 3.5) +
# --- E. KOSMETIK & TEMA ---
scale_fill_manual(values = c("Rentang Model (4 Musim)" = "gray50"), name = NULL) +
scale_color_manual(values = c("Simulasi Model (Rata-rata)" = "black",
"Observasi Aktual (BPS)" = "#E31A1C"), name = NULL) +
scale_y_continuous(limits = c(0, 30)) + # Sesuaikan batas Y agar lega
labs(
title = "Validasi Model: Simulasi vs Data Aktual BPS",
subtitle = "Perbandingan Produktivitas Timun di Kab. Bogor (2018-2022)",
y = "Produktivitas (Ton/Ha)",
x = "Tahun"
) +
# Menambahkan Kotak Statistik di Pojok
annotate("label", x = 2018, y = 28, label = label_stats,
hjust = 0, vjust = 1, family = "serif", size = 3.5,
fill = "white", alpha = 0.8, color = "black") +
theme_classic(base_size = 14) +
theme(
text = element_text(family = "serif"),
plot.title = element_text(face = "bold", hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5, color = "gray40"),
legend.position = "bottom",
legend.box = "vertical"
)
print(p_val)
COMPARE PRODUKTIVITAS
# ==============================================================================
# 1. PERSIAPAN DATA BPS (Manual Input dari Gambar Excel)
# ==============================================================================
# Input data BPS (2018-2022)
# Satuan Asli: Kuintal/Hektar. Kita konversi ke Ton/Ha (bagi 10)
df_bps## Tahun Prod_Kuintal Prod_Aktual Sumber
## 1 2018 157.04 15.704 Observasi (BPS)
## 2 2019 152.24 15.224 Observasi (BPS)
## 3 2020 151.57 15.157 Observasi (BPS)
## 4 2021 201.61 20.161 Observasi (BPS)
## 5 2022 243.83 24.383 Observasi (BPS)# ==============================================================================
# 2. AGREGASI DATA MODEL SIMULASI
# ==============================================================================
# Kita ambil dari df_simulasi_start2018
df_model_compare <- df_simulasi_start2018 %>%
# Filter tahun yang beririsan dengan data BPS
filter(Tahun >= 2018 & Tahun <= 2022) %>%
# Ambil nilai panen maksimum dari setiap musim tanam
group_by(Tahun, Musim_Ke) %>%
summarise(Panen_Musim = max(Prod1HAinTon, na.rm=TRUE), .groups = "drop") %>%
# Rata-ratakan 4 musim menjadi 1 nilai per tahun
group_by(Tahun) %>%
summarise(
Prod_Model = mean(Panen_Musim, na.rm=TRUE),
Prod_Model_SD = sd(Panen_Musim, na.rm=TRUE) # Standar Deviasi (Variasi musim)
) %>%
mutate(Sumber = "Model Simulasi (Potensial)")
# ==============================================================================
# 3. GABUNGKAN DATA (MERGE)
# ==============================================================================
df_final_compare <- left_join(df_bps, df_model_compare, by = "Tahun") %>%
mutate(
# Hitung Selisih (Gap)
Gap = Prod_Model - Prod_Aktual,
# Status: Apakah Model Overestimate (Lebih tinggi) atau Underestimate
Status = ifelse(Gap > 0, "Model > Aktual (Yield Gap)", "Model < Aktual")
)
print("Tabel Perbandingan Tahunan:")## [1] "Tabel Perbandingan Tahunan:"print(df_final_compare %>% select(Tahun, Prod_Aktual, Prod_Model, Gap, Status))## Tahun Prod_Aktual Prod_Model Gap Status
## 1 2018 15.704 21.74980 6.0457978 Model > Aktual (Yield Gap)
## 2 2019 15.224 21.10946 5.8854639 Model > Aktual (Yield Gap)
## 3 2020 15.157 19.65724 4.5002413 Model > Aktual (Yield Gap)
## 4 2021 20.161 19.50451 -0.6564874 Model < Aktual
## 5 2022 24.383 19.00766 -5.3753410 Model < Aktual# ==============================================================================
# 4. VISUALISASI PERBANDINGAN
# ==============================================================================
# Definisikan warna agar konsisten
cols <- c("Data BPS (Aktual)" = "#E31A1C", "Model Simulasi (Potensial)" = "#1f78b4")
ggplot(df_final_compare, aes(x = Tahun)) +
# --- 1. Area Error Bar Model (Variasi Musiman) ---
# Menunjukkan bahwa dalam 1 tahun, hasil panen bisa bervariasi tergantung musim
geom_errorbar(aes(ymin = Prod_Model - Prod_Model_SD,
ymax = Prod_Model + Prod_Model_SD,
color = "Model Simulasi (Potensial)"),
width = 0.1, linewidth = 0.8) +
# --- 2. Garis Model ---
geom_line(aes(y = Prod_Model, color = "Model Simulasi (Potensial)"),
linewidth = 1.2, linetype = "dashed") +
geom_point(aes(y = Prod_Model, color = "Model Simulasi (Potensial)"),
size = 3.5, shape = 16) +
# --- 3. Garis BPS ---
geom_line(aes(y = Prod_Aktual, color = "Data BPS (Aktual)"),
linewidth = 1.2) +
geom_point(aes(y = Prod_Aktual, color = "Data BPS (Aktual)"),
size = 3.5, shape = 17) +
# --- 4. Label Angka di Grafik ---
geom_text(aes(y = Prod_Aktual, label = round(Prod_Aktual, 1)),
vjust = -1, color = "#E31A1C", fontface = "bold", size = 3.5) +
geom_text(aes(y = Prod_Model, label = round(Prod_Model, 1)),
vjust = 1.5, color = "#1f78b4", fontface = "bold", size = 3.5) +
# --- 5. Kosmetik Grafik ---
scale_y_continuous(limits = c(0, max(df_final_compare$Prod_Aktual, df_final_compare$Prod_Model) + 5)) +
scale_color_manual(name = NULL, values = cols) +
labs(
title = "Komparasi Produktivitas: Model vs Aktual",
subtitle = "Analisis Yield Gap (Celah Hasil) di Kab. Bogor 2018-2022",
y = "Produktivitas (Ton/Ha)",
x = "Tahun",
caption = "Error Bar pada garis biru menunjukkan variasi hasil antar musim dalam tahun tersebut"
) +
theme_classic(base_size = 14) +
theme(
text = element_text(family = "serif"),
plot.title = element_text(face = "bold", hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5, color = "gray40"),
legend.position = "bottom"
)
BAR CHART VER.
# ==============================================================================
# 1. PERSIAPAN DATA
# ==============================================================================
# --- A. Data BPS (Aktual) ---
df_bps_bar <- data.frame(
Tahun = c(2018, 2019, 2020, 2021, 2022),
Nilai = c(157.04, 152.24, 151.57, 201.61, 243.83)
) %>%
mutate(
Nilai = Nilai / 10, # Konversi Kuintal ke Ton
Sumber = "Data Aktual (BPS)",
SD = NA # Tidak ada variasi musim untuk data tahunan BPS
)
# --- B. Data Model (Simulasi) ---
df_model_bar <- df_simulasi_start2018%>%
filter(Tahun >= 2018 & Tahun <= 2022) %>%
# Ambil max panen per musim, lalu rata-ratakan per tahun
group_by(Tahun, Musim_Ke) %>%
summarise(Panen_Musim = max(Prod1HAinTon, na.rm=TRUE), .groups = "drop") %>%
group_by(Tahun) %>%
summarise(
Nilai = mean(Panen_Musim, na.rm=TRUE),
SD = sd(Panen_Musim, na.rm=TRUE) # Standar Deviasi antar musim
) %>%
mutate(Sumber = "Model Simulasi")
# --- C. Gabung Data ---
df_gabung <- bind_rows(df_bps_bar, df_model_bar)
# Urutkan Level agar Hijau (Aktual) muncul duluan atau sesuai keinginan
df_gabung$Sumber <- factor(df_gabung$Sumber, levels = c("Data Aktual (BPS)", "Model Simulasi"))
# ==============================================================================
# 2. VISUALISASI (HIJAU & KUNING)
# ==============================================================================
p_compare <- ggplot(df_gabung, aes(x = factor(Tahun), y = Nilai, fill = Sumber)) +
# --- 1. Bar Chart Samping-Sampingan (Dodge) ---
geom_col(position = position_dodge(width = 0.8), width = 0.7, color = "black", alpha = 0.9) +
# --- 2. Error Bar (Untuk Model) ---
geom_errorbar(aes(ymin = Nilai - SD, ymax = Nilai + SD),
position = position_dodge(width = 0.8),
width = 0.25, size = 0.6, color = "black") +
# --- 3. Label Angka ---
geom_text(aes(label = round(Nilai, 1)),
position = position_dodge(width = 0.8),
vjust = -0.5, size = 4, fontface = "bold", family = "serif") +
# --- 4. Warna Kustom (Hijau & Kuning) ---
scale_fill_manual(values = c(
"Data Aktual (BPS)" = "#228B22", # ForestGreen (Hijau Tua yang enak dilihat)
"Model Simulasi" = "#FFD700" # Gold (Kuning Emas agar kontras dengan putih)
)) +
# --- 5. Pengaturan Sumbu & Skala ---
scale_y_continuous(limits = c(0, max(df_gabung$Nilai) * 1.15), expand = c(0,0)) +
# --- 6. Label & Tema ---
labs(
title = "Validasi Model: Aktual vs Simulasi (2018-2022)",
subtitle = "Hijau: Observasi Lapangan | Kuning: Potensi Model Klimatologis",
y = "Produktivitas (Ton/Ha)",
x = "Tahun",
fill = NULL
) +
theme_classic(base_size = 14) +
theme(
text = element_text(family = "serif"),
plot.title = element_text(face = "bold", hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5, color = "gray40"),
legend.position = "bottom",
axis.line = element_line(linewidth = 0.8),
panel.grid.major.y = element_line(color = "gray90", linetype = "dashed")
)
print(p_compare)
ADA RMSE-NYA
library(dplyr)
library(ggplot2)
# ==============================================================================
# 1. PERSIAPAN DATA
# ==============================================================================
# --- A. Data BPS (Aktual) ---
df_bps_bar <- data.frame(
Tahun = c(2018, 2019, 2020, 2021, 2022),
Nilai = c(157.04, 152.24, 151.57, 201.61, 243.83)
) %>%
mutate(
Nilai = Nilai / 10, # Konversi Kuintal ke Ton
Sumber = "Data Aktual (BPS)",
SD = NA
)
# --- B. Data Model (Simulasi) ---
df_model_bar <- df_simulasi_start2018 %>%
filter(Tahun >= 2018 & Tahun <= 2022) %>%
group_by(Tahun, Musim_Ke) %>%
summarise(Panen_Musim = max(Prod1HAinTon, na.rm=TRUE), .groups = "drop") %>%
group_by(Tahun) %>%
summarise(
Nilai = mean(Panen_Musim, na.rm=TRUE),
SD = sd(Panen_Musim, na.rm=TRUE)
) %>%
mutate(Sumber = "Model Simulasi")
# --- C. Hitung RMSE & MAPE (BAGIAN BARU) ---
# Kita gabungkan dulu kedua tabel untuk menghitung selisihnya
df_metric <- left_join(df_bps_bar, df_model_bar, by = "Tahun", suffix = c("_bps", "_model"))
# Rumus RMSE & MAPE
rmse_val <- sqrt(mean((df_metric$Nilai_model - df_metric$Nilai_bps)^2))
mape_val <- mean(abs((df_metric$Nilai_model - df_metric$Nilai_bps) / df_metric$Nilai_bps)) * 100
# Buat Label Teks untuk Grafik
label_validasi <- paste0(
"Statistik Validasi:\n",
"RMSE : ", round(rmse_val, 2), " Ton/Ha\n",
"MAPE : ", round(mape_val, 2), " %"
)
# --- D. Gabung Data untuk Plotting ---
df_gabung <- bind_rows(df_bps_bar, df_model_bar)
df_gabung$Sumber <- factor(df_gabung$Sumber, levels = c("Data Aktual (BPS)", "Model Simulasi"))
# ==============================================================================
# 2. VISUALISASI DENGAN LABEL METRIK
# ==============================================================================
p_compare <- ggplot(df_gabung, aes(x = factor(Tahun), y = Nilai, fill = Sumber)) +
# --- 1. Bar Chart ---
geom_col(position = position_dodge(width = 0.8), width = 0.7, color = "black", alpha = 0.9) +
# --- 2. Error Bar ---
geom_errorbar(aes(ymin = Nilai - SD, ymax = Nilai + SD),
position = position_dodge(width = 0.8),
width = 0.25, size = 0.6, color = "black") +
# --- 3. Label Angka di Bar ---
geom_text(aes(label = round(Nilai, 1)),
position = position_dodge(width = 0.8),
vjust = -0.5, size = 4, fontface = "bold", family = "serif") +
# --- 4. KOTAK RMSE & MAPE (BAGIAN BARU) ---
# Menempatkan kotak di pojok kiri atas
annotate("label",
x = 0.6, # Posisi X (sedikit di kiri tahun pertama)
y = max(df_gabung$Nilai) * 1.1, # Posisi Y (di atas bar tertinggi)
label = label_validasi,
hjust = 0, vjust = 1, # Rata kiri atas
family = "serif", size = 3.5,
fill = "white", alpha = 0.8, color = "black") +
# --- 5. Warna & Tema ---
scale_fill_manual(values = c(
"Data Aktual (BPS)" = "#228B22",
"Model Simulasi" = "#FFD700"
)) +
scale_y_continuous(limits = c(0, max(df_gabung$Nilai) * 1.15), expand = c(0,0)) +
labs(
title = "Validasi Model: Aktual vs Simulasi (2018-2022)",
subtitle = "Perbandingan Hasil Lapangan dengan Potensi Iklim",
y = "Produktivitas (Ton/Ha)",
x = "Tahun",
fill = NULL
) +
theme_classic(base_size = 14) +
theme(
text = element_text(family = "serif"),
plot.title = element_text(face = "bold", hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5, color = "gray40"),
legend.position = "bottom",
axis.line = element_line(linewidth = 0.8),
panel.grid.major.y = element_line(color = "gray90", linetype = "dashed")
)
print(p_compare)
R Markdown
This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.
When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:
summary(cars)## speed dist
## Min. : 4.0 Min. : 2.00
## 1st Qu.:12.0 1st Qu.: 26.00
## Median :15.0 Median : 36.00
## Mean :15.4 Mean : 42.98
## 3rd Qu.:19.0 3rd Qu.: 56.00
## Max. :25.0 Max. :120.00Including Plots
You can also embed plots, for example:
Note that the echo = FALSE parameter was added to the
code chunk to prevent printing of the R code that generated the
plot.