Appendix 2.
R scripts for the short RP and RV estimation
## set working directory
## setwd("C:/Users/CHY/Desktop/바탕_임시폴더/Desktop_folder_arrangement_2024_0802")
library("lubridate")
library("pracma")
library("eva")
idata <- read.csv("DD_wave_data_2012_2023.csv") ## Data format: date-time, Hs(m)
ndata <- nrow(idata) ## before removing of the NA data
## Sample data : 2012 - 2023, 12 years data
NYEAR <- 12.0; nleap <- 3
ndays <- 366*nleap + 365*(NYEAR-nleap) ## 3 leap years
HS <- idata$HS; digit <- 0.1
HS <- HS + digit*runif(ndata,-0.5, 0.5) ## rounding jittering
HS[which(HS < 0.0)] <- 0.0
DT <- ymd_hm(idata$DT)
dt <- 1.0 ## Unit: hours
mday <- 24/dt ## No. of data in a day
ncdata <- ndays*(24./dt)
MR0 <- ndata*100/ncdata ## Missing ratio before removing the NA data
## Basic and essential input parameters
cday <- 30.0 ## De-correlation time (unit: days; return period sub-scale)
TRP <- 1.0 ## Target return period & tolerance
RP <- TRP
maxiter <- 20 ## max. iteration number of the bisection method
epsilon <- 0.01 ## tolerance level
##++++++++++++++++++++++++++++++++++++++++++++++++++++
plot(DT,HS)
midx <- which(is.na(HS) == TRUE | HS == 0)
if(length(midx) > 0) {
HS <- HS[-midx]
DT <- DT[-midx]
ndata <- ndata-length(midx)
}
MR1 <- ndata*100/ncdata ## Missing ratio after removing the NA data
c(MR0, MR1)
corrected_NYEAR <- MR1*NYEAR/100
## Recording period is corrected using the missing ratio
summary(HS)
## Systematic grid-searching method:
Hseq <- seq(0,max(HS), 0.1)
ncat <- length(Hseq)
fresult <- matrix(0,nrow=ncat, ncol=3)
##
for (ii in 1:ncat) {
ncut <- findpeaks(HS, minpeakheight=Hseq[ii], minpeakdistance = mday*cday)
fresult[ii,1:2] <- c(Hseq[ii], nrow(ncut))
fresult[ii,3] <- corrected_NYEAR/nrow(ncut)
}
plot(fresult[,1], fresult[,3], ylim=c(0,NYEAR/3), type="o",
xlab="Hs(m)", ylab="Return Period (years)", cex.lab=1.3)
TRP <- 1.0
abline(h=TRP, col="red", lwd=2)
## Estimation using interpolation...
tidx <- which(fresult[,3] > TRP)[1]
c(fresult[tidx-1, 1], fresult[tidx,1])
inc_value <- (TRP- fresult[tidx-1,3])/
((fresult[tidx,3]-fresult[tidx-1,3])/(fresult[tidx,1]-fresult[tidx-1,1]))
est_value <- fresult[tidx-1,1] + inc_value
vstr <- as.character(signif(est_value, 3))
lines(c(est_value,est_value), c(0, TRP), col="blue", lty=1, lwd=3)
text(est_value-0.1, TRP+0.3, vstr, cex=1.2)
abline(h=0)
grid(lty=3, col="blue")
## Return value estimation fotr the given target return period using bi-section method
## initial guess, VM, using the lower and upper limits
VU <- max(HS); VL <- min(HS)
VM <- (VL+VU)/2
fresult2 <- matrix(0, nrow=maxiter, ncol=4)
fresult2[1,] <- c(VL,VM,VU,0)
for (ii in 1:maxiter) {
ncut <- findpeaks(HS, minpeakheight=VM, minpeakdistance = mday*cday)
ERP <- corrected_NYEAR/nrow(ncut)
if(abs(ERP-TRP) < epsilon) break
if(abs(VU-VL) < 0.005) break
if(ERP >= TRP) {
fresult2[ii,] <- c(VL, VM, VU, ERP)
VU <- VM; VM <- (VM+VL)/2
} else {
fresult2[ii,] <- c(VL, VM, VU, ERP)
VL <- VM; VM <- (VM+VU)/2
}
}
fresult2[ii,] <- c(VL, VM, VU, ERP)
fresult2 <- fresult2[1:ii,]
plot(fresult[,1], fresult[,3], ylim=c(0,NYEAR/3), type="o",
xlab="Hs(m)", ylab="Return Period (years)", cex.lab=1.3)
abline(h=TRP, col="red", lwd=2)
abline(v=fresult2[,2], col=1:ii, lwd=2)
text(fresult2[,2]-0.1, fresult2[,4]+1, as.character(1:ii), cex=0.8)
est_value <- as.character(signif(fresult2[ii,2],3))
text(fresult2[,2]-0.1, fresult2[,4]+1, as.character(1:ii), cex=0.8)
text(fresult2[ii,2]-0.05, fresult2[ii,4]+0.5, est_value, col="red", cex=1.2)