Skip to content

Commit

Permalink
Merge pull request #3 from rqthomas/salinity
Browse files Browse the repository at this point in the history 
Salinity
  • Loading branch information
@rqthomas
rqthomas authored Jul 23, 2022
2 parents 1e3fdb5 + a0eb41b commit ea165c9
Show file tree
Hide file tree
Showing 9 changed files with 156 additions and 100 deletions.
13 changes: 10 additions & 3 deletions R/edit_nml_functions.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -121,6 +121,7 @@ update_nml <- function(var_list,var_name_list,working_directory, nml){
#' add(10, 1)
get_glm_nc_var_all_wq <- function(ncFile,working_dir, z_out,vars_depth, vars_no_depth, diagnostic_vars){
glm_nc <- ncdf4::nc_open(paste0(working_dir, ncFile))
glm_vars <- names(glm_nc$var)
tallest_layer <- ncdf4::ncvar_get(glm_nc, "NS")
final_time_step <- length(tallest_layer)
tallest_layer <- tallest_layer[final_time_step]
Expand All @@ -129,9 +130,15 @@ get_glm_nc_var_all_wq <- function(ncFile,working_dir, z_out,vars_depth, vars_no_
heights <- heights[1:tallest_layer, final_time_step]
heights_out <- heights_surf - z_out

snow <- matrix(ncdf4::ncvar_get(glm_nc, "hsnow"), ncol = final_time_step)[final_time_step]
ice_white <- matrix(ncdf4::ncvar_get(glm_nc, "hwice"), ncol = final_time_step)[final_time_step]
ice_blue <- matrix(ncdf4::ncvar_get(glm_nc, "hice"), ncol = final_time_step)[final_time_step]
if("hsnow" %in% glm_vars){
snow <- matrix(ncdf4::ncvar_get(glm_nc, "hsnow"), ncol = final_time_step)[final_time_step]
ice_white <- matrix(ncdf4::ncvar_get(glm_nc, "hwice"), ncol = final_time_step)[final_time_step]
ice_blue <- matrix(ncdf4::ncvar_get(glm_nc, "hice"), ncol = final_time_step)[final_time_step]
}else{
snow <- matrix(ncdf4::ncvar_get(glm_nc, "snow_thickness"), ncol = final_time_step)[final_time_step]
ice_white <- matrix(ncdf4::ncvar_get(glm_nc, "white_ice_thickness"), ncol = final_time_step)[final_time_step]
ice_blue <- matrix(ncdf4::ncvar_get(glm_nc, "blue_ice_thickness"), ncol = final_time_step)[final_time_step]
}
avg_surf_temp <- matrix(ncdf4::ncvar_get(glm_nc, "avg_surf_temp"), ncol = final_time_step)[final_time_step]


Expand Down
28 changes: 15 additions & 13 deletions R/generate_forecast_score.R
View file
Original file line number Diff line number Diff line change
@@ -1,9 +1,3 @@
#' Generate forecast score from csv target and forecast files
#'
#' @param targets_file
#' @param forecast_file
#' @param output_directory
#'
#' @return
#' @export
#'
Expand All @@ -12,16 +6,22 @@ generate_forecast_score <- function(targets_file,
forecast_file,
output_directory){

target <- read_csv(targets_file, show_col_types = FALSE) |>
rename(z = depth) |>
mutate(x = NA,
y = NA,
target_id = "fcre")
target <- readr::read_csv(targets_file, show_col_types = FALSE) |>
dplyr::rename(z = depth) |>
dplyr::mutate(x = NA,
y = NA,
target_id = "fcre",
site_id = paste0(site_id,"-",z))

tools::file_path_sans_ext(basename(forecast_file))

file_name <- file.path(output_directory,paste0("score-",tools::file_path_sans_ext(basename(forecast_file)),".csv.gz"))

forecast_file %>%
read4cast::read_forecast(grouping_variables = c("time", "depth"),
target_variables = "temperature") %>%
dplyr::mutate(filename = forecast_file) %>%
dplyr::mutate(filename = forecast_file,
site_id = paste0(site_id,"-",depth)) %>%
#score4cast::select_forecasts() %>%
score4cast::pivot_forecast() %>%
score4cast::crps_logs_score(target) %>%
Expand All @@ -30,5 +30,7 @@ generate_forecast_score <- function(targets_file,
units = "seconds"),
horizon = horizon / 86400) |>
#score4cast::include_horizon() %>%
readr::write_csv(file.path(output_directory,paste0("score-",da_forecast_output$save_file_name_short,".csv.gz")))
readr::write_csv(file_name)

invisible(file_name)
}
9 changes: 4 additions & 5 deletions R/generate_initial_conditions.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,14 +21,13 @@ generate_initial_conditions <- function(states_config,
obs,
config,
historical_met_error = FALSE){

pars_config <- pars_config[pars_config$model == config$model_settings$model, ]

if(is.na(config$run_config$restart_file)){

init <- list()
if(!is.null(pars_config) & any(pars_config$model == config$model_settings$model)){
pars_config <- pars_config[pars_config$model == config$model_settings$model, ]
if(!is.null(pars_config)){
if("model" %in% names(pars_config)){
pars_config <- pars_config[pars_config$model == config$model_settings$model, ]
}
npars <- nrow(pars_config)
}else{
npars <- 0
Expand Down
4 changes: 2 additions & 2 deletions R/plotting_general.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -261,7 +261,7 @@ plotting_general <- function(file_name,
print(p)
}

if("extc_coef" %in% diagnostics_names){
if(("extc_coef" %in% diagnostics_names) | ("extc" %in% diagnostics_names)){

message("secchi")

Expand All @@ -278,7 +278,7 @@ plotting_general <- function(file_name,
obs_curr <- rep(NA, length(obs_date))
}

i <- which(diagnostics_names == "extc_coef")
i <- which(diagnostics_names == "extc_coef" | diagnostics_names == "extc")
ii <- which.min(abs(depths-1.0))
curr_var <- diagnostic_list[[i]]

Expand Down
4 changes: 2 additions & 2 deletions R/plotting_general_2.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -260,7 +260,7 @@ plotting_general_2 <- function(file_name,
print(p)
}

if("extc_coef" %in% diagnostics_names){
if(("extc_coef" %in% diagnostics_names) | ("extc" %in% diagnostics_names)){

message("secchi")

Expand All @@ -279,7 +279,7 @@ plotting_general_2 <- function(file_name,
obs_curr <- rep(NA, length(obs_date))
}

i <- which(diagnostics_names == "extc_coef")
i <- which(diagnostics_names == "extc_coef" | diagnostics_names == "extc")
ii <- which.min(abs(depths-1.0))
curr_var <- diagnostic_list[[i]]

Expand Down
6 changes: 3 additions & 3 deletions R/run_da_forecast.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -459,7 +459,7 @@ run_da_forecast <- function(states_init,
da_qc_flag[i] <- 0
}

if(npars > 0){


x[i, , , ] <- x_corr
if(npars > 0) pars[i, , ] <- pars_star
Expand All @@ -472,14 +472,14 @@ run_da_forecast <- function(states_init,
}

if(i == (hist_days + 1) & config$uncertainty$initial_condition == FALSE){
pars[i, , ] <- pars_star
if(npars > 0) pars[i, , ] <- pars_star
for(s in 1:nstates){
for(k in 1:ndepths_modeled){
x[i, s, k , ] <- mean(x_star[s, k, ])
}
}
}
}


for(s in 1:nstates){
for(m in 1:nmembers){
Expand Down
2 changes: 1 addition & 1 deletion R/set_up_model.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,7 +40,7 @@ set_up_model <- function(config,
to = file.path(ens_working_directory, "aed2.nml"), overwrite = TRUE)

file.copy(from = file.path(config$file_path$configuration_directory, config$model_settings$base_AED_phyto_pars_nml),
to = file.path(ens_working_directory, "aed2_phyto_pars.nml"), overwrite = TRUE)
to = file.path(ens_working_directory, "aed_phyto_pars.csv"), overwrite = TRUE)

file.copy(from = file.path(config$file_path$configuration_directory, config$model_settings$base_AED_zoop_pars_nml),
to = file.path(ens_working_directory, "aed2_zoop_pars.nml"), overwrite = TRUE)
Expand Down
186 changes: 117 additions & 69 deletions R/write_forecast_csv.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,116 +36,164 @@ write_forecast_csv <- function(da_forecast_output,

forecast_flag[which(is.na(forecast_flag))] <- 0

output_list <- NULL
for(i in 1:dim(x)[1]){
for(j in 1:dim(x)[2]){
for(k in 1:dim(x)[3]){
tmp <- tibble(predicted = x[i, j, k, ],
time = full_time[i],
depth = config$model_settings$modeled_depths[k],
variable = states_config$state_names[j],
forecast = forecast_flag[i],
ensemble = 1:dim(x)[4],
variable_type = "state")
output_list <- rbind(output_list, tmp)
indexes <- expand.grid(1:dim(x)[1], 1:dim(x)[2], 1:dim(x)[3])
ensembles <- 1:dim(x)[4]

output_list <- map_dfr(1:nrow(indexes), function(i, indexes){
var1 <- indexes$Var1[i]
var2 <- indexes$Var2[i]
var3 <- indexes$Var3[i]
tibble::tibble(predicted = x[var1, var2, var3, ],
time = full_time[var1],
depth = config$model_settings$modeled_depths[var3],
variable = states_config$state_names[var2],
forecast = forecast_flag[var1],
ensemble = ensembles,
variable_type = "state")
},
indexes = indexes
)

tmp_index <- 0
for(s in 1:length(obs_config$state_names_obs)){
if(!obs_config$state_names_obs[s] %in% states_config$state_names){
tmp_index <- tmp_index + 1
first_index <- 1
for(ii in 1:length(states_config$state_names)){
if(s %in% states_config$states_to_obs[[ii]]){
temp_index <- which(states_config$states_to_obs[[ii]] == s)
if(first_index == 1){
temp_var <- x[ ,ii , , ] * states_config$states_to_obs_mapping[[ii]][temp_index]
first_index <- 2
}else{
temp_var <- temp_var + x[, ii, , ] * states_config$states_to_obs_mapping[[ii]][temp_index]
}
}
}

indexes <- expand.grid(1:dim(temp_var)[1], 1:dim(temp_var)[2])

output_list_tmp <- map_dfr(1:nrow(indexes), function(i, indexes){
var1 <- indexes$Var1[i]
var3 <- indexes$Var2[i]
tibble::tibble(predicted = temp_var[var1, var3, ],
time = full_time[var1],
depth = config$model_settings$modeled_depths[var3],
variable = obs_config$target_variable[s],
forecast = forecast_flag[var1],
ensemble = ensembles,
variable_type = "state")
},
indexes = indexes
)
output_list <- bind_rows(output_list, output_list_tmp)
}
}


if(length(config$output_settings$diagnostics_names) > 0){
for(i in 1:dim(diagnostics)[1]){
for(j in 1:dim(diagnostics)[2]){
for(k in 1:dim(diagnostics)[3]){
tmp <- tibble(predicted = diagnostics[i, j, k, ],
time = full_time[j],
variable = states_config$state_names[i],
depth = config$model_settings$modeled_depths[k],
forecast = forecast_flag[j],
ensemble = 1:dim(diagnostics)[4],
variable_type = "diagnostic")
output_list <- rbind(output_list, tmp)
}
}
}

indexes <- expand.grid(1:dim(diagnostics)[1], 1:dim(diagnostics)[2], 1:dim(diagnostics)[3])

output_list2 <- map_dfr(1:nrow(indexes), function(i, indexes){
var1 <- indexes$Var1[i]
var2 <- indexes$Var2[i]
var3 <- indexes$Var3[i]
tibble::tibble(predicted = diagnostics[var1, var2, var3, ],
time = full_time[var1],
depth = config$model_settings$modeled_depths[var3],
variable = states_config$state_names[var2],
forecast = forecast_flag[var1],
ensemble = ensembles,
variable_type = "diagnostic")
},
indexes = indexes
)
output_list <- dplyr::bind_rows(output_list, output_list2)
}

for(i in 1:dim(pars)[1]){
for(j in 1:dim(pars)[2]){
tmp <- tibble(predicted = pars[i, j, ],
time = full_time[i],
variable = pars_config$par_names_save[j],
depth = NA,
forecast = forecast_flag[i],
ensemble = 1:dim(pars)[3],
variable_type = "parameter")
output_list <- rbind(output_list, tmp)
}
if(!is.null(pars)){
indexes <- expand.grid(1:dim(pars)[1], 1:dim(pars)[2])

output_list3 <- map_dfr(1:nrow(indexes), function(i, indexes){
var1 <- indexes$Var1[i]
var2 <- indexes$Var2[i]
tibble::tibble(predicted = pars[var1, var2, ],
time = full_time[var1],
depth = NA,
variable = pars_config$par_names_save[var2],
forecast = forecast_flag[var1],
ensemble = ensembles,
variable_type = "parameter")
},
indexes = indexes
)
output_list <- dplyr::bind_rows(output_list, output_list3)
}

if(!is.null(da_forecast_output$restart_list)){
lake_depth <- da_forecast_output$restart_list$lake_depth
}
for(i in 1:dim(lake_depth)[1]){
tmp <- tibble(predicted = lake_depth[i, ],
time = full_time[i],
variable = "depth",
depth = NA,
forecast = forecast_flag[i],
ensemble = 1:dim(lake_depth)[2],
variable_type = "state")
output_list <- rbind(output_list, tmp)
tmp <- tibble::tibble(predicted = lake_depth[i, ],
time = full_time[i],
variable = "depth",
depth = NA,
forecast = forecast_flag[i],
ensemble = 1:dim(lake_depth)[2],
variable_type = "state")
output_list <- dplyr::bind_rows(output_list, tmp)
}


if(length(config$output_settings$diagnostics_names) > 0){
for(i in 1:dim(diagnostics)[2]){
tmp <- tibble(predicted = 1.7 / diagnostics[1, i, which.min(abs(config$model_settings$modeled_depths-1.0)), ],
time = full_time[i],
variable = "secchi",
depth = NA,
forecast = forecast_flag[i],
ensemble = 1:dim(diagnostics)[4],
variable_type = "state")
output_list <- rbind(output_list, tmp)
tmp <- tibble::tibble(predicted = 1.7 / diagnostics[1, i, which.min(abs(config$model_settings$modeled_depths-1.0)), ],
time = full_time[i],
variable = "secchi",
depth = NA,
forecast = forecast_flag[i],
ensemble = 1:dim(diagnostics)[4],
variable_type = "state")
output_list <- dplyr::bind_rows(output_list, tmp)
}

}

for(i in 1:dim(snow_ice_thickness)[1]){
tmp <- tibble(predicted = apply(snow_ice_thickness[2:3, i, ], 2, sum),
time = full_time[i],
variable = "ice_thickness",
depth = NA,
forecast = forecast_flag[i],
ensemble = 1:dim(snow_ice_thickness)[3],
variable_type = "state")
output_list <- rbind(output_list, tmp)
tmp <- tibble::tibble(predicted = apply(snow_ice_thickness[2:3, i, ], 2, sum),
time = full_time[i],
variable = "ice_thickness",
depth = NA,
forecast = forecast_flag[i],
ensemble = 1:dim(snow_ice_thickness)[3],
variable_type = "state")
output_list <- dplyr::bind_rows(output_list, tmp)
}

time_of_forecast <- lubridate::with_tz(da_forecast_output$time_of_forecast, tzone = "UTC")

output_list <- output_list |>
mutate(pub_time = time_of_forecast,
start_time = forecast_start_datetime,
site_id = config$location$site_id,
model_id = config$run_config$sim_name) |>
select(start_time, pub_time, model_id, site_id, depth, time, ensemble, variable, predicted, forecast, variable_type)
dplyr::mutate(pub_time = time_of_forecast,
start_time = forecast_start_datetime,
site_id = config$location$site_id,
model_id = config$run_config$sim_name) |>
dplyr::select(start_time, pub_time, model_id, site_id, depth, time, ensemble, variable, predicted, forecast, variable_type)

if(!use_short_filename | is.na(da_forecast_output$save_file_name_short) | length(which(forecast_flag == 1)) == 0){
fname <- file.path(forecast_output_directory, paste0(da_forecast_output$save_file_name,".csv.gz"))
}else{
fname <- file.path(forecast_output_directory, paste0(da_forecast_output$save_file_name_short,".csv.gz"))
}

#Convert to target variable name
for(i in 1:length(states_config$state_names)){
if(length(which(obs_config$state_names_obs == states_config$state_names[i])) >0){
obs_name <- obs_config$target_variable[which(obs_config$state_names_obs == states_config$state_names[i])]
output_list <- output_list %>%
mutate(variable = ifelse(variable == states_config$state_names[i], obs_name, variable))
dplyr::mutate(variable = ifelse(variable == states_config$state_names[i], obs_name, variable))
}
}

write_csv(output_list, fname)
invisible(fname)

}
Loading

0 comments on commit ea165c9

Please sign in to comment.