R untuk Hidrologi (1): hydroTSM

Instalasi paket

Menginstal versi paket stabil terakhir dari CRAN:

install.packages("hydroTSM")

Sebagai alternatif, teman-teman bisa juga menginstal versi yang sedang dikembangkan dari Github:

if (!require(devtools)) install.packages("devtools")
library(devtools)
install_github("hzambran/hydroTSM")

Mengatur lingkungan kerja

Memuat pustaka hydroTSM, yangmana berisi data dan fungsi-fungsi yang akan digunakan dalam analisis ini.

library(hydroTSM)

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

## Loading required package: xts

Memuat data hujan harian dari stasiun pengamatan San Martino di Castrozza, Provinsi Trento, Italia. Data sejak 01 Jan 1921 hingga 31 Des 1990.

data(SanMartinoPPts)

Memilih hanya potongan data dengan durasi selama 6 tahun untuk analisis.

x <- window(SanMartinoPPts, start=as.Date("1985-01-01"))

Mengubah/menjumlah nilai hujan harian menjadi bulanan.

( m <- daily2monthly(x, FUN=sum) )

## 1985-01-01 1985-02-01 1985-03-01 1985-04-01 1985-05-01 1985-06-01 1985-07-01 
##      141.2        7.0      140.6       72.0      175.6      131.4       85.4 
## 1985-08-01 1985-09-01 1985-10-01 1985-11-01 1985-12-01 1986-01-01 1986-02-01 
##      159.4       27.2       58.4      101.8       54.8       75.8      131.6 
## 1986-03-01 1986-04-01 1986-05-01 1986-06-01 1986-07-01 1986-08-01 1986-09-01 
##       59.6      237.8      108.2      144.8       81.2      141.0       69.8 
## 1986-10-01 1986-11-01 1986-12-01 1987-01-01 1987-02-01 1987-03-01 1987-04-01 
##       38.2       44.4       20.4       46.8      111.0       45.6       98.4 
## 1987-05-01 1987-06-01 1987-07-01 1987-08-01 1987-09-01 1987-10-01 1987-11-01 
##      212.0      153.8      221.8      175.0       90.6      278.8      164.8 
## 1987-12-01 1988-01-01 1988-02-01 1988-03-01 1988-04-01 1988-05-01 1988-06-01 
##       29.8      118.0       49.8       22.4      100.6      187.4      193.0 
## 1988-07-01 1988-08-01 1988-09-01 1988-10-01 1988-11-01 1988-12-01 1989-01-01 
##      120.4      149.2       61.2      136.4       10.0       59.4        0.0 
## 1989-02-01 1989-03-01 1989-04-01 1989-05-01 1989-06-01 1989-07-01 1989-08-01 
##      152.6       46.2      365.4       77.4      241.6      302.8      114.4 
## 1989-09-01 1989-10-01 1989-11-01 1989-12-01 1990-01-01 1990-02-01 1990-03-01 
##       65.4       12.8      145.0      110.6       51.6       12.4       65.8 
## 1990-04-01 1990-05-01 1990-06-01 1990-07-01 1990-08-01 1990-09-01 1990-10-01 
##      127.0       74.4      175.0      143.8       90.8      106.0      153.0 
## 1990-11-01 1990-12-01 
##      326.6      106.0

Mengekstraksi nilai tanggal dari data harian x.

dates <- time(x)

Banyaknya tahun dalam x (dibutuhkan untuk perhitungan)

( nyears <- yip(from=start(x), to=end(x), out.type="nmbr" ) )

## [1] 6

Analisis data eksplorasi dasar (exploratory data analysis / EDA)

1. Ringkasan statistika

smry(x)

##               Index         x
## Min.     1985-01-01    0.0000
## 1st Qu.  1986-07-02    0.0000
## Median   1988-01-01    0.0000
## Mean     1988-01-01    3.7470
## 3rd Qu.  1989-07-01    2.6000
## Max.     1990-12-31  122.0000
## IQR            <NA>    2.6000
## sd             <NA>   10.0428
## cv             <NA>    2.6800
## Skewness       <NA>    5.3512
## Kurtosis       <NA>   39.1619
## NA's           <NA>    0.0000
## n              <NA> 2191.0000

Menggunakan fungsi hydroplot, yangmana (secara standar) akan memplot 9 grafik berbeda, yakni: 3 plot seri waktu, 3 boxplot dan 3 histogram yang meringkaas x. Pada contoh ini, hanya plot-plot harian dan bulanan yang akan dihasilkan, dan hanya data yang dimulai pada 01 Jan 1987 yang akan diplot.

hydroplot(x, var.type="Precipitation", main="at San Martino", 
          pfreq = "dm", from="1987-01-01")

2. Banyaknya hari dengan informasi hujan (termasuk 0, bukan NA) per tahun

dwi(x)

## 1985 1986 1987 1988 1989 1990 
##  365  365  365  366  365  365

3. Banyaknya hari dengan informasi hujan (termasuk 0, bukan NA) per bulan per tahun

dwi(x, out.unit="mpy")

##      Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 1985  31  28  31  30  31  30  31  31  30  31  30  31
## 1986  31  28  31  30  31  30  31  31  30  31  30  31
## 1987  31  28  31  30  31  30  31  31  30  31  30  31
## 1988  31  29  31  30  31  30  31  31  30  31  30  31
## 1989  31  28  31  30  31  30  31  31  30  31  30  31
## 1990  31  28  31  30  31  30  31  31  30  31  30  31

4. Memplot nilai hujan bulanan tiap tahun, berguna untuk mengidentifikasi bulan-bulan basah/kering.

# kelas data zoo (ts) harian menjadi bulanan
m <- daily2monthly(x, FUN=sum, na.rm=TRUE)
     
# membuat sebuah matrix dengan nilai bulanan per tahun dalam tiap kolom
M <- matrix(m, ncol=12, byrow=TRUE)
colnames(M) <- month.abb
rownames(M) <- unique(format(time(m), "%Y"))
     
# memplot nilai hujan bulanan
require(lattice)

## Loading required package: lattice

print(matrixplot(M, ColorRamp="Precipitation", 
           main="Monthly precipitation at San Martino st., [mm/month]"))

Analisis tahunan

Nilai-nilai hujan tahunan (annual)

daily2annual(x, FUN=sum, na.rm=TRUE)

## 1985-01-01 1986-01-01 1987-01-01 1988-01-01 1989-01-01 1990-01-01 
##     1154.8     1152.8     1628.4     1207.8     1634.2     1432.4

Rata-rata hujan tahunan

Cara yang mudah dipahami:

mean( daily2annual(x, FUN=sum, na.rm=TRUE) )

## [1] 1368.4

Cara lain (lebih berguna untuk analisis aliran, dimana FUN=mean):

Fungsi annualfunction menerapkan FUN dua kali pada x:

( i) pertama, pada seluruh elemen x termasuk ke tahun yang sama, untuk mendapatkan nilai tahunan yang sesuai, dan

2. kedua, pada seluruh nilai tahunan x yang sebelumnya diperoleh, untuk mendapatkan nilai tahunan tunggal.

annualfunction(x, FUN=sum, na.rm=TRUE) / nyears

##  value 
## 1368.4

Analisis bulanan

Nilai tengah (median) dari hujan bulanan pada station x. Tidak diperlukan, hanya untuk melihat nilai-nilainya melalui boxplot.

monthlyfunction(m, FUN=median, na.rm=TRUE)

##   Jan   Feb   Mar   Apr   May   Jun   Jul   Aug   Sep   Oct   Nov   Dec 
##  63.7  80.4  52.9 113.8 141.9 164.4 132.1 145.1  67.6  97.4 123.4  57.1

Vektor dengan nama bulan yang memiliki format singkatan tiga huruf.

cmonth <- format(time(m), "%b")

Membuat tipe data bulan menjadi factor yang memiliki bobot terurut.

months <- factor(cmonth, levels=unique(cmonth), ordered=TRUE)

Boxplot nilai hujan bulanan

boxplot( coredata(m) ~ months, col="lightblue", main="Monthly Precipitation", 
         ylab="Precipitation, [mm]", xlab="Month")

Analisis musiman

Rata-rata nilai hujan musiman.

seasonalfunction(x, FUN=sum, na.rm=TRUE) / nyears

##      DJF      MAM      JJA      SON 
## 213.1333 369.4000 470.8000 315.0667

Mengekstrak nilai musiman untuk tiap tahun.

( DJF <- dm2seasonal(x, season="DJF", FUN=sum) )

##  1985  1986  1987  1988  1989  1990 
## 148.2 262.2 178.2 197.6 212.0 174.6

( MAM <- dm2seasonal(m, season="MAM", FUN=sum) )

##  1985  1986  1987  1988  1989  1990 
## 388.2 405.6 356.0 310.4 489.0 267.2

( JJA <- dm2seasonal(m, season="JJA", FUN=sum) )

##  1985  1986  1987  1988  1989  1990 
## 376.2 367.0 550.6 462.6 658.8 409.6

( SON <- dm2seasonal(m, season="SON", FUN=sum) )

##  1985  1986  1987  1988  1989  1990 
## 187.4 152.4 534.2 207.6 223.2 585.6

Memplot evolusi waktu dari nilai hujan musiman.

hydroplot(x, pfreq="seasonal", FUN=sum, stype="default")

Beberapa indikasi ekstrim

Langkah yang umum/biasa untuk bagian analisis ini:

Memuat data hujan harian dari stasiun pengamatan San Martino di Castrozza, Provinsi Trento, Italia, dengan data sejak 01 Jan 1921 hingga 31 Des 1990.

data(SanMartinoPPts)

Memilih hanya potongan data dengan durasi selama 3 tahun untuk analisis.

x <- window(SanMartinoPPts, start=as.Date("1988-01-01"))

Memplot seri waktu terpilih:

hydroplot(x,  ptype="ts", pfreq="o", var.unit="mm")

( R10mm <- length( x[x>10] ) )

## [1] 127

wet.index <- which(x >= 1)

( PRwn95 <- quantile(x[wet.index], probs=0.95, na.rm=TRUE) )

##   95% 
## 39.75

(very.wet.index <- which(x >= PRwn95))

##  [1]   30   92  234  287  422  423  461  550  551  674  676  719  939  950  998
## [16] 1058 1061 1075

( R95p <- sum(x[very.wet.index]) )

## [1] 1196.4

x.5max <- rollapply(data=x, width=5, FUN=sum, fill=NA, partial= TRUE, 
                    align="center")

hydroplot(x.5max,  ptype="ts+boxplot", pfreq="o", var.unit="mm")

## [Note: pfreq='o' => ptype has been changed to 'ts']

(x.5max.annual <- daily2annual(x.5max, FUN=max, na.rm=TRUE))

## 1988-01-01 1989-01-01 1990-01-01 
##      113.2      170.8      237.2

Klimograf

Sejak versi 0.5-0, hydroTSM menyertakan sebuah fungsi untuk memplot klimograf, tidak hanya data hujan/presipitasi tetapi juga data suhu udara seperti:

# memuat seri waktu harian dari hujan, max dan min suhu
data(MaquehueTemuco)

# mengekstrak seri waktu individual dari hujan, max dan min suhu
pcp <- MaquehueTemuco[, 1]
tmx <- MaquehueTemuco[, 2]
tmn <- MaquehueTemuco[, 3]
     
# memplot klimograf
m <- climograph(pcp=pcp, tmx=tmx, tmn=tmn, na.rm=TRUE)

Detail piranti lunak

Tutorial ini dibuat dengan spesifikasi komputer seperti berikut (versi penerjemah):

## [1] "x86_64-pc-linux-gnu (64-bit)"

## [1] "R version 3.6.3 (2020-02-29)"

## [1] "hydroTSM 0.6-0"

Lampiran

Agar mempermudah penggunaan hydroTSM bagi teman-teman yang belum familiar dengan R, penulis menyediakan sekumpulan informasi sebagai panduan untuk teman-teman menggunakan R.