.Rhistory

fill = variable
), alpha = 0.5) +
geom_line(aes(y = value, color = variable)) +
geom_line(data = observations_data, aes(y = value, color = variable)) +
labs(y = "°C", x = "") +
scale_x_datetime(
limits = xlimits,
expand = c(0, 0),
date_breaks = "1 days",
date_labels = "%b %d"
) +
theme_light() + guides(fill = "none") +
theme(
plot.title = element_text(hjust = 0.5, size = 18),
legend.text = element_text(size = 10),
axis.text.x = element_text(size = 10),
axis.title.x = element_text(size = 10),
axis.text.y = element_text(size = 10),
axis.title.y = element_text(size = 10),
legend.box.background = element_rect(fill = 'transparent')
) +
scale_alpha(guide = 'none')
wrf <- readRDS(paste0(base_dir, "WRF_bothregions.RDS")) %>%
dplyr::filter(ID == ejemplo) %>%
dplyr::select(c(time, reftime, lead, WRF_T, ID))
all_wrf <- dplyr::left_join(all, wrf, by = c("reftime", "time", "lead", "ID"))
prediction <- tibble();
for (l in unique_leads){
train <- all_wrf %>% filter(reftime < as.POSIXct(train_date, tz = 'UTC') & lead == l)
val <- all_wrf %>% filter(reftime >= as.POSIXct(train_date, tz = 'UTC') & lead == l)
if (l <= 72) { #
emos_formula   <- as.character("OBS ~ IFS_T_mea + WRF_T | IFS_T_std ")
} else {
emos_formula   <- as.character("OBS ~ IFS_T_mea | IFS_T_std ")
}
emosfit <- crch(formula = emos_formula,
dat     = train,
dist    = "gaussian",
type    = "crps")
predres <- predict(object = emosfit,
newdata  = val,
type = "parameter")
prediction <- predres %>%
bind_cols(lead=val$lead,reftime=val$reftime,time=val$time) %>%
rename(T_mea_p = location, T_std_p = scale) %>%
bind_rows(prediction)
}
all_emos_multi <- all_emos %>%
filter(reftime >= as.POSIXct(train_date, tz = 'UTC')) %>%
right_join(prediction) %>%
arrange(reftime,lead)
all_emos_multi %>%
dplyr::mutate(
CRPS = SpecsVerification::GaussCrps(IFS_T_mea, IFS_T_std, OBS),
CRPS_p = SpecsVerification::GaussCrps(IFS_T_mea_p, IFS_T_std_p, OBS),
CRPS_pm = SpecsVerification::GaussCrps(T_mea_p, T_std_p, OBS),
) %>%
dplyr::group_by(lead) %>%
dplyr::summarize(DMO = mean(CRPS,na.rm=T), EMOS_IFS = mean(CRPS_p,na.rm=T), EMOS_MULTI = mean(CRPS_pm,na.rm=T)) %>%
tidyr::gather(variable, value, -lead) %>%
ggplot(aes(x=lead, y=value, colour=variable)) +
ggtitle("EMOS postprocesamiento multimodelo") +
scale_color_manual(name = '',values = c("red","skyblue", "darkblue")) +
theme_bw() +
labs(y = "CRPS", x= "lead time [h]", colour= "Forecast") +
theme(plot.title = element_text(hjust = 0.5)) +
geom_line()
gemos_formula <- as.character("OBS ~ IFS_T_mea + ELEV + TPI5 + TPI20 + TPI100 | IFS_T_std")
gemos_predictors <- readRDS(paste0(base_dir, "preproc_train.RDS")) %>%
as_tibble()
gemos <- ifs %>%
dplyr::left_join(obs, by = c("time", "ID")) %>%
dplyr::select(c(IFS_T_std, IFS_T_mea, time, reftime, lead, ID, OBS)) %>%
dplyr::left_join(gemos_predictors) %>%
dplyr::filter(ID %in% gemos_predictors$ID) %>%
arrange(ID, reftime,lead)
prediction <- tibble();
for (l in unique_leads) {
print(paste0("Training gEMOS IFS data for all stations and lead time ", l))
train <- gemos %>%
filter(reftime < as.POSIXct(train_date, tz = 'UTC') & lead == l)
gemosfit <- crch(
formula = stats::as.formula(gemos_formula),
dat     = train,
dist    = "gaussian",
type    = "crps"
)
val <- gemos %>%
filter(reftime > as.POSIXct(train_date, tz = 'UTC') & lead == l)
predres <- predict(object = gemosfit,
newdata  = val,
type = "parameter")
prediction <- predres %>%
bind_cols(lead=val$lead, reftime=val$reftime, time=val$time, ID = val$ID) %>%
rename(IFS_T_mea_p = location, IFS_T_std_p = scale) %>%
bind_rows(prediction)
}
all_gemos <- gemos %>%
dplyr::filter(reftime >= as.POSIXct(train_date, tz = 'UTC')) %>%
right_join(prediction) %>%
dplyr::select(IFS_T_mea, IFS_T_std, IFS_T_mea_p, IFS_T_std_p, reftime, lead, ID, OBS) %>%
arrange(ID, reftime,lead)
all_gemos %>%
dplyr::mutate(
CRPS = SpecsVerification::GaussCrps(IFS_T_mea, IFS_T_std, OBS),
CRPS_p = SpecsVerification::GaussCrps(IFS_T_mea_p, IFS_T_std_p, OBS),
) %>%
dplyr::group_by(lead) %>%
dplyr::summarize(DMO = mean(CRPS,na.rm=T), gEMOS = mean(CRPS_p,na.rm=T)) %>%
tidyr::gather(variable, value, -lead) %>%
ggplot(aes(x=lead, y=value, colour=variable)) +
ggtitle("EMOS postprocesamiento global") +
scale_color_manual(name = '',values = c("red","skyblue", "darkblue")) +
theme_bw() +
labs(y = "CRPS", x= "lead time [h]", colour= "Forecast") +
theme(plot.title = element_text(hjust = 0.5)) +
geom_line()
# Disponible para la region Cuyo:
# 87418 87436 87311 87416 87506 87509 87420
ejemplo2 <- 87047
#87418
obs_station <- obs %>%
dplyr::filter(ID == ejemplo)
ifs_station <- ifs %>%
dplyr::filter(ID == ejemplo)%>%
dplyr::select(-c(Lat, Lon))
all <- dplyr::right_join(ifs_station, obs_station, by = c("time", "ID"))
prediction <- tibble();
for (l in unique_leads) {
print(paste0("Training EMOS and gEMOS IFS data lead time ", l))
train_gemos <- gemos %>%
dplyr::filter(reftime < as.POSIXct(train_date, tz = 'UTC') & lead == l)
train_emos <- gemos %>%
dplyr::filter(reftime < as.POSIXct(train_date, tz = 'UTC') & lead == l & ID == ejemplo2)
val <- gemos %>%
dplyr::filter(reftime > as.POSIXct(train_date, tz = 'UTC') & lead == l & ID == ejemplo2)
emosfit <- crch(
formula = stats::as.formula(emos_formula),
dat     = train_emos,
dist    = "gaussian",
type    = "crps"
)
predres_emos <- predict(object = emosfit,
newdata  = val,
type = "parameter")
gemosfit <- crch(
formula = stats::as.formula(gemos_formula),
dat     = train_gemos,
dist    = "gaussian",
type    = "crps"
)
predres_gemos <- predict(object = gemosfit,
newdata  = val,
type = "parameter")
prediction <- predres_emos %>%
dplyr::rename(IFS_T_mea_emos = location, IFS_T_std_emos = scale) %>%
bind_cols(IFS_T_mea=val$IFS_T_mea,
IFS_T_std=val$IFS_T_std,
lead=val$lead,
reftime=val$reftime,
time=val$time,
OBS=val$OBS
) %>%
bind_cols(predres_gemos) %>%
dplyr::rename(IFS_T_mea_gemos = location, IFS_T_std_gemos = scale) %>%
bind_rows(prediction)
}
emos_formula   <- as.character("OBS ~ IFS_T_mea  | IFS_T_std ")
train_date <- "2023-06-30" #elegir período de training
prediction <- tibble();
unique_leads <- unique(all$lead) %>%
na.omit()
for (l in unique_leads){
print(paste0("Training EMOS IFS data for lead time ", l))
train <- all %>% filter(reftime < as.POSIXct(train_date, tz = 'UTC') & lead == l)
val <- all %>% filter(reftime >= as.POSIXct(train_date, tz = 'UTC') & lead == l)
emosfit <- crch(formula = emos_formula,
dat     = train,
dist    = "gaussian",
type    = "crps")
predres <- predict(object = emosfit,
newdata  = val,
type = "parameter")
prediction <- predres %>%
bind_cols(lead=val$lead,reftime=val$reftime,time=val$time) %>%
rename(IFS_T_mea_p = location, IFS_T_std_p = scale) %>%
bind_rows(prediction)
}
all_emos <- all %>%
filter(reftime >= as.POSIXct(train_date, tz = 'UTC')) %>%
right_join(prediction) %>% arrange(reftime,lead)
print(all_emos)
# Este gráfico muestra la comparación de los CRPS del modelo IFS (salida directa del modelo DMO) y espués de haber aplicado el método EMOS.
all_emos %>%
dplyr::mutate(
CRPS = SpecsVerification::GaussCrps(IFS_T_mea, IFS_T_std, OBS),
CRPS_p = SpecsVerification::GaussCrps(IFS_T_mea_p, IFS_T_std_p, OBS),
) %>%
dplyr::group_by(lead) %>%
dplyr::summarize(DMO = mean(CRPS,na.rm=T), EMOS = mean(CRPS_p,na.rm=T)) %>%
tidyr::gather(variable, value, -lead) %>%
ggplot(aes(x=lead, y=value, colour=variable)) + ggtitle("EMOS postprocesamiento") +
scale_color_manual(name = '',values = c("red","skyblue")) +
theme_bw() + labs(y = "CRPS", x= "lead time [h]", colour= "Forecast") +
theme(plot.title = element_text(hjust = 0.5)) + geom_line()
all_emos
sale <- ggplot(aes(x=lead, y=value, colour=variable)) + ggtitle("EMOS postprocesamiento") +
scale_color_manual(name = '',values = c("red","skyblue")) +
theme_bw() + labs(y = "CRPS", x= "lead time [h]", colour= "Forecast") +
theme(plot.title = element_text(hjust = 0.5)) + geom_line()
sale <- ggplot(data = all_emos, aes(x=lead, y=value, colour=variable)) + ggtitle("EMOS postprocesamiento") +
scale_color_manual(name = '',values = c("red","skyblue")) +
theme_bw() + labs(y = "CRPS", x= "lead time [h]", colour= "Forecast") +
theme(plot.title = element_text(hjust = 0.5)) + geom_line()
print(sale)
all_emos %>%
dplyr::mutate(
CRPS = SpecsVerification::GaussCrps(IFS_T_mea, IFS_T_std, OBS),
CRPS_p = SpecsVerification::GaussCrps(IFS_T_mea_p, IFS_T_std_p, OBS),
) %>%
dplyr::group_by(lead) %>%
dplyr::summarize(DMO = mean(CRPS,na.rm=T), EMOS = mean(CRPS_p,na.rm=T)) %>%
tidyr::gather(variable, value, -lead)
sale <- ggplot(data = all_emos, aes(x=lead, y=value, colour=variable)) + ggtitle("EMOS postprocesamiento") +
scale_color_manual(name = '',values = c("red","skyblue")) +
theme_bw() + labs(y = "CRPS", x= "lead time [h]", colour= "Forecast") +
theme(plot.title = element_text(hjust = 0.5)) + geom_line()
print(sale)
q()
library(aws.wrfsmn)
knitr::opts_chunk$set(echo = TRUE)
get.wrf.files(year = 2023, month = 1, day = 1, cycle = 12, time = "01H")
for (i in 1:5)
get.wrf.files(year = 2023, month = 1, day = 1, cycle = 12, time = "01H")
}
for (i in 1:5)
{get.wrf.files(year = 2023, month = 1, day = 1, cycle = 12, time = "01H")}
for (i in 1:5)
{get.wrf.files(year = 2023, month = 1, day = i, cycle = 12, time = "01H")}
filenames <- c()
for (i in 1:5)
{aux.filenames <- get.wrf.files(year = 2023, month = 1, day = i, cycle = 12, time = "01H")
filenames <- c(filenames, aux.filenames)
}
73*5
filenames
q()
install.packages("parsnip")
library("parsnip")
library("tinymodels")
install.packages("tinymodels")
install.packages("modeldata")
library("modeldata")
data(ames)
ames
tidymodels_prefer()
install.packages("tidymodels")
library("tidymodels")
tidymodels_prefer()
ames
linear_reg %>% set_engine)"lm"
linear_reg %>% set_engine("lm")
linear_reg() %>% set_engine("lm")
linear_reg() %>% set_engine("glmnet")
linear_reg()
lm_model <- linear_reg() %>%
set_engine("lm")
lm_model
lm_form_fit <- lm_model %>%
fit(Sale_Price ~ Longitude + Latitude, data = ames)
lm_form_fit
rand_forest(trees = 1000, min_n = 5) %>%
set_engine("ranger") %>%
set_mode("regression") %>%
translate()
rand_forest(trees = 1000, min_n = 5) %>%
set_engine("ranger", verbose = TRUE) %>%
set_mode("regression")
?set_engine
show_engines()
f <- rand_forest(trees = 1000, min_n = 5) %>%
set_engine("ranger", verbose = TRUE) %>%
set_mode("regression")
f
show_engines(f)
show_engines("rand_forest")
show_engines("mars")
show_engines("linear_reg")
lm_form_fit
lm_form_fit %>% extract_fit_engine()
lm_form_fit
lm_form_fit$fit
predict(lm_form_fit$fit)
predict(lm_form_fit)
predict(lm_form_fit)
parsnio::predict(lm_form_fit)
parsnip::predict(lm_form_fit)
model_res <- lm_form_fit %>%
extract_fit_engine() %>%
summary()
model_res
coef(model_res)
lm_form_fit
tidy(lm_form_fit)
predict(lm_form_fit, new_data = salida)
salida <- ames %>% slice(1:5)
salida
salida <- ames %>% slice(1:4)
salida
salida <- ames %>% slice(1:5)
predict(lm_form_fit, new_data = salida)
salida <- ames$Sale_Price %>% slice(1:5)
ames$Sale_Price[1:5]
predict(lm_form_fit, new_data = salida)
salida
salida <- ames %>% slice(1:5)
predict(lm_form_fit, new_data = salida)
parsnip_addin()
parsnip_addin()
install.packages("shiny")
library("shiny")
parsnip_addin()
library("parsnip")
parsnip_addin()
install.packages("miniUI")
library("miniUI")
install.packages("rstudioapi")
library("miniUI")
library("shiny")
library("rstudioapi")
library("parsnip")
parsnip_addin()
parsnip_addin()
q()
install.packages("aws.wrfsmn")
library("aws.wrfsmn")
?ith
q()
q()
setwd("./Documentos/aws.wrfsmn/")
devtools::load_all()
?aws.wrfsmn::find.nearest.point()
lista <- list.files(path = "..", pattern = ".nc")
data <- rast(lista)
lista <- list.files(path = "..", pattern = ".nc", full.names = T)
data <- rast(lista)
data()
data
data.transform <- project(data, "+proj=longlat +datum=WGS84", method = "bilinear")
lista
load.by.variable(lista, "T2", TRUE, "bilinear")
data.by.ver <- load.by.variable(lista, "T2", TRUE, "bilinear")
data.by.ver
names(data)
unique(names(data))
usethis::use_r("daily_data_fields")
data.by.ver
time(data.by.ver)
times(data.by.ver)
time(data.by.ver)
data.by.ver
terra::time(data.by.ver)
library(lubridate)
day(terra::time(data.by.ver))
as.Date(terra::time(data.by.ver))
unique(as.Date(terra::time(data.by.ver)))
date.string <- unique(as.Date(terra::time(data.by.ver)))
which(as.Date(terra::time(data.by.ver)) == date.string[1])
which(as.Date(terra::time(data.by.ver)) == date.string[2])
s <- which(as.Date(terra::time(data.by.ver)) == date.string[2])
data.by.ver[[s]]
aggregate(data.by.ver[[s]], fun = "mean")
aggregate(data.by.ver[[s]], fun = mean)
?aggregate
terra::aggregate(data.by.ver[[s]], fun = "mean")
s
data.by.ver
terra::aggregate(data.by.ver[[s]], fun = "mean")
plot(terra::aggregate(data.by.ver[[s]], fun = "mean"))
s
data.by.ver[[s]]
terra::aggregate(data.by.ver[[s]], fun = "sum")
terra::aggregate(data.by.ver[[s]], fun = sum)
sum(data.by.ver[[s]])
plot(sum(data.by.ver[[s]]))
terra::plot(sum(data.by.ver[[s]]))
terra::plot(mean(data.by.ver[[s]]))
terra::plot(mean(data.by.ver[[s]]))
terra::plot(sum(data.by.ver[[s]]))
mean(data.by.ver[[s]])
mean(data.by.ver[[s]], na.rm = T)
data.by.ver[[s]]
mean(data.by.ver[[s]])
?mean
max(data.by.ver[[s]])
plot(max(data.by.ver[[s]]))
terra::plot(max(data.by.ver[[s]]))
mean(data.by.ver[[s]])
typeof(values(data.by.ver))
library(raster)
typeof(values(data.by.ver))
mean(data.by.ver[[s]], na.rm = T)
terra::plot(mean(data.by.ver[[s]], na.rm = T))
summarize
summarize(data.by.ver)
devtools::load_all()
daily.data.fields(data.by.ver, "sum")
sum(data.by.ver[[s]])
devtools::load_all()
daily.data.fields(data.by.ver, "sum")
unique(as.Date(time(data.by.ver)))
data.by.ver
unique(as.Date(time(data.by.ver)))
as.Date(time(data.by.ver))
data.by.ver
time(data.by.ver)
devtools::load_all()
daily.data.fields(data.by.ver, "sum")
devtools::load_all()
daily.data.fields(data.by.ver, "sum")
salida <- daily.data.fields(data.by.ver, "sum")
salida
devtools::load_all()
salida <- daily.data.fields(data.by.ver, "sum")
salida
salida
plot(salida)
terra::plot(salida)
unlist(salida)
devtools::load_all()
salida <- daily.data.fields(data.by.ver, "sum")
salida
plot(output)
plot(salida)
names(salida)
names(salida) <- paste0(names(salida), "_", date.string)
names(salida)
devtools::load_all()
salida <- daily.data.fields(data.by.ver, "sum")
devtools::load_all()
salida <- daily.data.fields(data.by.ver, "sum")
plot(salida)
salida <- daily.data.fields(data.by.ver, "mean")
salida <- daily.data.fields(data.by.ver, "mean")
salida
sum(data.by.ver[[s]])
mean(data.by.ver[[s]])
class(data.by.ver[[s]])
typeof(values(data.by.ver))
mean(data.by.ver[[s]])
devtools::load_all()
salida <- daily.data.fields(data.by.ver, "mean")
devtools::load_all()
salida <- daily.data.fields(data.by.ver, "mean")
devtools::load_all()
salida <- daily.data.fields(data.by.ver, "mean")
values(data.by.ver)
values(data.by.ver)[1:10]
values(data.by.ver)[100:150]
which(data.by.ver == NaN)
data.by.ver[[data.by.ver == NaN]] <- NA
values(data.by.ver)[is.nan(values(data.by.ver))] <- NA
values(data.by.ver)[1:10]
?values
raster::values(data.by.ver)[1:10]
terra::values(data.by.ver)[1:10]
devtools::load_all()
salida <- daily.data.fields(data.by.ver, "mean")
salida
devtools::load_all()
salida <- daily.data.fields(data.by.ver, "mean")
devtools::load_all()
salida <- daily.data.fields(data.by.ver, "mean")
mean(data.by.ver)
terra::values(data.by.ver)[1:10]
devtools::load_all()
salida <- daily.data.fields(data.by.ver, "mean")
salida <- daily.data.fields(data.by.ver, "sum")
devtools::load_all()
salida <- daily.data.fields(data.by.ver, "mean")
salida
data.by.ver
plot(values(data.by.ver))
min(values(data.by.ver))
max(values(data.by.ver))
mean(values(data.by.ver), na.rm = T)
app(data.by.ver, fun = mean)
terra::app(data.by.ver, fun = mean)
plot(terra::app(data.by.ver, fun = mean))
terra::app(data.by.ver, fun = sum)
daily.data.fields(data.by.ver, "sum")
devtools::load_all()
salida <- daily.data.fields(data.by.ver, "mean")
library(terra)
salida <- daily.data.fields(data.by.ver, "mean")
salida
plot(salida)
salida <- daily.data.fields(data.by.ver, "sum")
salida
plot(salida)
devtools::load_all()
salida <- daily.data.fields(data.by.ver, "sum")
salida <- daily.data.fields(data.by.ver, "mean")
devtools::load_all()
salida <- daily.data.fields(data.by.ver, "mean")
?daily.data.fields
devtools::document()
?daily.data.fields
salida <- daily.data.fields(data.by.ver, "min")
salida
plot(salida)
q()