theta = .75, cross.val = "loo")
scoreComputation <- function(obs, series, index, custom_function = NULL, weights = NULL,
methods){
index.obs <- c()
for (i in c(1:length(index))) {
index.obs <- c(index.obs, valueIndex(obs, index.code = index[i])$Index$Data)
}
aux.obs <- setTimeResolution(obs, time_resolution = "DD")
aux.obs <- setGridDates.asPOSIXlt(aux.obs, tz = "UTC")
index.list <- list()
for (j in c(1:length(methods))) {
index.cal <- c()
for (i in c(1:length(index))) {
index.cal <- c(index.cal, valueIndex(grid = series[[j]], index.code = index[i])$Index$Data)
}
aux.obs$Dates$start <- as.POSIXct(aux.obs$Dates$start, tz = "CEST")
aux.obs$Dates$end <- as.POSIXct(aux.obs$Dates$end, tz = "CEST")
if (length(custom_function) > 0) {
for (i in c(1:length(custom_function))) {
index.obs <- c(index.obs, custom_function[[i]](aux.obs))
index.cal <- c(index.cal, custom_function[[i]](series[[j]]))
}
}
index.list[[j]] <- index.cal
}
names(index.list) <- methods
normalization <- function(measure){
measure.norm <- c()
for (i in c(1:length(measure))) {
measure.norm <- c(measure.norm, 1-((measure[i]-min(measure))/(max(measure)-min(measure))))
}
return(measure.norm)
}
measures <- list()
for (i in c(1:length(index.cal))) {
aux <- c()
for (j in c(1:length(methods))) {
aux <- c(aux, abs(index.list[[j]][i]-index.obs[i]))
}
measures[[length(measures)+1]] <- aux
}
aux <- NULL
norm.vector <- list()
for (i in c(1:length(measures))) {
norm.vector[[length(norm.vector)+1]] <- normalization(measures[[i]])
}
scores <- c()
for (j in c(1:(length(methods)))) {
score <- c()
for (i in c(1:length(measures))) {
score <- c(score,norm.vector[[i]][j])
}
if (length(weights > 0)) {
score <- weighted.mean(score, w = weights)
}else{
score <- mean(score)
}
scores <- c(scores, score)
}
names(scores) <- methods
return(scores)
}
score <- scoreComputation(obs = obs, series = list(scaling, eqm, pqm, gpqm95, gpqm75, cal),
index = c("Skewness","Mean","SDII","R10","R10p","R20","R20p","P98Wet","P98WetAmount"),
custom_function = list(MaxReturnValue), weights = NULL,
methods = c("scaling", "eqm", "pqm", "gpqm95", "gpqm75", "adaptive"))
print(score)
sessionInfo()
help(subsetGrid)
obs$Data
obs$Metadata
help("intersectGrid")
help(biasCorrection)
.libPaths("/home/meteo/miniconda3/envs/oscarmeteo/lib/R/library")
install.packages('rlang')
install.packages("rlang")
library(devtools)
load_all("/home/meteo/transformeR")
load_all("/home/meteo/Escritorio/downscaleR")
load_all("/home/meteo/climate4R.value")
library(loadeR)
temp
library(knitr)
getwd()
knit2pdf('/home/meteo/Escritorio/adaptive_notebook.Rmd')
temp
unzip(temp)
temp
ex <- unzip(temp)
ex
rm(ex)
loadStationData(temp, var = 'pr')
showConnections(all = FALSE)
showConnections()
help("download.file")
help("tempdir")
tempdir()
tempdir()
temp <- tempdir()
download.file("https://zenodo.org/record/7014397/files/South_Pacific_precipitation.zip",temp)
temp <- tempfile(pattern = "file", tmpdir = tempdir(), fileext = "")
download.file("https://zenodo.org/record/7014397/files/South_Pacific_precipitation.zip",temp)
loadStationData(dataset = temp, var = "pr")
knitr::opts_chunk$set(echo = TRUE,
fig.align = "center",
fig.width = 7.5,
fig.asp = 1,
message = FALSE,
tidy = FALSE,
cache = TRUE)
temp <- tempfile(fileext = ".zip")
download.file("https://zenodo.org/record/7014397/files/South_Pacific_precipitation.zip",temp)
loadStationData(dataset = temp, var = "pr")
loadStationData(dataset = "https://zenodo.org/record/7014397/files/South_Pacific_precipitation.zip", var = "pr")
temp <- tempfile()
download.file("https://zenodo.org/record/7014397/files/South_Pacific_precipitation.zip",temp)
obs <- loadStationData(dataset = temp, var = "pr")
loadStationData(dataset = temp, var = "pr")
temp <- tempfile(fileext = ".zip")
download.file("https://zenodo.org/record/7014397/files/South_Pacific_precipitation.zip",temp)
obs <- loadStationData(dataset = temp, var = "pr")
trmm <- loadStationData(dataset = temp, var = "pp_trmm")
wt <- loadStationData(dataset = temp, var = "wt")
temp <- tempfile(fileext = ".zip")
download.file("https://zenodo.org/record/7014397/files/South_Pacific_precipitation.zip",temp, mode = "wb")
obs <- loadStationData(dataset = temp, var = "pr")
trmm <- loadStationData(dataset = temp, var = "pp_trmm")
wt <- loadStationData(dataset = temp, var = "wt")
temp <- "/home/meteo/Descargas/South_Pacific_precipitation.zip"
loadStationData(dataset = temp, var = "pr")
temp
dataInventory(temp)
loadStationData(dataset = temp, var = "pr")
help("loadStationData")
.libPaths("/home/meteo/miniconda3/envs/oscarmeteo/lib/R/library")
temp
knitr::opts_chunk$set(echo = TRUE,
fig.align = "center",
fig.width = 7.5,
fig.asp = 1,
message = FALSE,
tidy = FALSE,
cache = TRUE)
temp <- "/home/meteo/Descargas/South_Pacific_precipitation.zip"
closeAllConnections()
rm(list=ls())
unlink("Escritorio/adaptive_notebook_cache", recursive = TRUE)
temp <- "/home/meteo/Descargas/South_Pacific_precipitation.zip"
loadStationData(dataset = temp, var = "pr")
install.packages('rlang')
.libPaths("/home/meteo/miniconda3/envs/oscarmeteo/lib/R/library")
install.packages('rlang')
install.packages("rlang")
library(devtools)
load_all("/home/meteo/transformeR")
load_all("/home/meteo/Escritorio/downscaleR")
load_all("/home/meteo/climate4R.value")
library(loadeR)
knitr::opts_chunk$set(echo = TRUE,
fig.align = "center",
fig.width = 7.5,
fig.asp = 1,
message = FALSE,
tidy = FALSE,
cache = TRUE)
temp <- "/home/meteo/Descargas/South_Pacific_precipitation.zip"
loadGridData(temp, var = "pr")
loadStationData(temp, var = "pr")
temp
unzip(temp)
obs <- loadStationData(dataset = temp, var = "pr")
trmm <- loadStationData(dataset = temp, var = "pp_trmm")
wt <- loadStationData(dataset = temp, var = "wt")
getwd()
setwd(/home/meteo/Escritorio)
setwd("/home/meteo/Escritorio")
getwd()
save("obs.RData",obs)
save(obs, "obs.RData")
save(file = "obs.RData",obs)
save(file = "trmm.RData",obs)
save(file = "trmm.RData",trmm)
save(file = "wt.RData",wt)
library(VALUE)
library(loadeR)
library(transformeR)
library(geosphere)
library(scales)
library(abind)
library(loadeR)
library(transformeR)
library(downscaleR)
library(magrittr)
library(downscaleR.keras)
library(tfprobability)
setwd("/home/meteo/Escritorio/cluster/WORK/ARTICULOS/2022_FWI_deep/data")
y <- loadStationData(dataset = "predictands",
var = "fwi13",
years = 1985:2011,
tz = "UTC")
#y <- filterNA(y)
preds <- list.files("predictions/")[c(1,5,4,2,3)]
colors <- c("gold", "red", "blue", "seagreen", "grey30")
region.list <- list()
for (ind_region in unique(y$Metadata$climatic_region_code)) {
print(paste0("#####Computing Boxplots for ",ind_region, " code#####"))
st.codes <- y$Metadata$station_id[which(y$Metadata$climatic_region_code == ind_region)]
y.aux <- subsetStation(y, station.id = st.codes)
predictions.list <- list()
for (ind_preds in c(1:length(preds))) {
x <- loadStationData(dataset = "predictions",
var = gsub("\\.txt$", "", preds[ind_preds]),
years = 1985:2011,
tz = "UTC")
x <- subsetStation(x, station.id = st.codes)
rmse.values <- c()
for (loc in c(1:length(st.codes))) {
x.st <- subsetStation(x, station.id = st.codes[loc])
y.st <- subsetStation(y, station.id = st.codes[loc])
q_x <- quantile(as.vector(x.st$Data), na.rm = TRUE, probs = .90)
q_y <- quantile(as.vector(y.st$Data), na.rm = TRUE, probs = .90)
x.st <- x.st$Data[which(x.st$Data>= q_x)]
y.st <- y.st$Data[which(y.st$Data>= q_y)]
rmse <- sqrt(mean((x.st - y.st)^2, na.rm = T))
rmse.values <- c(rmse.values, rmse)
}
predictions.list[[gsub("\\.txt$", "", preds[ind_preds])]] <- rmse.values
}
region.list[[ind_region]] <- predictions.list
}
stripchart(region.list$COASM,
method = "stack",
xlab = "",
ylab = "RMSE FWI90",
main = "",
col = colors,
pch = 19,
vertical = T,
ylim = range(c(region.list$ATL, region.list$COASM, region.list$CONTM)))
grid()
points_CONT <- list(
FWI_analogs = region.list$CONTM$`FWI-analogs`,
FWI_GLM = region.list$CONTM$`FWI-GLM`,
FWI_CNN_MS = region.list$CONTM$`FWI-CNN-MS`,
FWI_CNN_MS_G = region.list$CONTM$`FWI-CNN-MS-G`,
FWI_CNN_MS_MG = region.list$CONTM$`FWI-CNN-MS-MG`
)
points_ATL <- list(
FWI_analogs = region.list$ATL$`FWI-analogs`,
FWI_GLM = region.list$ATL$`FWI-GLM`,
FWI_CNN_MS = region.list$ATL$`FWI-CNN-MS`,
FWI_CNN_MS_G = region.list$ATL$`FWI-CNN-MS-G`,
FWI_CNN_MS_MG = region.list$ATL$`FWI-CNN-MS-MG`
)
for (i in 1:length(points_CONT)) {
points(points_CONT[[i]],
col = colors[i],
pch = 6,
vertical = TRUE)
}
for (i in 1:length(points_ATL)) {
points(points_ATL[[i]],
y = "FWI-analogs",
col = colors[i],
pch = 4,
vertical = TRUE)
}
#Add legend
legend("topleft", legend = c("ATL", "COASTM", "CONTM"),  pch = c(19, 6, 4))
par(mfrow = c(1,1))
stripchart(region.list$CONTM,
method = "stack",
xlab = "",
ylab = "RMSE FWI90",
main = "",
col = colors,
pch = 19,
vertical = T,
ylim = range(c(region.list$ATL, region.list$COASM, region.list$CONTM)))
stripchart(region.list$COASM,
method = "stack",
xlab = "",
ylab = "RMSE FWI90",
main = "Coastal Mediterranean",
col = colors,
pch = 6,
vertical = T,
add = T)
stripchart(region.list$ATL,
method = "stack",
xlab = "",
ylab = "RMSE FWI90",
main = "Atlantic",
col = colors,
pch = 4,
vertical = T,
add = T)
grid()
gc()
setwd("/home/meteo/Escritorio/cluster/WORK/ARTICULOS/2022_FWI_deep/data")
y <- loadStationData(dataset = "predictands",
var = "fwi13",
years = 1985:2011,
tz = "UTC")
#y <- filterNA(y)
preds <- list.files("predictions/")[c(1,5,4,2,3)]
colors <- c("gold", "red", "blue", "seagreen", "grey30")
region.list <- list()
for (ind_region in unique(y$Metadata$climatic_region_code)) {
print(paste0("#####Computing Boxplots for ",ind_region, " code#####"))
st.codes <- y$Metadata$station_id[which(y$Metadata$climatic_region_code == ind_region)]
y.aux <- subsetStation(y, station.id = st.codes)
predictions.list <- list()
for (ind_preds in c(1:length(preds))) {
x <- loadStationData(dataset = "predictions",
var = gsub("\\.txt$", "", preds[ind_preds]),
years = 1985:2011,
tz = "UTC")
x <- subsetStation(x, station.id = st.codes)
rmse.values <- c()
for (loc in c(1:length(st.codes))) {
x.st <- subsetStation(x, station.id = st.codes[loc])
y.st <- subsetStation(y, station.id = st.codes[loc])
# q_x <- quantile(as.vector(x.st$Data), na.rm = TRUE, probs = .90)
# q_y <- quantile(as.vector(y.st$Data), na.rm = TRUE, probs = .90)
#
# x.st <- x.st$Data[which(x.st$Data>= q_x)]
# y.st <- y.st$Data[which(y.st$Data>= q_y)]
rmse <- sqrt(mean((x.st - y.st)^2, na.rm = T))
rmse.values <- c(rmse.values, rmse)
}
predictions.list[[gsub("\\.txt$", "", preds[ind_preds])]] <- rmse.values
}
region.list[[ind_region]] <- predictions.list
}
gc()
setwd("/home/meteo/Escritorio/cluster/WORK/ARTICULOS/2022_FWI_deep/data")
y <- loadStationData(dataset = "predictands",
var = "fwi13",
years = 1985:2011,
tz = "UTC")
#y <- filterNA(y)
preds <- list.files("predictions/")[c(1,5,4,2,3)]
colors <- c("gold", "red", "blue", "seagreen", "grey30")
region.list <- list()
for (ind_region in unique(y$Metadata$climatic_region_code)) {
print(paste0("#####Computing Boxplots for ",ind_region, " code#####"))
st.codes <- y$Metadata$station_id[which(y$Metadata$climatic_region_code == ind_region)]
y.aux <- subsetStation(y, station.id = st.codes)
predictions.list <- list()
for (ind_preds in c(1:length(preds))) {
x <- loadStationData(dataset = "predictions",
var = gsub("\\.txt$", "", preds[ind_preds]),
years = 1985:2011,
tz = "UTC")
x <- subsetStation(x, station.id = st.codes)
rmse.values <- c()
for (loc in c(1:length(st.codes))) {
x.st <- subsetStation(x, station.id = st.codes[loc])
y.st <- subsetStation(y, station.id = st.codes[loc])
# q_x <- quantile(as.vector(x.st$Data), na.rm = TRUE, probs = .90)
# q_y <- quantile(as.vector(y.st$Data), na.rm = TRUE, probs = .90)
#
# x.st <- x.st$Data[which(x.st$Data>= q_x)]
# y.st <- y.st$Data[which(y.st$Data>= q_y)]
rmse <- sqrt(mean((x.st - y.st)^2, na.rm = T))
rmse.values <- c(rmse.values, rmse)
}
predictions.list[[gsub("\\.txt$", "", preds[ind_preds])]] <- rmse.values
}
region.list[[ind_region]] <- predictions.list
}
gc()
y <- loadStationData(dataset = "predictands",
var = "fwi13",
years = 1985:2011,
tz = "UTC")
#y <- filterNA(y)
preds <- list.files("predictions/")[c(1,5,4,2,3)]
colors <- c("gold", "red", "blue", "seagreen", "grey30")
region.list <- list()
for (ind_region in unique(y$Metadata$climatic_region_code)) {
print(paste0("#####Computing Boxplots for ",ind_region, " code#####"))
st.codes <- y$Metadata$station_id[which(y$Metadata$climatic_region_code == ind_region)]
y.aux <- subsetStation(y, station.id = st.codes)
predictions.list <- list()
for (ind_preds in c(1:length(preds))) {
x <- loadStationData(dataset = "predictions",
var = gsub("\\.txt$", "", preds[ind_preds]),
years = 1985:2011,
tz = "UTC")
x <- subsetStation(x, station.id = st.codes)
rmse.values <- c()
for (loc in c(1:length(st.codes))) {
x.st <- subsetStation(x, station.id = st.codes[loc])
y.st <- subsetStation(y, station.id = st.codes[loc])
######FWI90#######
# q_x <- quantile(as.vector(x.st$Data), na.rm = TRUE, probs = .90)
# q_y <- quantile(as.vector(y.st$Data), na.rm = TRUE, probs = .90)
#
# x.st <- x.st$Data[which(x.st$Data>= q_x)]
# y.st <- y.st$Data[which(y.st$Data>= q_y)]
#######FWI#######
x.st <- x.st$Data
y.st <- y.st$Data
rmse <- sqrt(mean((x.st - y.st)^2, na.rm = T))
rmse.values <- c(rmse.values, rmse)
}
predictions.list[[gsub("\\.txt$", "", preds[ind_preds])]] <- rmse.values
}
region.list[[ind_region]] <- predictions.list
}
par(mfrow = c(1,1))
stripchart(region.list$CONTM,
method = "stack",
xlab = "",
ylab = "RMSE FWI",
main = "",
col = colors,
pch = 19,
vertical = T,
ylim = range(c(region.list$ATL, region.list$COASM, region.list$CONTM)))
stripchart(region.list$COASM,
method = "stack",
xlab = "",
ylab = "RMSE FWI",
main = "Coastal Mediterranean",
col = colors,
pch = 6,
vertical = T,
add = T)
stripchart(region.list$ATL,
method = "stack",
xlab = "",
ylab = "RMSE FWI",
main = "Atlantic",
col = colors,
pch = 4,
vertical = T,
add = T)
grid()
legend("topleft", legend = c("ATL", "COASTM", "CONTM"),  pch = c(19, 6, 4))
gc()
setwd("/home/meteo/Escritorio/cluster/WORK/ARTICULOS/2022_FWI_deep/data")
y <- loadStationData(dataset = "predictands",
var = "fwi13",
years = 1985:2011,
tz = "UTC")
#y <- filterNA(y)
preds <- list.files("predictions/")[c(1,5,4,2,3)]
colors <- c("gold", "red", "blue", "seagreen", "grey30")
region.list <- list()
for (ind_region in unique(y$Metadata$climatic_region_code)) {
print(paste0("#####Computing Boxplots for ",ind_region, " code#####"))
st.codes <- y$Metadata$station_id[which(y$Metadata$climatic_region_code == ind_region)]
y.aux <- subsetStation(y, station.id = st.codes)
predictions.list <- list()
for (ind_preds in c(1:length(preds))) {
x <- loadStationData(dataset = "predictions",
var = gsub("\\.txt$", "", preds[ind_preds]),
years = 1985:2011,
tz = "UTC")
x <- subsetStation(x, station.id = st.codes)
rmse.values <- c()
for (loc in c(1:length(st.codes))) {
x.st <- subsetStation(x, station.id = st.codes[loc])
y.st <- subsetStation(y, station.id = st.codes[loc])
######FWI90#######
q_x <- quantile(as.vector(x.st$Data), na.rm = TRUE, probs = .90)
q_y <- quantile(as.vector(y.st$Data), na.rm = TRUE, probs = .90)
x.st <- x.st$Data[which(x.st$Data>= q_x)]
y.st <- y.st$Data[which(y.st$Data>= q_y)]
#######FWI#######
# x.st <- x.st$Data
# y.st <- y.st$Data
rmse <- sqrt(mean((x.st - y.st)^2, na.rm = T))
rmse.values <- c(rmse.values, rmse)
}
predictions.list[[gsub("\\.txt$", "", preds[ind_preds])]] <- rmse.values
}
region.list[[ind_region]] <- predictions.list
}
par(mfrow = c(1,1))
stripchart(region.list$CONTM,
method = "stack",
xlab = "",
ylab = "RMSE FWI90",
main = "",
col = colors,
pch = 19,
vertical = T,
ylim = range(c(region.list$ATL, region.list$COASM, region.list$CONTM)))
stripchart(region.list$COASM,
method = "stack",
xlab = "",
ylab = "RMSE FWI90",
main = "Coastal Mediterranean",
col = colors,
pch = 6,
vertical = T,
add = T)
stripchart(region.list$ATL,
method = "stack",
xlab = "",
ylab = "RMSE FWI90",
main = "Atlantic",
col = colors,
pch = 4,
vertical = T,
add = T)
grid()
legend("topleft", legend = c("ATL", "COASTM", "CONTM"),  pch = c(19, 6, 4))
