November Rain merupakan satu dari banyak lagu yang akrab di telingaku. Tak hanya itu, lagu ini juga biasanya aku gunakan sebagai pengingat waktu. November biasanya menjadi awal musim penghujan. Musim penghujan tersebut pada umumnya terjadi pada bulan November hingga Maret karena adanya angin muson barat yang membawa banyak uap air ke sebagian besar wilayah Indonesia.
Musim kemarau biasanya terjadi pada bulan Mei hingga September. Hal ini dikarenakan adanya angin muson timur yang sifatnya kering dan berhembus ke sebagian besar wilayah Indonesia. Meskipun demikian, pola musim penghujan dan kemarau tersebut dapat berbeda karena banyak hal, misalnya adalah posisi wilayahnya dan adanya pengaruh angin lokal.
Untuk itu, melalui artikel ini aku mengajak pembaca untuk melihat dan menganalisis curah hujan di Indonesia. Nanti kita akan menganalisis curah hujan di Indonesia berdasarkan waktu dan posisi geografisnya dengan menggunakan R. Untuk mengawalinya, kita perlu memanggil beberapa paket yang diperlukan.
Paket {pRecipe} akan kita gunakan untuk mengunduh dan mengimpor data presipitasi atau curah hujan. Paket {giscoR} kita pakai untuk mengunduh data dari pangkalan data Eurostat GISCO (Geographic Information System of the Commission). Secara khusus, kita menggunakan paket ini untuk mendapatkan data peta Indonesia. Kita gunakan paket {terra} untuk melakukan analisis data spasial. Seperti biasa, paket {tidyverse} kita gunakan untuk menganalisis dan memvisualisasikan data. Paket {rayshader} dapat membantu kita dalam pembuatan peta dan visualisasi data dalam 2D dan 3D. Terakhir, paket {sf} memungkinkan kita untuk mengolah data spasial yang berupa vektor.
Mempersiapkan Data
Data yang perlu kita miliki adalah data curah hujan. Data tersebut dapat kita peroleh dengan menggunakan fungsi download_data() dari {pRecipe}.
Ada empat argumen yang kita inputkan dalam fungsi download_data() di atas. Argumen dataset = "mswep" mengindikasikan bahwa set data yang kita unduh adalah data MSWEP (Multi-Source Weighted-Ensemble Precipitation). Argumen path = getwd() menunjukkan tempat penyimpanan hasil unduhannya. getwd() mengindikasikan bahwa tempat penyimpanannya adalah folder kerja kita. Argumen domain = "raw" menunjukkan domain set data yang kita inginkan. Nilai domain tersebut juga dapat diisi dengan “global”, “land”, atau “ocean”. Terakhir, timestep = "monthly" menerangkan bahwa set data yang kita inginkan memiliki interval waktu bulanan.
Setelah kita jalankan kode di atas, kita mendapatkan sebuah fail dengan nama “mswep_tp_mm_global_197902_202301_025_monthly.nc”. Fail tersebut berisi data curah hujan untuk semua wilayah di dunia. Padahal, kita hanya memerlukan data untuk Indonesia saja. Oleh karena itu, kita perlu memotong data tersebut. Untuk melakukan hal ini, kita memerlukan poligon wilayah Indonesia. Bagaimana mendapatkannya? Kita dapat memerolehnya dengan menggunakan fungsi gisco_get_countries() dari paket {giscoR}.
Sekarang kita buat sebuah SpatRaster dari fail “mswep_tp_mm_global_197902_202301_025_monthly.nc” dengan menggunakan fungsi rast() dari {terra}. Setelah itu, kita potong hasilnya dengan sf_negara menggunakan fungsi crop() dari {terra}. Kita namai hasilnya sebagai data_mswep.
Jika kita menjalankan perintah print(data_mswep), kita mengetahui bahwa objek data_mswep tersebut memiliki 68 baris, 184 kolom, dan 528 lapisan/layer. Baris dan kolom tersebut menyatakan koordinat \(x\) dan \(y\) wilayah Indonesia sedangkan lapisannya merupakan nilai-nilai variabel precipitation atau curah hujan di tiap-tiap koordinat tersebut. Banyaknya lapisan adalah 528 karena nilai precipitation untuk setiap koordinat wilayahnya merentang dari 1979-02-01 sampai 2023-01-01 (528 bulan).
Pada bagian names, kita dapat melihat nama-nama variabelnya adalah precipitation_1, precipitation_2, precipitation_3, dan seterusnya sampai precipitation_528. Nama-nama ini kurang informatif. Karena nama-nama itu merepresentasikan curah hujan dari 1979-02-01 sampai 2023-01-01, kita dapat menggantinya dengan barisan tanggal itu.
barisan_tanggal <- as.character(
seq(as.Date("1979-02-01"), as.Date("2023-01-01"), by = "month")
)
names(data_mswep) <- barisan_tanggalAgar memudahkan kita dalam mengolah data, kita ubah objek data_mswep menjadi sebuah data frame. Untuk melakukannya, kita dapat menggunakan fungsi as.data.frame() dari paket {terra}. Kita inputkan xy = TRUE ke dalam fungsi tersebut agar mempertahankan koordinat \(x\) dan \(y\). Setelah itu, kita juga dapat mengubah hasilnya menjadi sebuah tibble dengan fungsi as_tibble() agar mudah terbaca. Kita namai hasilnya dengan df_mswep.
Jika kita lihat hasilnya dengan print(df_mswep), ternyata df_mswep memiliki 12.512 baris dan 530 kolom! Kita buat data tersebut menjadi data yang memanjang dengan fungsi pivot_longer(). Kolom-kolom yang kita putar adalah semua kolom kecuali kolom x dan y sehingga kita inputkan argumen cols = !c("x", "y"). Kita namai variabel yang memuat kolom-kolom dari data sebelumnya dengan tanggal_bawah dan kita namai variabel yang berisi nilai-nilainya dengan curah_hujan. Agar nilai-nilai dalam variabel tanggal_bawah memiliki kelas date, kita modifikasi variabel ini dengan mutate() dan as.Date(). Untuk melihat hasilnya, kita dapat memanggil fungsi print() terhadap df_mswep.
df_mswep <- df_mswep |>
pivot_longer(
!c("x", "y"),
names_to = "tanggal_bawah",
values_to = "curah_hujan"
) |>
mutate(
tanggal_bawah = as.Date(tanggal_bawah)
)Data df_mswep telah rapi. Data tersebut memuat 6.606.336 baris dan 4 variabel/kolom. Keempat variabel tersebut adalah x, y, tanggal_bawah, dan curah_hujan. Data ini telah siap untuk kita eksplorasi pada bagian berikutnya.
Bagian 1 di atas menghasilkan sf_negara, data_mswep, dan df_mswep. Ketiga data tersebut sudah aku siapkan di awan. Untuk itu, kamu dapat memuatnya langsung ke lembar kerja R dengan kode berikut.
Menjelajah Data
Data yang telah kita peroleh pada bagian Bagian 1 dapat kita gunakan untuk melihat tren rerata curah hujan per bulannya, mulai Februari 1979 sampai Januari 2023. Untuk melakukannya, kita perlu mengelompokkan baris-baris dalam df_mswep berdasarkan nilai dalam variabel tanggal_bawah. Pengelompokkan ini dapat kita lakukan dengan menggunakan fungsi group_by(). Setelah itu, kita hitung rerata curah_hujan pada setiap kelompoknya dengan menggunakan fungsi summarise().
Hasil pada proses sebelumnya kemudian dapat gunakan untuk membuat diagram garis. Perhatikan Gambar 1!
df_mswep |>
group_by(tanggal_bawah) |>
summarise(
rerata_curah_hujan = mean(curah_hujan, na.rm = TRUE),
.groups = "drop"
) |>
ggplot(
aes(x = tanggal_bawah, y = rerata_curah_hujan)
) +
geom_line() +
theme_minimal() +
labs(
x = "Waktu",
y = "Rerata curah hujan (mm)"
)
Apakah kita dapat membandingkan rerata curah hujan setiap bulannya? Pertanyaan ini dapat kita jawab dengan membandingkan distribusi curah hujan untuk bulan Januari sampai Desember. Hal ini dapat kita lakukan dengan diagram kotak. Perhatikan Gambar 2!
vek_nama_bulan <- c("Januari", "Februari", "Maret", "April", "Mei", "Juni", "Juli", "Agustus", "September", "Oktober", "November", "Desember")
df_mswep |>
mutate(
bulan = month(tanggal_bawah),
nama_bulan = vek_nama_bulan[month(tanggal_bawah)]
) |>
ggplot(
aes(x = fct_reorder(nama_bulan, bulan), y = curah_hujan)
) +
geom_boxplot(outliers = FALSE) +
theme_minimal() +
labs(
x = "Bulan",
y = "Curah hujan (mm)"
)
Berdasarkan Gambar 2, kita dapat melihat bahwa bulan Januari dan Desember memiliki nilai tengah (median) curah hujan yang paling tinggi dibandingkan dengan bulan-bulan lainnya. Sebaliknya, bulan Agustus dan September memiliki curah hujan yang relatif lebih rendah dibandingkan dengan yang lain.
Pada Bagian 2 ini kita telah menjelajah data curah hujan di Indonesia. Melalui penjelajahan tersebut, kita telah melihat tren rerata curah hujan di Indonesia dari waktu ke waktu. Selain itu, kita juga telah dapat membandingkan rerata curah hujan setiap bulannya. Meskipun demikian, dalam penjelajahan tersebut kita belum mempertimbangkan perbedaan curah hujan antara satu wilayah dengan wilayah lainnya. Kita belum memanfaatkan variabel x dan y dalam data df_mswep.
Visualisasi Data dengan Peta 2D
Curah hujan tergantung dari wilayahnya. Suatu wilayah dapat memiliki curah hujan yang relatif lebih tinggi dibandingkan wilayah lainnya. Demikian juga sebaliknya. Untuk itu, kita perlu melakukan visualisasi data yang berbeda dari apa yang kita lakukan pada Bagian 2. Kita akan memanfaatkan peta.
Peta pertama yang kita buat adalah peta yang menunjukkan rerata curah hujan di setiap wilayah Indonesia empat tahun terakhir, yaitu tahun 2019 – 2022. (Data pada tahun 2023 kita abaikan karena hanya memuat data bulan Januari.) Untuk itu, kita siapkan terlebih dahulu datanya. Data tersebut kita beri nama df_mswep_19_22.
# Data curah hujan tahunan
df_mswep_tahunan <- df_mswep |>
filter(
tanggal_bawah >= as.Date("1980-01-01"),
tanggal_bawah <= as.Date("2022-12-01")
) |>
mutate(tahun = year(tanggal_bawah)) |>
group_by(x, y, tahun) |>
summarise(
curah_hujan_tahunan = sum(curah_hujan, na.rm = TRUE),
.groups = "drop"
)
# Data curah hujan tahunan pada tahun 2019 -- 2022
df_mswep_19_22 <- df_mswep_tahunan |>
filter(
tahun %in% c(2019, 2020, 2021, 2022)
)
print(df_mswep_19_22)# A tibble: 50,048 × 4
x y tahun curah_hujan_tahunan
<dbl> <dbl> <dbl> <dbl>
1 95.1 -10.9 2019 1056.
2 95.1 -10.9 2020 1966.
3 95.1 -10.9 2021 3029.
4 95.1 -10.9 2022 2868.
5 95.1 -10.6 2019 1071.
6 95.1 -10.6 2020 2034.
7 95.1 -10.6 2021 2997.
8 95.1 -10.6 2022 3136.
9 95.1 -10.4 2019 1029.
10 95.1 -10.4 2020 2097.
# ℹ 50,038 more rowsKita gunakan df_mswep_19_22 untuk membuat peta. Terdapat tiga fungsi utama yang perlu kita gunakan:
geom_raster()untuk memplot kotak-kotak yang warnanya mengikuti nilai rerata curah hujan di setiap koordinat kotak-kotak tersebut;geom_contour()untuk memvisualisasikan kontur (garis bentuk); dangeom_sf()untuk memvisualisasikan objek-objek sf (misalnya titik, garis, atau poligon).
Peta pertama kita dapat dibuat dengan menggunakan baris kode seperti berikut. Hasilnya ditunjukkan pada Gambar 3.
df_mswep_19_22 |>
ggplot() +
geom_raster(
aes(
x = x,
y = y,
fill = curah_hujan_tahunan
)
) +
geom_contour(
aes(
x = x,
y = y,
z = curah_hujan_tahunan
),
color = "white"
) +
geom_sf(
data = sf_negara,
fill = "transparent",
color = "grey10",
size = .5
) +
facet_wrap(vars(tahun), ncol = 2) +
theme_minimal()
Kita dapat memperbaiki tampilan pada Gambar 3. Perbaikan tersebut dapat kita lakukan pada tema, interval, dan warnanya. Kita siapkan ketiga hal tersebut.
tema_peta <- function(){
theme_minimal() +
theme(
axis.line = element_blank(),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
legend.position = "bottom",
panel.grid.major = element_line(
color = NA
),
panel.grid.minor = element_line(
color = NA
),
plot.background = element_rect(
fill = NA, color = NA
),
legend.background = element_rect(
fill = "white", color = NA
),
panel.border = element_rect(
fill = NA, color = NA
),
plot.margin = unit(
c(
t = 0, r = 0,
b = 0, l = 0
), "lines"
)
)
}
warna <- hcl.colors(
n = 5,
palette = "Temps",
rev = TRUE
)
# Curah hujan bulanan (mm per bulan)
batas_int_bulanan <- c(0, 100, 300, 400, ceiling(max(df_mswep$curah_hujan)))
batas_int_bulanan_baku <- (batas_int_bulanan - min(batas_int_bulanan)) / (max(batas_int_bulanan) - min(batas_int_bulanan))
# Curah hujan tahunan (mm per tahun)
batas_int_tahunan <- c(0, 1500, 3000, 4500, ceiling(max(df_mswep_tahunan$curah_hujan_tahunan)))
batas_int_tahunan_baku <- (batas_int_tahunan - min(batas_int_tahunan)) / (max(batas_int_tahunan) - min(batas_int_tahunan))Sebagai catatan, interval pada baris kode di atas memiliki batas-batas batas_int_tahunan dan batas_int_bulanan. Batas-batas interval bulanan maupun tahunan tersebut dibuat dengan berdasarkan kategori curah hujan BMKG.
Sekarang, kita terapkan tema, interval, dan warna yang sudah kita siapkan untuk membuat peta. Hasilnya disajikan pada Gambar 4.
df_mswep_19_22 |>
ggplot() +
geom_raster(
aes(
x = x,
y = y,
fill = curah_hujan_tahunan
)
) +
geom_contour(
aes(
x = x,
y = y,
z = curah_hujan_tahunan
),
color = "white",
linewidth = .25
) +
geom_sf(
data = sf_negara,
fill = "transparent",
color = "grey10",
linewidth = .5
) +
scale_fill_gradientn(
name = "Curah hujan\n(mm/tahun)",
colors = warna,
values = batas_int_tahunan_baku,
breaks = c(0, 1500, 3000, 4500),
labels = c(0, 1500, 3000, 4500)
) +
facet_wrap(vars(tahun), ncol = 2) +
tema_peta()
Berdasarkan peta pada Gambar 4, kita dapat melihat bahwa tahun 2019 curah hujannya relatif lebih rendah dibandingkan dengan tahun 2020 – 2022. Tak hanya itu, kita juga dapat membandingkan curah hujan antarwilayah di Indonesia. Tampak bahwa wilayah pesisir barat Pulau Sumatera, Kalimantan bagian utara, dan Papua bagian timur memiliki curah hujan yang relatif lebih tinggi.
Sekarang kita akan membuat sebuah peta curah hujan di Indonesia pada bulan dan tahun tertentu. Apakah kamu punya usulan bulan dan tahun berapa yang perlu kita gambar? Kita pilih bulan dan tahun yang terakhir saja, yaitu Januari 2023. Perhatikan Gambar 5!
peta_2023_01 <- df_mswep |>
filter(
tanggal_bawah == as.Date("2023-01-01")
) |>
ggplot() +
geom_raster(
aes(
x = x,
y = y,
fill = curah_hujan
)
) +
geom_contour(
aes(
x = x,
y = y,
z = curah_hujan
),
breaks = c(0, 100, 300, 400),
color = "white",
linewidth = .25
) +
geom_sf(
data = sf_negara,
fill = "transparent",
color = "grey10",
linewidth = .5
) +
scale_fill_gradientn(
name = "Curah hujan\n(mm/bulan)",
colors = warna,
values = batas_int_bulanan_baku,
breaks = c(0, 100, 300, 400),
labels = c(0, 100, 300, 400)
) +
tema_peta()
print(peta_2023_01)
Berdasarkan Gambar 5, kita dapat melihat bahwa sebagian besar wilayah Indonesia memiliki curah hujan yang cukup tinggi. Beberapa wilayah saja yang curah hujannya relatif rendah, yaitu sebagian wilayah Sulawesi, bagian selatan Kepulauan Maluku, dan wilayah Nusa Tenggara. Untuk lebih jelasnya, kita dapat melihat distribusi curah hujannya pada Januari 2023, seperti yang ditunjukkan pada Gambar 6!
df_mswep |>
filter(
tanggal_bawah == as.Date("2023-01-01")
) |>
ggplot(aes(x = curah_hujan)) +
geom_density() +
theme_minimal() +
theme(
axis.title.y = element_blank(),
axis.text.y = element_blank()
) +
labs(
x = "Curah hujan (mm/bulan)"
)
Visualisasi Data dengan Peta 3D
Peta yang kita hasilkan pada Bagian 3 merupakan peta 2D meskipun terbantu oleh kontur sehingga tampak seperti 3D. Apakah kita dapat membuat peta 3D? Ya, kita dapat melakukannya dengan menggunakan paket {rayshader}.
Untuk membuat peta 3D dari peta_2023_01 yang telah kita buat sebelumnya, kita gunakan baris kode berikut.
plot_gg(
ggobj = peta_2023_01,
width = 7,
height = 7,
scale = 250,
solid = FALSE,
shadow = TRUE,
shadowcolor = "white",
shadowwidth = 0,
shadow_intensity = 1,
zoom = .7,
phi = 85,
theta = 0
)Fungsi plot_gg() tersebut gunanya untuk menjadikan objek peta_2023_01 menjadi 3D dengan memetakan warnanya sebagai ketinggian dalam gambar 3D tersebut. Setelah baris perintah di atas dijalankan, dalam beberapa saat kita akan mendapatkan jendela baru yang muncul. Jendela tersebut seperti yang ditunjukkan pada Gambar 7.
Ketika jendela kamera tersebut masih muncul, kita dapat mengatur posisi dan propertinya dengan menggunakan fungsi render_camera(). Kita dapat mengaturnya seperti pada baris kode berikut.
Argumen theta = 30 menyatakan bahwa sudut rotasinya adalah 30 derajat. Argumen phi = 75 menyatakan bahwa sudut azimutnya sebesar 75 derajat. Setelah itu, kita dapat merender peta 3D tersebut dengan menggunakan baris kode sebagai berikut.
u <- "https://dl.polyhaven.org/file/ph-assets/HDRIs/hdr/4k/air_museum_playground_4k.hdr"
fail_hdri <- basename(u)
download.file(
url = u,
destfile = fail_hdri,
mode = "wb"
)
render_highquality(
filename = "curah_hujan_indonesia.png",
preview = TRUE,
interactive = FALSE,
parallel = TRUE,
light = TRUE,
environment_light = fail_hdri,
intensity = .45,
rotate_env = 90,
width = 1280,
height = 720
)Hasil render tersebut ditunjukkan pada Gambar 8.
Catatan Akhir
Kita telah menganalisis curah hujan di Indonesia dari waktu ke waktu. Secara lebih spesifik, kita telah melihat bulan-bulan apa saja yang curah hujannya relatif lebih tinggi. Selain itu, kita juga telah melihat bagaimana peta, baik 2D maupun 3D, dapat kita gunakan untuk memvisualisasikan data curah hujan dari berbagai wilayah Indonesia. Dengan cara ini, kita dapat mengidentifikasi daerah-daerah mana yang curah hujannya lebih tinggi dan demikian juga sebaliknya.