From caa0cdba55222a3b97cc40d96b7dae370fa4084e Mon Sep 17 00:00:00 2001 From: dschlaep Date: Sat, 24 Sep 2016 15:55:29 +0200 Subject: [PATCH] Faster input/output MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Improvements: - use package ‘iotools' when available for reading files - use ‘cat’ instead of wrapper function ‘write’ to write to files - new functions ‘swsf_read_csv’, ‘swsf_read_inputfile’, and ‘reconstitute_inputfile’ to handle reading and writing of data input files - all objects ‘sw_input_XXX_use’ are now named logical vectors - ‘aon’, ‘exinfo’, and ‘places’ are now named logical lists - part 1 gains new setting ‘eta.estimate’: if TRUE, then ‘estimate time of arrival’ of simulations is calculated as mean and (new) 95% prediction interval - re-wrote code to read in data files - re-wrote code of ‘do_OneSite’ in section ‘create’ for ‘LookupTranspCoeffFromTable_XXX’: loop instead of copy paste - re-wrote code of ‘do_OneSite’ in section ‘create’ for ‘adding vegetation information’ - new function ‘update_biomass’ replaces ‘biomassComponents’ Test projects: - now deletes also log files when running with MPI - script ‘Run_all_test_projects.R’ can now accept command line inputs when run non-interactively (run ‘r -f Run_all_test_projects.R --arg —help’ for explanations) Bug fixes: - fixing bug in ‘export_objects_to_workers’ for downscaling climate scenarios and extracting soil information - fixing bug in ‘dw_NCEPCFSR_Global’: incorrect year information (copy-paste mistake) Former-commit-id: 0af8751c023587315232f857da3ee5ad7e3d8a78 [formerly 2f7eb0c3d36f96cbe25f02381097fbaafe3995a5] Former-commit-id: bc5057df996bd4af295db6c4d6695ca6f7227249 --- 2_SWSF_p1of5_Settings_v51.R | 1 + 2_SWSF_p2of5_CreateDB_Tables_v51.R | 2 +- 2_SWSF_p3of5_ExternalDataExtractions_v51.R | 88 +- 2_SWSF_p4of5_Code_v51.R | 920 +++++++++------------ 2_SWSF_p5of5_Functions_v51.R | 66 +- DESCRIPTION | 5 +- 6 files changed, 531 insertions(+), 551 deletions(-) diff --git a/2_SWSF_p1of5_Settings_v51.R b/2_SWSF_p1of5_Settings_v51.R index 798aeb2b..5d8eeab1 100644 --- a/2_SWSF_p1of5_Settings_v51.R +++ b/2_SWSF_p1of5_Settings_v51.R @@ -35,6 +35,7 @@ rm(list=ls(all=TRUE)) #------Overall timing t.overall <- Sys.time() be.quiet <- FALSE +eta.estimate <- interactive() print.debug <- interactive() debug.warn.level <- sum(c(print.debug, interactive())) debug.dump.objects <- interactive() diff --git a/2_SWSF_p2of5_CreateDB_Tables_v51.R b/2_SWSF_p2of5_CreateDB_Tables_v51.R index 3a3ccd80..0801c205 100644 --- a/2_SWSF_p2of5_CreateDB_Tables_v51.R +++ b/2_SWSF_p2of5_CreateDB_Tables_v51.R @@ -240,7 +240,7 @@ if (length(Tables) == 0 || do.clean) { #first add any from the experimentals table if its turned on #next add any from the treatments table if its turned on treatments_lookupweatherfolders <- character(0) - if(any(names(sw_input_treatments_use[-1][which(sw_input_treatments_use[-1] > 0 & is.finite(as.numeric(sw_input_treatments_use[-1])))])=="LookupWeatherFolder")) { + if(any(names(sw_input_treatments_use[sw_input_treatments_use])=="LookupWeatherFolder")) { treatments_lookupweatherfolders <- c(treatments_lookupweatherfolders, sw_input_treatments$LookupWeatherFolder[runIDs_sites]) } if(any(create_experimentals=="LookupWeatherFolder")) { diff --git a/2_SWSF_p3of5_ExternalDataExtractions_v51.R b/2_SWSF_p3of5_ExternalDataExtractions_v51.R index 0d651e1a..f255775b 100644 --- a/2_SWSF_p3of5_ExternalDataExtractions_v51.R +++ b/2_SWSF_p3of5_ExternalDataExtractions_v51.R @@ -2116,7 +2116,9 @@ if (exinfo$GDODCPUCLLNL || exinfo$ExtractClimateChangeScenarios_CMIP5_BCSD_NEX_U # extract the GCM data depending on parallel backend if (identical(parallel_backend, "mpi")) { - export_objects_to_workers(list.export, list(parent = parent.frame()), "mpi") + export_objects_to_workers(list.export, + list(local = environment(), parent = parent.frame(), global = .GlobalEnv), + "mpi") if (is_NEX && useRCurl && !saveNEXtempfiles) Rmpi::mpi.bcast.cmd(library("RCurl", quietly=TRUE)) if (is_GDODCPUCLLNL) @@ -2147,7 +2149,9 @@ if (exinfo$GDODCPUCLLNL || exinfo$ExtractClimateChangeScenarios_CMIP5_BCSD_NEX_U Rmpi::mpi.bcast.cmd(gc()) } else if (identical(parallel_backend, "snow")) { - export_objects_to_workers(list.export, list(parent = parent.frame()), "snow", cl) + export_objects_to_workers(list.export, + list(local = environment(), parent = parent.frame(), global = .GlobalEnv), + "snow", cl) if (is_NEX && useRCurl && !saveNEXtempfiles) snow::clusterEvalQ(cl, library("RCurl", quietly = TRUE)) @@ -2637,18 +2641,22 @@ if ( exinfo$ExtractClimateChangeScenarios_CMIP3_ClimateWizardEnsembles_Global || #temp value in C #ppt value in mm #add data to sw_input_climscen and set the use flags - sw_input_climscen_values_use[i.temp <- match(paste("PPTmm_m", st_mo, "_sc", formatC(sc, width=2,format="d", flag="0"), sep=""), colnames(sw_input_climscen_values_use))] <- 1 + i.temp <- paste0("PPTmm_m", st_mo, "_sc", formatC(sc, width=2,format="d", flag="0")) + sw_input_climscen_values_use[i.temp] <- TRUE sw_input_climscen_values[, i.temp] <- sc.ppt - sw_input_climscen_values_use[i.temp <- match(paste("TempC_m", st_mo, "_sc", formatC(sc, width=2,format="d", flag="0"), sep=""), colnames(sw_input_climscen_values_use))] <- 1 + i.temp <- paste0("TempC_m", st_mo, "_sc", formatC(sc, width=2,format="d", flag="0")) + sw_input_climscen_values_use[i.temp] <- TRUE sw_input_climscen_values[, i.temp] <- sc.temp } if (exinfo$ExtractClimateChangeScenarios_CMIP3_ClimateWizardEnsembles_USA) { sc.temp <- sc.temp * 5/9 #temp addand in C sc.ppt <- 1 + sc.ppt/100 #ppt change as factor #add data to sw_input_climscen and set the use flags - sw_input_climscen_use[i.temp <- match(paste("PPTfactor_m", st_mo, "_sc", formatC(sc, width=2,format="d", flag="0"), sep=""), colnames(sw_input_climscen_use))] <- 1 + i.temp <- paste0("PPTfactor_m", st_mo, "_sc", formatC(sc, width=2,format="d", flag="0")) + sw_input_climscen_use[i.temp] <- TRUE sw_input_climscen[, i.temp] <- sc.ppt - sw_input_climscen_use[i.temp <- match(paste("deltaTempC_m", st_mo, "_sc", formatC(sc, width=2,format="d", flag="0"), sep=""), colnames(sw_input_climscen_use))] <- 1 + i.temp <- paste0("deltaTempC_m", st_mo, "_sc", formatC(sc, width=2,format="d", flag="0")) + sw_input_climscen_use[i.temp] <- TRUE sw_input_climscen[, i.temp] <- sc.temp } } @@ -2657,11 +2665,11 @@ if ( exinfo$ExtractClimateChangeScenarios_CMIP3_ClimateWizardEnsembles_Global || if (res > 0) print(paste(res, "sites didn't extract climate scenario information by 'ExtractClimateChangeScenarios_CMIP3_ClimateWizardEnsembles'")) #write data to datafile.climatescenarios_values - tempdat <- rbind(sw_input_climscen_values_use, sw_input_climscen_values) - write.csv(tempdat, file=file.path(dir.sw.dat, datafile.climatescenarios_values), row.names=FALSE) + write.csv(reconstitute_inputfile(sw_input_climscen_values_use, sw_input_climscen_values), + file = file.path(dir.sw.dat, datafile.climatescenarios_values), row.names = FALSE) unlink(file.path(dir.in, datafile.SWRWinputs_preprocessed)) - rm(list.scenarios.datafile, list.scenarios.external, tempdat, sc.temp, sc.ppt, res, locations) + rm(list.scenarios.datafile, list.scenarios.external, sc.temp, sc.ppt, res, locations) } else { print("Not all scenarios requested in 'datafile.SWRunInformation' are available in with 'ExtractClimateChangeScenarios_CMIP3_ClimateWizardEnsembles'") } @@ -2850,23 +2858,27 @@ if (exinfo$ExtractSoilDataFromCONUSSOILFromSTATSGO_USA || exinfo$ExtractSoilData #add data to sw_input_soils and set the use flags i.temp <- grep("Matricd_L", names(sw_input_soils_use)) sw_input_soils[runIDs_sites[i_Done], i.temp[lys]] <- soil_data[i_good, lys, "matricd"] - sw_input_soils_use[i.temp][lys] <- 1 + sw_input_soils_use[i.temp[lys]] <- TRUE + sw_input_soils_use[i.temp[-lys]] <- FALSE i.temp <- grep("GravelContent_L", names(sw_input_soils_use)) sw_input_soils[runIDs_sites[i_Done], i.temp[lys]] <- soil_data[i_good, lys, "rockvol"] - sw_input_soils_use[i.temp][lys] <- 1 + sw_input_soils_use[i.temp[lys]] <- TRUE + sw_input_soils_use[i.temp[-lys]] <- FALSE i.temp <- grep("Sand_L", names(sw_input_soils_use)) sw_input_soils[runIDs_sites[i_Done], i.temp[lys]] <- soil_data[i_good, lys, "sand"] - sw_input_soils_use[i.temp][lys] <- 1 + sw_input_soils_use[i.temp[lys]] <- TRUE + sw_input_soils_use[i.temp[-lys]] <- FALSE i.temp <- grep("Clay_L", names(sw_input_soils_use)) sw_input_soils[runIDs_sites[i_Done], i.temp[lys]] <- soil_data[i_good, lys, "clay"] - sw_input_soils_use[i.temp][lys] <- 1 + sw_input_soils_use[i.temp[lys]] <- TRUE + sw_input_soils_use[i.temp[-lys]] <- FALSE #write data to datafile.soils - tempdat <- rbind(sw_input_soils_use, sw_input_soils) - write.csv(tempdat, file=file.path(dir.sw.dat, datafile.soils), row.names=FALSE) + write.csv(reconstitute_inputfile(sw_input_soils_use, sw_input_soils), + file = file.path(dir.sw.dat, datafile.soils), row.names = FALSE) unlink(file.path(dir.in, datafile.SWRWinputs_preprocessed)) - rm(tempdat, i.temp, i_Done) + rm(i.temp, i_Done) } if (!be.quiet) @@ -2966,7 +2978,9 @@ if (exinfo$ExtractSoilDataFromCONUSSOILFromSTATSGO_USA || exinfo$ExtractSoilData #call the simulations depending on parallel backend if (identical(parallel_backend, "mpi")) { - export_objects_to_workers(list.export, list(parent = parent.frame()), "mpi") + export_objects_to_workers(list.export, + list(local = environment(), parent = parent.frame(), global = .GlobalEnv), + "mpi") Rmpi::mpi.bcast.cmd(library(raster, quietly=TRUE)) sim_cells_SUIDs <- Rmpi::mpi.applyLB(x=is_ToDo, fun=extract_SUIDs, res = cell_res_wise, grid = grid_wise, sp_sites = run_sites_wise) @@ -2976,7 +2990,9 @@ if (exinfo$ExtractSoilDataFromCONUSSOILFromSTATSGO_USA || exinfo$ExtractSoilData Rmpi::mpi.bcast.cmd(gc()) } else if (identical(parallel_backend, "snow")) { - export_objects_to_workers(list.export, list(parent = parent.frame()), "snow", cl) + export_objects_to_workers(list.export, + list(local = environment(), parent = parent.frame(), global = .GlobalEnv), + "snow", cl) snow::clusterEvalQ(cl, library(raster, quietly = TRUE)) sim_cells_SUIDs <- snow::clusterApplyLB(cl, x=is_ToDo, fun=extract_SUIDs, res = cell_res_wise, grid = grid_wise, sp_sites = run_sites_wise) @@ -3183,23 +3199,27 @@ if (exinfo$ExtractSoilDataFromCONUSSOILFromSTATSGO_USA || exinfo$ExtractSoilData #add data to sw_input_soils and set the use flags i.temp <- grep("Matricd_L", names(sw_input_soils_use)) sw_input_soils[runIDs_sites[i_Done], i.temp[lys]] <- round(sim_cells_soils[i_good, paste0("bulk_L", lys)], 2) - sw_input_soils_use[i.temp][lys] <- 1 + sw_input_soils_use[i.temp[lys]] <- TRUE + sw_input_soils_use[i.temp[-lys]] <- FALSE i.temp <- grep("GravelContent_L", names(sw_input_soils_use)) sw_input_soils[runIDs_sites[i_Done], i.temp[lys]] <- round(sim_cells_soils[i_good, paste0("cfrag_L", lys)]) / 100 - sw_input_soils_use[i.temp][lys] <- 1 + sw_input_soils_use[i.temp[lys]] <- TRUE + sw_input_soils_use[i.temp[-lys]] <- FALSE i.temp <- grep("Sand_L", names(sw_input_soils_use)) sw_input_soils[runIDs_sites[i_Done], i.temp[lys]] <- round(sim_cells_soils[i_good, paste0("sand_L", lys)]) / 100 - sw_input_soils_use[i.temp][lys] <- 1 + sw_input_soils_use[i.temp[lys]] <- TRUE + sw_input_soils_use[i.temp[-lys]] <- FALSE i.temp <- grep("Clay_L", names(sw_input_soils_use)) sw_input_soils[runIDs_sites[i_Done], i.temp[lys]] <- round(sim_cells_soils[i_good, paste0("clay_L", lys)]) / 100 - sw_input_soils_use[i.temp][lys] <- 1 + sw_input_soils_use[i.temp[lys]] <- TRUE + sw_input_soils_use[i.temp[-lys]] <- FALSE #write data to datafile.soils - tempdat <- rbind(sw_input_soils_use, sw_input_soils) - write.csv(tempdat, file=file.path(dir.sw.dat, datafile.soils), row.names=FALSE) + write.csv(reconstitute_inputfile(sw_input_soils_use, sw_input_soils), + file = file.path(dir.sw.dat, datafile.soils), row.names = FALSE) unlink(file.path(dir.in, datafile.SWRWinputs_preprocessed)) - rm(lys, tempdat, i.temp, i_Done) + rm(lys, i.temp, i_Done) } if (!be.quiet) print(paste("'ExtractSoilDataFromISRICWISEv12_Global' was extracted for n =", sum(i_good), "out of", sum(do_extract[[2]]), "sites")) @@ -3590,17 +3610,18 @@ if (exinfo$ExtractSkyDataFromNOAAClimateAtlas_USA || exinfo$ExtractSkyDataFromNC #add data to sw_input_cloud and set the use flags i.temp <- grep("RH", names(sw_input_cloud_use)) - sw_input_cloud_use[i.temp] <- 1 + sw_input_cloud_use[i.temp] <- TRUE sw_input_cloud[runIDs_sites[i_good], i.temp[st_mo]] <- round(monthlyclim[i_good, "RH", ], 2) i.temp <- grep("SkyC", names(sw_input_cloud_use)) - sw_input_cloud_use[i.temp] <- 1 + sw_input_cloud_use[i.temp] <- TRUE sw_input_cloud[runIDs_sites[i_good], i.temp[st_mo]] <- round(monthlyclim[i_good, "cover", ], 2) i.temp <- grep("wind", names(sw_input_cloud_use)) - sw_input_cloud_use[i.temp] <- 1 + sw_input_cloud_use[i.temp] <- TRUE sw_input_cloud[runIDs_sites[i_good], i.temp[st_mo]] <- round(monthlyclim[i_good, "wind", ], 2) #write data to datafile.cloud - write.csv(rbind(sw_input_cloud_use, sw_input_cloud), file = file.path(dir.sw.dat, datafile.cloud), row.names = FALSE) + write.csv(reconstitute_inputfile(sw_input_cloud_use, sw_input_cloud), + file = file.path(dir.sw.dat, datafile.cloud), row.names = FALSE) unlink(file.path(dir.in, datafile.SWRWinputs_preprocessed)) rm(i.temp) @@ -3670,17 +3691,18 @@ if (exinfo$ExtractSkyDataFromNOAAClimateAtlas_USA || exinfo$ExtractSkyDataFromNC #add data to sw_input_cloud and set the use flags i.temp <- grep("RH", names(sw_input_cloud_use)) - sw_input_cloud_use[i.temp] <- 1 + sw_input_cloud_use[i.temp] <- TRUE sw_input_cloud[runIDs_sites[i_good], i.temp][, st_mo] <- round(monthlyclim[i_good, "RH", ], 2) i.temp <- grep("SkyC", names(sw_input_cloud_use)) - sw_input_cloud_use[i.temp] <- 1 + sw_input_cloud_use[i.temp] <- TRUE sw_input_cloud[runIDs_sites[i_good], i.temp][, st_mo] <- round(monthlyclim[i_good, "cover", ], 2) i.temp <- grep("wind", names(sw_input_cloud_use)) - sw_input_cloud_use[i.temp] <- 1 + sw_input_cloud_use[i.temp] <- TRUE sw_input_cloud[runIDs_sites[i_good], i.temp][, st_mo] <- round(monthlyclim[i_good, "wind", ], 2) #write data to datafile.cloud - write.csv(rbind(sw_input_cloud_use, sw_input_cloud), file=file.path(dir.sw.dat, datafile.cloud), row.names=FALSE) + write.csv(reconstitute_inputfile(sw_input_cloud_use, sw_input_cloud), + file = file.path(dir.sw.dat, datafile.cloud), row.names = FALSE) unlink(file.path(dir.in, datafile.SWRWinputs_preprocessed)) rm(i.temp) diff --git a/2_SWSF_p4of5_Code_v51.R b/2_SWSF_p4of5_Code_v51.R index d8dec0c8..316d8598 100644 --- a/2_SWSF_p4of5_Code_v51.R +++ b/2_SWSF_p4of5_Code_v51.R @@ -58,22 +58,18 @@ dirname.sw.runs.weather <- "WeatherData" #timing: basis for estimated time of arrival, ETA timerfile <- "temp_timer.csv" -temp <- file.path(dir.out, timerfile) -if (file.exists(temp) && !continueAfterAbort) { - try(file.remove(temp), silent=TRUE) -} -if (!file.exists(temp)) { - write.table(t(c(0,NA)), file=temp, append=TRUE, sep=",", dec=".", col.names=FALSE, row.names=FALSE) - runIDs_done <- NULL -} else { - ttemp <- read.csv(file = temp, header = FALSE) - runIDs_done <- if (nrow(ttemp) > 1) sort(ttemp[-1, 1]) else NULL - rm(ttemp) -} +ftemp <- file.path(dir.out, timerfile) +times <- if (!file.exists(ftemp) || !continueAfterAbort) { + cat("0,NA", file = ftemp, sep = "\n") + } else { + swsf_read_csv(file = ftemp, header = FALSE, colClasses = c("NULL", "numeric"))[[1]] + } +runIDs_done <- if (length(times) > 1) sort(times[-1]) else NULL + #timing: output for overall timing information timerfile2 <- "Timing_Simulation.csv" -if(file.exists(temp <- file.path(dir.out, timerfile2))) try(file.remove(temp), silent=TRUE) -write.table(t(c("", "Time_s", "Number")), file=temp, append=TRUE, sep=",", dec=".", col.names=FALSE, row.names=FALSE) +ftemp <- file.path(dir.out, timerfile2) +cat(",Time_s,Number", file = ftemp, sep = "\n") #concatenate file keeps track of sql files inserted into data concatFile <- "sqlFilesInserted.csv" @@ -158,108 +154,127 @@ if(!parallel_runs) { #if(print.debug) trace(what=circular:::SdCircularRad, tracer=quote({print(x); print(sys.calls()[[6]]); print(paste(rbar, circsd))}), at=4) #------prepare output -aon.help <- matrix(data=output_aggregates, ncol=2, nrow=length(output_aggregates)/2, byrow=TRUE) -aon <- data.frame(t(as.numeric(aon.help[,-1]))) -names(aon) <- aon.help[,1] +temp <- matrix(data=output_aggregates, ncol=2, nrow=length(output_aggregates)/2, byrow=TRUE) +aon <- lapply(temp[, 2], function(x) as.logical(as.numeric(x))) +names(aon) <- temp[, 1] #------import data if(!be.quiet) print(paste("SWSF reads input data: started at", t1 <- Sys.time())) +# Read data from files do_check_include <- FALSE if (usePreProcessedInput && file.exists(file.path(dir.in, datafile.SWRWinputs_preprocessed))) { load(file = file.path(dir.in, datafile.SWRWinputs_preprocessed), envir = .GlobalEnv) # This however annihilates all objects in .GlobalEnv with the same names ! } else { - # Read data from files - SWRunInformation <- tryCatch(read.csv(file.path(dir.in, datafile.SWRunInformation), as.is=TRUE),error=function(e) { print("datafile.SWRunInformation: Bad Path"); print(e)}) - stopifnot(sapply(c("Label", "site_id", "WeatherFolder", "X_WGS84", "Y_WGS84", "ELEV_m", "Include_YN"), - function(x) x %in% colnames(SWRunInformation)), # required columns - all(SWRunInformation$site_id == seq_len(nrow(SWRunInformation))), # consecutive site_id - !grepl("[[:space:]]", SWRunInformation$Label), # no space-characters in label - !grepl("[[:space:]]", SWRunInformation$WeatherFolder) # no space-characters in weather-data names - ) - include_YN <- SWRunInformation$Include_YN - labels <- SWRunInformation$Label - - sw_input_soillayers <- tryCatch(read.csv(file.path(dir.in, datafile.soillayers)),error=function(e) { print("datafile.soillayers: Bad Path"); print(e)}) - - sw_input_treatments_use <- tryCatch(read.csv(temp <- file.path(dir.in, datafile.treatments), nrows=1),error=function(e) { print("datafile.treatments: Bad Path"); print(e)}) - sw_input_treatments <- read.csv(temp, skip=1, as.is=TRUE) - colnames(sw_input_treatments) <- colnames(sw_input_treatments_use) - stopifnot( - !grepl("[[:space:]]", sw_input_treatments$LookupWeatherFolder) # no space-characters in weather-data names - ) - - sw_input_experimentals_use <- tryCatch(read.csv(temp <- file.path(dir.in, datafile.Experimentals), nrows=1),error=function(e) { print("datafile.Experimentals: Bad Path"); print(e)}) - sw_input_experimentals <- read.csv(temp, skip=1, as.is=TRUE) - colnames(sw_input_experimentals) <- colnames(sw_input_experimentals_use) - create_experimentals <- names(sw_input_experimentals_use[-1][which(sw_input_experimentals_use[-1] > 0 & is.finite(as.numeric(sw_input_experimentals_use[-1])))]) - stopifnot( - !grepl("[[:space:]]", sw_input_experimentals$LookupWeatherFolder) # no space-characters in weather-data names - ) - - #update treatment specifications based on experimental design - sw_input_treatments_use_combined <- ifelse(sw_input_treatments_use[-1] == 1 | names(sw_input_treatments_use[-1]) %in% create_experimentals, 1, 0) - temp<-which(!(create_experimentals %in% names(sw_input_treatments_use[-1]))) - if(length(temp) != 0) sw_input_treatments_use_combined <- cbind(sw_input_treatments_use_combined, matrix(data=1,nrow=1,ncol=length(temp),dimnames=list(NA, c(create_experimentals[temp])))) - create_treatments <- names(sw_input_treatments_use_combined[,which(sw_input_treatments_use_combined > 0 & is.finite(as.numeric(sw_input_treatments_use_combined)))]) - - if(dim(SWRunInformation)[2] == 1) stop("SWRunInformation might be tab separated instead of comma.") - if(dim(sw_input_soillayers)[2] == 1) stop("SoilLayers might be tab separated instead of comma.") - if(dim(sw_input_treatments_use)[2] == 1) stop("Treatments might be tab separated instead of comma.") - if(dim(sw_input_experimentals_use)[2] == 1) stop("Experimentals might be tab separated instead of comma.") - - sw_input_cloud_use <- sw_input_cloud <- list() - sw_input_prod_use <- sw_input_prod <- list() - sw_input_site_use <- sw_input_site <- list() - sw_input_soils_use <- sw_input_soils <- list() - sw_input_weather_use <- sw_input_weather <- list() - sw_input_climscen_use <- sw_input_climscen <- list() - sw_input_climscen_values_use <- sw_input_climscen_values <- list() - tr_files <- tr_prod <- tr_site <- tr_soil <- tr_weather <- tr_cloud <- list() - tr_input_climPPT <- tr_input_climTemp <- tr_input_shiftedPPT <- tr_input_EvapCoeff <- tr_input_TranspCoeff_Code <- tr_input_TranspCoeff <- tr_input_TranspRegions <- tr_input_SnowD <- tr_VegetationComposition <- list() - - sw_input_cloud_use <- tryCatch(read.csv(temp <- file.path(dir.sw.dat, datafile.cloud), nrows=1),error=function(e) { print("datafile.cloud: Bad Path"); print(e)}) - sw_input_cloud <- read.csv(temp, skip=1) - colnames(sw_input_cloud) <- colnames(sw_input_cloud_use) - - sw_input_prod_use <- tryCatch(read.csv(temp <- file.path(dir.sw.dat, datafile.prod), nrows=1),error=function(e) { print("datafile.prod: Bad Path"); print(e)}) - sw_input_prod <- read.csv(temp, skip=1) - colnames(sw_input_prod) <- colnames(sw_input_prod_use) - sw_input_prod_use[-1] <- ifelse(sw_input_prod_use[-1] == 1 | names(sw_input_prod_use[-1]) %in% create_experimentals, 1, 0) #update specifications based on experimental design - - sw_input_site_use <- tryCatch(read.csv(temp <- file.path(dir.sw.dat, datafile.siteparam), nrows=1),error=function(e) { print("datafile.siteparam: Bad Path"); print(e)}) - sw_input_site <- read.csv(temp, skip=1) - colnames(sw_input_site) <- colnames(sw_input_site_use) - sw_input_site_use[-1] <- ifelse(sw_input_site_use[-1] == 1 | names(sw_input_site_use[-1]) %in% create_experimentals, 1, 0) #update specifications based on experimental design - - sw_input_soils_use <- tryCatch(read.csv(temp <- file.path(dir.sw.dat, datafile.soils), nrows=1),error=function(e) { print("datafile.soils: Bad Path"); print(e)}) - sw_input_soils <- read.csv(temp, skip=1) - colnames(sw_input_soils) <- colnames(sw_input_soils_use) - sw_input_soils_use[-1] <- ifelse(sw_input_soils_use[-1] == 1 | names(sw_input_soils_use[-1]) %in% create_experimentals, 1, 0) #update specifications based on experimental design - - sw_input_weather_use <- tryCatch(read.csv(temp <- file.path(dir.sw.dat, datafile.weathersetup), nrows=1),error=function(e) { print("datafile.weathersetup: Bad Path"); print(e)}) - sw_input_weather <- read.csv(temp, skip=1) - colnames(sw_input_weather) <- colnames(sw_input_weather_use) - - sw_input_climscen_use <- tryCatch(read.csv(temp <- file.path(dir.sw.dat, datafile.climatescenarios), nrows=1),error=function(e) { print("datafile.climatescenarios: Bad Path"); print(e)}) - sw_input_climscen <- read.csv(temp, skip=1) - colnames(sw_input_climscen) <- colnames(sw_input_climscen_use) - - sw_input_climscen_values_use <- tryCatch(read.csv(temp <- file.path(dir.sw.dat, datafile.climatescenarios_values), nrows=1),error=function(e) { print("datafile.climatescenarios_values: Bad Path"); print(e)}) - sw_input_climscen_values <- read.csv(temp, skip=1) - colnames(sw_input_climscen_values) <- colnames(sw_input_climscen_values_use) - - if(dim(sw_input_cloud_use)[2] == 1) stop("Cloud datafile might be tab separated instead of comma.") - if(dim(sw_input_prod_use)[2] == 1) stop("Prod datafile might be tab separated instead of comma.") - if(dim(sw_input_site_use)[2] == 1) stop("Site datafile might be tab separated instead of comma.") - if(dim(sw_input_soils_use)[2] == 1) stop("Soils datafile might be tab separated instead of comma.") - if(dim(sw_input_weather_use)[2] == 1) stop("Weather datafile might be tab separated instead of comma.") - if(dim(sw_input_climscen_use)[2] == 1) stop("Climate Use datafile datafile might be tab separated instead of comma.") - if(dim(sw_input_climscen_values_use)[2] == 1) stop("Climate Values datafile datafile might be tab separated instead of comma.") - - #Create a list of possible treatment files with data. + SWRunInformation <- tryCatch(swsf_read_csv(file.path(dir.in, datafile.SWRunInformation)), + error = print) + stopifnot(sapply(c("Label", "site_id", "WeatherFolder", "X_WGS84", "Y_WGS84", "ELEV_m", "Include_YN"), + function(x) x %in% names(SWRunInformation)), # required columns + all(SWRunInformation$site_id == seq_len(nrow(SWRunInformation))), # consecutive site_id + !grepl("[[:space:]]", SWRunInformation$Label), # no space-characters in label + !grepl("[[:space:]]", SWRunInformation$WeatherFolder) # no space-characters in weather-data names + ) + include_YN <- SWRunInformation$Include_YN + labels <- SWRunInformation$Label + + sw_input_soillayers <- tryCatch(swsf_read_csv(file.path(dir.in, datafile.soillayers)), + error = print) + + temp <- tryCatch(swsf_read_inputfile(file.path(dir.in, datafile.treatments)), + error = print) + sw_input_treatments_use <- temp[["use"]] + sw_input_treatments <- temp[["data"]] + stopifnot( + !grepl("[[:space:]]", sw_input_treatments$LookupWeatherFolder) # no space-characters in weather-data names + ) + + temp <- tryCatch(swsf_read_inputfile(file.path(dir.in, datafile.Experimentals)), + error = print) + sw_input_experimentals_use <- temp[["use"]] + sw_input_experimentals <- temp[["data"]] + create_experimentals <- names(sw_input_experimentals_use[sw_input_experimentals_use]) + stopifnot( + !grepl("[[:space:]]", sw_input_experimentals$LookupWeatherFolder) # no space-characters in weather-data names + ) + + #update treatment specifications based on experimental design + sw_input_treatments_use_combined <- sw_input_treatments_use | + names(sw_input_treatments_use) %in% create_experimentals + temp <- which(!(create_experimentals %in% names(sw_input_treatments_use))) + if (length(temp) > 0) { + sw_input_treatments_use_combined <- cbind(sw_input_treatments_use_combined, + matrix(1, nrow = 1, ncol = length(temp), dimnames = list(NA, c(create_experimentals[temp])))) + } + create_treatments <- names(sw_input_treatments_use_combined)[sw_input_treatments_use_combined] + + if (dim(SWRunInformation)[2] == 1) + stop("SWRunInformation might be tab separated instead of comma.") + if (dim(sw_input_soillayers)[2] == 1) + stop("SoilLayers might be tab separated instead of comma.") + if (dim(sw_input_treatments)[2] == 1) + stop("Treatments might be tab separated instead of comma.") + if (dim(sw_input_experimentals)[2] == 1) + stop("Experimentals might be tab separated instead of comma.") + + temp <- tryCatch(swsf_read_inputfile(file.path(dir.sw.dat, datafile.cloud)), + error = print) + sw_input_cloud_use <- temp[["use"]] + sw_input_cloud <- temp[["data"]] + + temp <- tryCatch(swsf_read_inputfile(file.path(dir.sw.dat, datafile.prod)), + error = print) + sw_input_prod <- temp[["data"]] + sw_input_prod_use <- temp[["use"]] + sw_input_prod_use <- sw_input_prod_use | names(sw_input_prod_use) %in% create_experimentals #update specifications based on experimental design + + temp <- tryCatch(swsf_read_inputfile(file.path(dir.sw.dat, datafile.siteparam)), + error = print) + sw_input_site <- temp[["data"]] + sw_input_site_use <- temp[["use"]] + sw_input_site_use <- sw_input_site_use | names(sw_input_site_use) %in% create_experimentals #update specifications based on experimental design + + temp <- tryCatch(swsf_read_inputfile(file.path(dir.sw.dat, datafile.soils)), + error = print) + sw_input_soils <- temp[["data"]] + sw_input_soils_use <- temp[["use"]] + sw_input_soils_use <- sw_input_soils_use | names(sw_input_soils_use) %in% create_experimentals #update specifications based on experimental design + + temp <- tryCatch(swsf_read_inputfile(file.path(dir.sw.dat, datafile.weathersetup)), + error = print) + sw_input_weather_use <- temp[["use"]] + sw_input_weather <- temp[["data"]] + + temp <- tryCatch(swsf_read_inputfile(file.path(dir.sw.dat, datafile.climatescenarios)), + error = print) + sw_input_climscen_use <- temp[["use"]] + sw_input_climscen <- temp[["data"]] + + temp <- tryCatch(swsf_read_inputfile(file.path(dir.sw.dat, datafile.climatescenarios_values)), + error = print) + sw_input_climscen_values_use <- temp[["use"]] + sw_input_climscen_values <- temp[["data"]] + + if (dim(sw_input_cloud)[2] == 1) + stop("Cloud datafile might be tab separated instead of comma.") + if (dim(sw_input_prod)[2] == 1) + stop("Prod datafile might be tab separated instead of comma.") + if (dim(sw_input_site)[2] == 1) + stop("Site datafile might be tab separated instead of comma.") + if (dim(sw_input_soils)[2] == 1) + stop("Soils datafile might be tab separated instead of comma.") + if (dim(sw_input_weather)[2] == 1) + stop("Weather datafile might be tab separated instead of comma.") + if (dim(sw_input_climscen)[2] == 1) + stop("Climate Use datafile datafile might be tab separated instead of comma.") + if (dim(sw_input_climscen_values)[2] == 1) + stop("Climate Values datafile datafile might be tab separated instead of comma.") + + #Create a list of possible treatment files with data. + tr_files <- tr_prod <- tr_site <- tr_soil <- tr_weather <- tr_cloud <- list() + tr_input_climPPT <- tr_input_climTemp <- tr_input_shiftedPPT <- tr_input_EvapCoeff <- tr_input_TranspCoeff_Code <- tr_input_TranspCoeff <- tr_input_TranspRegions <- tr_input_SnowD <- tr_VegetationComposition <- list() + if(any(create_treatments=="sw")) print("SW treatment is not used because library Rsoilwat only uses one version of soilwat. Sorry") if(any(create_treatments=="filesin")) { @@ -287,11 +302,15 @@ if (usePreProcessedInput && file.exists(file.path(dir.in, datafile.SWRWinputs_pr tr_cloud[basename(temp)] <-unlist(lapply(temp,FUN=function(x) return(swReadLines(swClear(new("swCloud")),x)))) } - if(any(create_treatments == "LookupClimatePPTScenarios")) tr_input_climPPT <- read.csv( file.path(dir.sw.in.tr, "LookupClimatePPTScenarios", trfile.LookupClimatePPTScenarios)) - if(any(create_treatments == "LookupClimateTempScenarios")) tr_input_climTemp <- read.csv( file.path(dir.sw.in.tr, "LookupClimateTempScenarios", trfile.LookupClimateTempScenarios)) - if(any(create_treatments == "LookupShiftedPPTScenarios")) tr_input_shiftedPPT <- read.csv( file.path(dir.sw.in.tr, "LookupShiftedPPTScenarios", trfile.LookupShiftedPPTScenarios), row.names=1) - if(any(create_treatments == "LookupEvapCoeffFromTable")) tr_input_EvapCoeff <- read.csv( file.path(dir.sw.in.tr, "LookupEvapCoeffFromTable", trfile.LookupEvapCoeffFromTable), row.names=1) - if(any(create_treatments == "LookupTranspCoeffFromTable_Grass", create_treatments == "LookupTranspCoeffFromTable_Shrub", create_treatments == "LookupTranspCoeffFromTable_Tree", create_treatments == "LookupTranspCoeffFromTable_Forb", create_treatments == "AdjRootProfile")){ + if (any(create_treatments == "LookupClimatePPTScenarios")) + tr_input_climPPT <- swsf_read_csv(file.path(dir.sw.in.tr, "LookupClimatePPTScenarios", trfile.LookupClimatePPTScenarios)) + if (any(create_treatments == "LookupClimateTempScenarios")) + tr_input_climTemp <- swsf_read_csv(file.path(dir.sw.in.tr, "LookupClimateTempScenarios", trfile.LookupClimateTempScenarios)) + if (any(create_treatments == "LookupShiftedPPTScenarios")) + tr_input_shiftedPPT <- swsf_read_csv(file.path(dir.sw.in.tr, "LookupShiftedPPTScenarios", trfile.LookupShiftedPPTScenarios), row.names=1) + if (any(create_treatments == "LookupEvapCoeffFromTable")) + tr_input_EvapCoeff <- swsf_read_csv(file.path(dir.sw.in.tr, "LookupEvapCoeffFromTable", trfile.LookupEvapCoeffFromTable), row.names=1) + if (any(create_treatments == "LookupTranspCoeffFromTable_Grass", create_treatments == "LookupTranspCoeffFromTable_Shrub", create_treatments == "LookupTranspCoeffFromTable_Tree", create_treatments == "LookupTranspCoeffFromTable_Forb", create_treatments == "AdjRootProfile")) { tr_input_TranspCoeff_Code <- tryCatch(read.csv(temp <- file.path(dir.sw.in.tr, "LookupTranspCoeffFromTable", trfile.LookupTranspCoeffFromTable), nrows=2), error=function(e) { print("LookupTranspCoeffFromTable.csv: Bad Path"); print(e)}) tr_input_TranspCoeff_Code <- tr_input_TranspCoeff_Code[-2,] tr_input_TranspCoeff <- read.csv(temp, skip=2) @@ -323,7 +342,7 @@ if (usePreProcessedInput && file.exists(file.path(dir.in, datafile.SWRWinputs_pr } do_check_include <- TRUE - save(SWRunInformation, include_YN, labels, sw_input_soillayers, sw_input_treatments_use, sw_input_treatments, sw_input_experimentals_use, sw_input_experimentals, create_experimentals, sw_input_treatments_use_combined, create_treatments, sw_input_cloud_use, sw_input_cloud, sw_input_prod_use, sw_input_prod, sw_input_site_use, sw_input_site, sw_input_soils_use, sw_input_soils, sw_input_weather_use, sw_input_weather, sw_input_climscen_use, sw_input_climscen, sw_input_climscen_values_use, sw_input_climscen_values, tr_files, tr_prod, tr_site, tr_soil, tr_weather, tr_cloud, tr_input_climPPT, tr_input_climTemp, tr_input_shiftedPPT, tr_input_EvapCoeff, tr_input_TranspCoeff_Code, tr_input_TranspCoeff, tr_input_TranspRegions, tr_input_SnowD, tr_VegetationComposition, param.species_regeneration, no.species_regeneration, + save(SWRunInformation, include_YN, labels, sw_input_soillayers, sw_input_treatments_use, sw_input_treatments, sw_input_experimentals_use, sw_input_experimentals, create_experimentals, create_treatments, sw_input_cloud_use, sw_input_cloud, sw_input_prod_use, sw_input_prod, sw_input_site_use, sw_input_site, sw_input_soils_use, sw_input_soils, sw_input_weather_use, sw_input_weather, sw_input_climscen_use, sw_input_climscen, sw_input_climscen_values_use, sw_input_climscen_values, tr_files, tr_prod, tr_site, tr_soil, tr_weather, tr_cloud, tr_input_climPPT, tr_input_climTemp, tr_input_shiftedPPT, tr_input_EvapCoeff, tr_input_TranspCoeff_Code, tr_input_TranspCoeff, tr_input_TranspRegions, tr_input_SnowD, tr_VegetationComposition, param.species_regeneration, no.species_regeneration, file = file.path(dir.in, datafile.SWRWinputs_preprocessed), compress = FALSE) # No compression for fast access; RDS may be slightly faster, but would require loop over assign(, envir = .GlobalEnv) } @@ -414,8 +433,9 @@ if(any(simulation_timescales=="daily")){ #------------------------FLAGS FOR EXTERNAL DATA temp <- matrix(data=do.ExtractExternalDatasets, ncol=2, nrow=length(do.ExtractExternalDatasets)/2, byrow=TRUE) -exinfo <- data.frame(t(as.numeric(temp[,-1]))) -names(exinfo) <- temp[,1] +exinfo <- lapply(temp[, 2], function(x) as.logical(as.numeric(x))) +names(exinfo) <- temp[, 1] + exinfo$use_sim_spatial <- exinfo$ExtractSoilDataFromCONUSSOILFromSTATSGO_USA || @@ -553,9 +573,9 @@ if (extract_determine_database == "SWRunInformation" && "dailyweather_source" %i weather.digits <- 2 -lwf_cond1 <- sw_input_treatments_use$LookupWeatherFolder && sum(is.na(sw_input_treatments$LookupWeatherFolder[runIDs_sites])) == 0 -lwf_cond2 <- (sum(is.na(SWRunInformation$WeatherFolder[runIDs_sites])) == 0) && !any(as.logical(c(exinfo$GriddedDailyWeatherFromMaurer2002_NorthAmerica, exinfo$GriddedDailyWeatherFromDayMet_USA, exinfo$GriddedDailyWeatherFromNRCan_10km_Canada, exinfo$GriddedDailyWeatherFromNCEPCFSR_Global))) -lwf_cond3 <- sw_input_experimentals_use$LookupWeatherFolder && sum(is.na(sw_input_experimentals$LookupWeatherFolder)) == 0 +lwf_cond1 <- sw_input_treatments_use["LookupWeatherFolder"] && sum(is.na(sw_input_treatments$LookupWeatherFolder[runIDs_sites])) == 0 +lwf_cond2 <- (sum(is.na(SWRunInformation$WeatherFolder[runIDs_sites])) == 0) && !any(exinfo$GriddedDailyWeatherFromMaurer2002_NorthAmerica, exinfo$GriddedDailyWeatherFromDayMet_USA, exinfo$GriddedDailyWeatherFromNRCan_10km_Canada, exinfo$GriddedDailyWeatherFromNCEPCFSR_Global) +lwf_cond3 <- sw_input_experimentals_use["LookupWeatherFolder"] && sum(is.na(sw_input_experimentals$LookupWeatherFolder)) == 0 lwf_cond4 <- any(create_treatments == "LookupWeatherFolder") @@ -691,10 +711,10 @@ if(do_weather_source){ } dw_NCEPCFSR_Global <- function(sites_dailyweather_source) { - if(exinfo$GriddedDailyWeatherFromNCEPCFSR_Global && (simstartyr >= 1979 && endyr <= 2010)){ + if(exinfo$GriddedDailyWeatherFromNCEPCFSR_Global){ # Check which of the NCEPCFSR_Global weather data are available # - Grids domain: 0E to 359.688E and 89.761N to 89.761S - there <- simstartyr >= 1950 && endyr <= 2013 + there <- simstartyr >= 1979 && endyr <= 2010 if (any(there)) { there <- (SWRunInformation[runIDs_sites, "X_WGS84"] >= 0 - 180 & SWRunInformation[runIDs_sites, "X_WGS84"] <= 360 - 180) & (SWRunInformation[runIDs_sites, "Y_WGS84"] >= -89.761 & SWRunInformation[runIDs_sites, "Y_WGS84"] <= 89.761) if (any(there)) { @@ -761,7 +781,7 @@ if(getCurrentWeatherDataFromDatabase){ if(!(createAndPopulateWeatherDatabase || file.exists(dbWeatherDataFile))) stop("Create or use existing Weather database with Scenario data inside.") } else { - if(!any(create_treatments == "LookupWeatherFolder", as.logical(exinfo$GriddedDailyWeatherFromMaurer2002_NorthAmerica), as.logical(exinfo$GriddedDailyWeatherFromDayMet_NorthAmerica))) + if(!any(create_treatments == "LookupWeatherFolder", exinfo$GriddedDailyWeatherFromMaurer2002_NorthAmerica, exinfo$GriddedDailyWeatherFromDayMet_NorthAmerica)) stop("Daily weather data must be provided through 'LookupWeatherFolder', 'Maurer2002_NorthAmerica', or 'DayMet_NorthAmerica' since no weather database is used") } @@ -912,9 +932,8 @@ if (any(actions == "create")) { for (pc in do_prior_lookup) { if (any(create_treatments == pc$flag)) { #lookup values per category for each simulation run and copy values to datafile - if (!(sw_input_experimentals_use[pc$flag]) || - (as.logical(sw_input_experimentals_use[pc$flag]) && - length(unique(sw_input_experimentals[, pc$flag])) == 1L)) { + temp <- sw_input_experimentals_use[pc$flag] + if (!temp || (temp && (length(unique(sw_input_experimentals[, pc$flag])) == 1L))) { # Lookup prior to do_OneSite() only if option is off in sw_input_experimentals or constant if (continueAfterAbort) { @@ -922,7 +941,7 @@ if (any(actions == "create")) { sw_input_use <- get(pc$sw_input_use) icols <- grep(pc$pattern, names(sw_input_use)) - icols <- icols[sw_input_use[icols] > 0] + icols <- icols[sw_input_use[icols]] temp <- get(pc$sw_input)[, icols, drop = FALSE] if (all(!apply(is.na(temp), 2, all))) { @@ -935,7 +954,7 @@ if (any(actions == "create")) { if (!be.quiet) print(paste(Sys.time(), ": performing", shQuote(pc$flag))) - trtype <- if (as.logical(sw_input_experimentals_use[pc$flag])) { + trtype <- if (sw_input_experimentals_use[pc$flag]) { unique(sw_input_experimentals[, pc$flag]) } else { sw_input_treatments[, pc$flag] @@ -963,7 +982,8 @@ if (any(actions == "create")) { assign(pc$sw_input, tempdat$sw_input, envir = .GlobalEnv) #write data to datafile - write.csv(rbind(tempdat$sw_input_use, tempdat$sw_input), file = pc$datafile, row.names = FALSE) + write.csv(reconstitute_inputfile(tempdat$sw_input_use, tempdat$sw_input), + file = pc$datafile, row.names = FALSE) unlink(file.path(dir.in, datafile.SWRWinputs_preprocessed)) } @@ -985,14 +1005,14 @@ if (any(actions == "create")) { #------flags temp <- matrix(data=do.PriorCalculations, ncol=2, nrow=length(do.PriorCalculations)/2, byrow=TRUE) -pcalcs <- data.frame(t(as.numeric(temp[,-1]))) -names(pcalcs) <- temp[,1] +pcalcs <- lapply(temp[, 2], function(x) as.logical(as.numeric(x))) +names(pcalcs) <- temp[, 1] if (actionWithSoilWat) { - do.GetClimateMeans <- (sum(sw_input_climscen_values_use[-1]) > 0) || + do.GetClimateMeans <- any(sw_input_climscen_values_use) || pcalcs$EstimateConstantSoilTemperatureAtUpperAndLowerBoundaryAsMeanAnnualAirTemperature || - sw_input_site_use$SoilTempC_atLowerBoundary || - sw_input_site_use$SoilTempC_atUpperBoundary || + sw_input_site_use["SoilTempC_atLowerBoundary"] || + sw_input_site_use["SoilTempC_atUpperBoundary"] || pcalcs$EstimateInitialSoilTemperatureForEachSoilLayer || any(create_treatments == "PotentialNaturalVegetation_CompositionShrubsC3C4_Paruelo1996") || any(create_treatments == "AdjMonthlyBioMass_Temperature") || @@ -1000,7 +1020,7 @@ if (actionWithSoilWat) { any(create_treatments == "Vegetation_Biomass_ScalingSeason_AllGrowingORNongrowing") } -if (any(as.logical(pcalcs))) { +if (any(unlist(pcalcs))) { if (!be.quiet) print(paste("SWSF makes calculations prior to simulation runs: started at", t1 <- Sys.time())) @@ -1017,7 +1037,7 @@ if (any(as.logical(pcalcs))) { req_sl_ids <- paste0(requested_soil_layers, collapse = "x") # Available layers - ids_depth <- strsplit(names(sw_input_soils_use)[-1][as.logical(sw_input_soils_use)[-1]], "_", fixed = TRUE) + ids_depth <- strsplit(names(sw_input_soils_use)[sw_input_soils_use], "_", fixed = TRUE) stopifnot(length(ids_depth) > 0) var_layers <- unique(sapply(ids_depth, function(x) paste0(x[-length(x)], collapse = "_"))) ids_depth2 <- unique(sapply(ids_depth, function(x) x[length(x)])) @@ -1073,7 +1093,7 @@ if (any(as.logical(pcalcs))) { i.temp <- grep(var_layers[iv], names(sw_input_soils_use))[lyrs] sw_input_soils[runIDs_sites_ws[il_set], i.temp] <- round(sw_input_soils_data2[[iv]][, lyrs], if (var_layers[iv] %in% sl_vars_sub) 4L else 2L) - sw_input_soils_use[i.temp] <- 1L + sw_input_soils_use[i.temp] <- TRUE } sw_input_soillayers[runIDs_sites_ws[il_set], @@ -1085,8 +1105,8 @@ if (any(as.logical(pcalcs))) { #write data to datafile.soillayers write.csv(sw_input_soillayers, file = file.path(dir.in, datafile.soillayers), row.names = FALSE) #write data to datafile.soils - tempdat <- rbind(sw_input_soils_use, sw_input_soils) - write.csv(tempdat, file = file.path(dir.sw.dat, datafile.soils), row.names = FALSE) + write.csv(reconstitute_inputfile(sw_input_soils_use, sw_input_soils), + file = file.path(dir.sw.dat, datafile.soils), row.names = FALSE) unlink(file.path(dir.in, datafile.SWRWinputs_preprocessed)) print("'InterpolateSoilDatafileToRequestedSoilLayers': don't forget to adjust lookup tables with per-layer values if applicable for this project") @@ -1110,13 +1130,13 @@ if (any(as.logical(pcalcs))) { icol_bsE <- grep("EvapCoeff", names(sw_input_soils_use)) icol_sand <- grep("Sand_L", names(sw_input_soils_use)) icol_clay <- grep("Clay_L", names(sw_input_soils_use)) - use_layers <- which(sw_input_soils_use[icol_sand] == 1 & sw_input_soils_use[icol_clay] == 1) + use_layers <- which(sw_input_soils_use[icol_sand] & sw_input_soils_use[icol_clay]) stopifnot(length(use_layers) > 0) do_calc <- TRUE if (continueAfterAbort) { temp <- icol_bsE[use_layers] - icols <- temp[sw_input_soils_use[temp] > 0] + icols <- temp[sw_input_soils_use[temp]] if (length(icols) > 0L) { do_calc <- !all(rowSums(sw_input_soils[runIDs_sites, icols, drop = FALSE], na.rm = TRUE) > 0) } @@ -1169,18 +1189,18 @@ if (any(as.logical(pcalcs))) { icols_bsE_used <- icol_bsE[icol] icols_bse_notused <- icol_bsE[-icol] - sw_input_soils_use[icols_bsE_used] <- 1 + sw_input_soils_use[icols_bsE_used] <- TRUE sw_input_soils[runIDs_sites, icols_bsE_used] <- round(coeff_bs_evap[, icol], 4) - sw_input_soils_use[icols_bse_notused] <- 0 + sw_input_soils_use[icols_bse_notused] <- FALSE sw_input_soils[runIDs_sites, icols_bse_notused] <- 0 #write data to datafile.soils - tempdat <- rbind(sw_input_soils_use, sw_input_soils) - write.csv(tempdat, file = file.path(dir.sw.dat, datafile.soils), row.names = FALSE) + write.csv(reconstitute_inputfile(sw_input_soils_use, sw_input_soils), + file = file.path(dir.sw.dat, datafile.soils), row.names = FALSE) unlink(file.path(dir.in, datafile.SWRWinputs_preprocessed)) - rm(tempdat, icol, icols_bsE_used, icols_bse_notused, coeff_bs_evap, temp_coeff, + rm(icol, icols_bsE_used, icols_bse_notused, coeff_bs_evap, temp_coeff, lyrs_bs_evap, depth_bs_evap, temp_depth, ldepth_max_bs_evap, sand, clay, sand_mean, clay_mean, depth_min_bs_evap, layers_depth) } @@ -1191,45 +1211,20 @@ if (any(as.logical(pcalcs))) { print(paste(Sys.time(), "completed 'CalculateBareSoilEvaporationCoefficientsFromSoilTexture'")) } -# SoilWat >=v31 calculates field capacity and wilting point internally; they are no longer required as inputs and this option has become obsolete -# if(pcalcs$CalculateFieldCapacityANDWiltingPointFromSoilTexture){ -# #lookup soil texture data from 'datafile.soils' for those that the use flag of sand and clay is set, and calculate field capacity and wilting point -# ld <- 1:SoilLayer_MaxNo -# use.layers <- which(sw_input_soils_use[match(paste("Sand_L", ld, sep=""), colnames(sw_input_soils_use))] == 1) -# sand <- sw_input_soils[, match(paste("Sand_L", use.layers, sep=""), colnames(sw_input_soils_use))] -# clay <- sw_input_soils[, match(paste("Clay_L", use.layers, sep=""), colnames(sw_input_soils_use))] -# -# fieldc <- sapply(1:ncol(sand), FUN=function(i) SWPtoVWC(-0.033, sand[, i], clay[, i])) -# wiltp <- sapply(1:ncol(sand), FUN=function(i) SWPtoVWC(-1.5, sand[, i], clay[, i])) -# -# #add data to sw_input_cloud and set the use flags -# sw_input_soils_use[i.fieldc <- match(paste("FieldC_L", use.layers, sep=""), colnames(sw_input_soils_use))] <- 1 -# sw_input_soils_use[i.wiltp <- match(paste("WiltP_L", use.layers, sep=""), colnames(sw_input_soils_use))] <- 1 -# sw_input_soils[, i.fieldc] <- fieldc -# sw_input_soils[, i.wiltp] <- wiltp -# -# #write data to datafile.soils -# tempdat <- rbind(sw_input_soils_use, sw_input_soils) -# write.csv(tempdat, file=file.path(dir.sw.dat, datafile.soils), row.names=FALSE) -# unlink(file.path(dir.in, datafile.SWRWinputs_preprocessed)) -# -# rm(use.layers, sand, clay, fieldc, wiltp, tempdat, i.fieldc, i.wiltp) -# } - #------used during each simulation run: define functions here if(any(actions == "create") && pcalcs$EstimateConstantSoilTemperatureAtUpperAndLowerBoundaryAsMeanAnnualAirTemperature){ - sw_input_site_use$SoilTempC_atLowerBoundary <- 1 #set use flag - sw_input_site_use$SoilTempC_atUpperBoundary <- 1 + sw_input_site_use["SoilTempC_atLowerBoundary"] <- TRUE #set use flag + sw_input_site_use["SoilTempC_atUpperBoundary"] <- TRUE #call function 'SiteClimate' in each SoilWat-run } if(any(actions == "create") && pcalcs$EstimateInitialSoilTemperatureForEachSoilLayer){ #set use flags - ld <- 1:SoilLayer_MaxNo - use.layers <- which(sw_input_soils_use[match(paste("Sand_L", ld, sep=""), colnames(sw_input_soils_use))] == 1) - index.soilTemp <- match(paste("SoilTemp_L", ld, sep=""), colnames(sw_input_soils_use))[use.layers] + ld <- seq_len(SoilLayer_MaxNo) + use.layers <- which(sw_input_soils_use[paste0("Sand_L", ld)]) + index.soilTemp <- paste0("SoilTemp_L", ld)[use.layers] soilTemp <- sw_input_soils[runIDs_sites, index.soilTemp, drop = FALSE] - sw_input_soils_use[index.soilTemp] <- 1 + sw_input_soils_use[index.soilTemp] <- TRUE } if(!be.quiet) print(paste("SWSF makes calculations prior to simulation runs: ended after", round(difftime(Sys.time(), t1, units="secs"), 2), "s")) @@ -1241,15 +1236,15 @@ if (any(actions == "map_input") && length(map_vars) > 0) { if (!be.quiet) print(paste("SWSF generates maps of input variables for quality control: started at", t1 <- Sys.time())) dir.create(dir.inmap <- file.path(dir.out, "Input_maps"), showWarnings = FALSE) - input_avail <- list(SWRunInformation = list(cols = colnames(SWRunInformation), use = rep(TRUE, ncol(SWRunInformation))), - sw_input_soillayers = list(cols = colnames(sw_input_soillayers), use = rep(TRUE, ncol(sw_input_soillayers))), - sw_input_cloud = list(cols = colnames(sw_input_cloud), use = as.logical(sw_input_cloud_use)), - sw_input_prod = list(cols = colnames(sw_input_prod), use = as.logical(sw_input_prod_use)), - sw_input_site = list(cols = colnames(sw_input_site), use = as.logical(sw_input_site_use)), - sw_input_soils = list(cols = colnames(sw_input_soils), use = as.logical(sw_input_soils_use)), - sw_input_weather = list(cols = colnames(sw_input_weather), use = as.logical(sw_input_weather_use)), - sw_input_climscen = list(cols = colnames(sw_input_climscen), use = as.logical(sw_input_climscen_use)), - sw_input_climscen_values = list(cols = colnames(sw_input_climscen_values), use = as.logical(sw_input_climscen_use)) + input_avail <- list(SWRunInformation = list(cols = names(SWRunInformation), use = rep(TRUE, ncol(SWRunInformation))), + sw_input_soillayers = list(cols = names(sw_input_soillayers), use = rep(TRUE, ncol(sw_input_soillayers))), + sw_input_cloud = list(cols = names(sw_input_cloud), use = sw_input_cloud_use), + sw_input_prod = list(cols = names(sw_input_prod), use = sw_input_prod_use), + sw_input_site = list(cols = names(sw_input_site), use = sw_input_site_use), + sw_input_soils = list(cols = names(sw_input_soils), use = sw_input_soils_use), + sw_input_weather = list(cols = names(sw_input_weather), use = sw_input_weather_use), + sw_input_climscen = list(cols = names(sw_input_climscen), use = sw_input_climscen_use), + sw_input_climscen_values = list(cols = names(sw_input_climscen_values), use = sw_input_climscen_use) ) for (iv in seq_along(map_vars)) { @@ -1441,7 +1436,7 @@ do_OneSite <- function(i_sim, i_labels, i_SWRunInformation, i_sw_input_soillayer } } else {# needs to be read from soilsin file if(tasks$create == -1L) stop("This currently doesn't work") #TODO make it work low PR - layers_depth <- swSoils_Layers(tr_soil[[soilsin]])[,1] + layers_depth <- swSoils_Layers(tr_soil[[soilsin]])[, 1] d <- length(layers_depth) soildepth <- max(layers_depth) } @@ -1507,8 +1502,8 @@ do_OneSite <- function(i_sim, i_labels, i_SWRunInformation, i_sw_input_soillayer #------1. Step: Information for this SoilWat-run from prepared SoilWat-run stored in dir.sw.in #Make a local copy of the swInput object do not want to destroy orignal - swRunScenariosData<-list(scenario_No) - swRunScenariosData[[1]]<-swDataFromFiles + swRunScenariosData <- list() + swRunScenariosData[[1]] <- swDataFromFiles #get folder and file names outsetupin <- swOutSetupIn @@ -1601,9 +1596,9 @@ do_OneSite <- function(i_sim, i_labels, i_SWRunInformation, i_sw_input_soillayer ) for (pc in do_lookup) { - if (as.logical(sw_input_experimentals_use[pc$flag]) && !done_prior[pc$flag]) { - if (any(is.na(i_sw_input_treatments[pc$flag])) || - !all(unique(i_sw_input_treatments[pc$flag]) %in% rownames(pc$tr_input))) { + if (sw_input_experimentals_use[pc$flag] && !done_prior[pc$flag]) { + if (any(is.na(i_sw_input_treatments[[pc$flag]])) || + !all(unique(i_sw_input_treatments[[pc$flag]]) %in% rownames(pc$tr_input))) { print(paste("ERROR:", shQuote(pc$flag), "column in expirementals cannot have any NAs or name is not in tr_input table.")) tasks$create <- 0L @@ -1632,109 +1627,51 @@ do_OneSite <- function(i_sim, i_labels, i_SWRunInformation, i_sw_input_soillayer } #Treatment chunks - if(print.debug) print("Start of LookupTranspCoeff") - if(print.debug) print(".........grass") - if(any(create_treatments == "LookupTranspCoeffFromTable_Grass")){ - if(temp<-is.na(i_sw_input_treatments$LookupTranspCoeffFromTable_Grass)) print("LookupTranspCoeffFromTable_Grass for this run cannot be NA.") - if(temp1<-!all(i_sw_input_treatments$LookupTranspCoeffFromTable_Grass %in% colnames(tr_input_TranspCoeff))) print("LookupTranspCoeffFromTable_Grass name for this run are not in tr_input_TranspCoeff table column names.") - if(temp || temp1) { - tasks$create <- 0L - } else { - trco <- TranspCoeffByVegType( - tr_input_code = tr_input_TranspCoeff_Code, tr_input_coeff = tr_input_TranspCoeff, - soillayer_no = d, - trco_type = i_sw_input_treatments[1, "LookupTranspCoeffFromTable_Grass"], - layers_depth = layers_depth, - adjustType = "positive") - - if(!any(is.na(trco)) || sum(trco,na.rm=T) > 0){#trco does not have NA and sum is greater than 0. - #set the use flags - i.temp <- grepl(pattern=paste("Grass", "_TranspCoeff", sep=""), x=names(sw_input_soils_use)) - sw_input_soils_use[i.temp][1:length(trco)] <- rep(1, times=length(trco)) - #add data to sw_input_soils - i_sw_input_soils[i.temp][1:length(trco)] <- trco - } else { - print("The function TranspCoeffByVegType returned NA or does not sum to greater than 0 for this run for type grass.") - tasks$create <- 0L - } - } - } - if(print.debug) print(".........shrub") - if(any(create_treatments == "LookupTranspCoeffFromTable_Shrub")){ - if(temp<-is.na(i_sw_input_treatments$LookupTranspCoeffFromTable_Shrub)) print("LookupTranspCoeffFromTable_Shrub for this run cannot be NA.") - if(temp1<-!all(i_sw_input_treatments$LookupTranspCoeffFromTable_Shrub %in% colnames(tr_input_TranspCoeff))) print("LookupTranspCoeffFromTable_Shrub name for this run are not in tr_input_TranspCoeff table column names.") - if(temp || temp1) { - tasks$create <- 0L - } else { - trco <- TranspCoeffByVegType( - tr_input_code = tr_input_TranspCoeff_Code, tr_input_coeff = tr_input_TranspCoeff, - soillayer_no = d, - trco_type = i_sw_input_treatments[1, "LookupTranspCoeffFromTable_Shrub"], - layers_depth = layers_depth, - adjustType = "inverse") - - #set the use flags - if(!any(is.na(trco)) || sum(trco,na.rm=T) > 0){ - i.temp <- grepl(pattern=paste("Shrub", "_TranspCoeff", sep=""), x=names(sw_input_soils_use)) - sw_input_soils_use[i.temp][1:length(trco)] <- rep(1, times=length(trco)) - #add data to sw_input_soils - i_sw_input_soils[i.temp][1:length(trco)] <- trco - } else { - print("The function TranspCoeffByVegType returned NA or does not sum to greater than 0 for this run for type shrub.") - tasks$create <- 0L - } - } - } - if(print.debug) print(".........tree") - if(any(create_treatments == "LookupTranspCoeffFromTable_Tree")){ - if(temp<-is.na(i_sw_input_treatments$LookupTranspCoeffFromTable_Tree)) print("LookupTranspCoeffFromTable_Tree for this run cannot be NA.") - if(temp1<-!all(i_sw_input_treatments$LookupTranspCoeffFromTable_Tree %in% colnames(tr_input_TranspCoeff))) print("LookupTranspCoeffFromTable_Tree name for this run are not in tr_input_TranspCoeff table column names.") - if(temp || temp1) { - tasks$create <- 0L - } else { - trco <- TranspCoeffByVegType( - tr_input_code = tr_input_TranspCoeff_Code, tr_input_coeff = tr_input_TranspCoeff, - soillayer_no = d, - trco_type = i_sw_input_treatments[1, "LookupTranspCoeffFromTable_Tree"], - layers_depth = layers_depth, - adjustType = "inverse") - - if(!any(is.na(trco)) || sum(trco,na.rm=T) > 0){ - i.temp <- grepl(pattern=paste("Tree", "_TranspCoeff", sep=""), x=names(sw_input_soils_use)) - sw_input_soils_use[i.temp][1:length(trco)] <- rep(1, times=length(trco)) - #add data to sw_input_soils - i_sw_input_soils[i.temp][1:length(trco)] <- trco - } else { - print("The function TranspCoeffByVegType returned NA or does not sum to greater than 0 for this run for type Tree.") - tasks$create <- 0L - } - } - } - if(print.debug) print(".........forb") - if(any(create_treatments == "LookupTranspCoeffFromTable_Forb")){ - if(temp<-is.na(i_sw_input_treatments$LookupTranspCoeffFromTable_Forb)) print("LookupTranspCoeffFromTable_Forb for this run cannot be NA.") - if(temp1<-!all(i_sw_input_treatments$LookupTranspCoeffFromTable_Forb %in% colnames(tr_input_TranspCoeff))) print("LookupTranspCoeffFromTable_Forb name for this run are not in tr_input_TranspCoeff table column names.") - if(temp || temp1) { - tasks$create <- 0L - } else { - trco <- TranspCoeffByVegType( - tr_input_code = tr_input_TranspCoeff_Code, tr_input_coeff = tr_input_TranspCoeff, - soillayer_no = d, - trco_type = i_sw_input_treatments[1, "LookupTranspCoeffFromTable_Forb"], - layers_depth = layers_depth, - adjustType = "inverse") - - if(!any(is.na(trco)) || sum(trco,na.rm=T) > 0){ - i.temp <- grepl(pattern=paste("Forb", "_TranspCoeff", sep=""), x=names(sw_input_soils_use)) - sw_input_soils_use[i.temp][1:length(trco)] <- rep(1, times=length(trco)) - #add data to sw_input_soils - i_sw_input_soils[i.temp][1:length(trco)] <- trco - } else { - print("The function TranspCoeffByVegType returned NA or does not sum to greater than 0 for this run for type Forb.") - tasks$create <- 0L - } - } - } + if (print.debug) + print("Start of LookupTranspCoeff") + + do_vegs <- list( + veg = c("Grass", "Shrub", "Tree", "Forb"), + flag = c("LookupTranspCoeffFromTable_Grass", "LookupTranspCoeffFromTable_Shrub", + "LookupTranspCoeffFromTable_Tree", "LookupTranspCoeffFromTable_Forb"), + adjustType = c("positive", "inverse", "inverse", "inverse")) + + for (k in seq_along(do_vegs[["veg"]])) { + if (print.debug) + print(paste0(".........", do_vegs[["veg"]][k])) + + if (any(create_treatments == do_vegs[["flag"]][k])) { + temp <- is.na(i_sw_input_treatments[1, do_vegs[["flag"]][k]]) + temp1 <- !all(i_sw_input_treatments[1, do_vegs[["flag"]][k]] %in% colnames(tr_input_TranspCoeff)) + if (temp || temp1) { + if (temp) + print(paste(do_vegs[["flag"]][k], "for this run cannot be NA.")) + if (temp1) + print(paste(do_vegs[["flag"]][k], "name for this run are not in tr_input_TranspCoeff table column names.")) + tasks$create <- 0L + } else { + trco <- TranspCoeffByVegType( + tr_input_code = tr_input_TranspCoeff_Code, tr_input_coeff = tr_input_TranspCoeff, + soillayer_no = d, + trco_type = i_sw_input_treatments[1, do_vegs[["flag"]][k]], + layers_depth = layers_depth, + adjustType = do_vegs[["adjustType"]][k]) + + if (!any(is.na(trco)) || sum(trco, na.rm = TRUE) > 0) {#trco does not have NA and sum is greater than 0. + #set the use flags + i.temp <- grep(paste0(do_vegs[["veg"]][k], "_TranspCoeff"), names(sw_input_soils_use)) + sw_input_soils_use[i.temp[seq_along(trco)]] <- TRUE + if (length(i.temp) > length(trco)) + sw_input_soils_use[i.temp[(length(trco) + 1):length(i.temp)]] <- FALSE + #add data to sw_input_soils + i_sw_input_soils[i.temp[seq_along(trco)]] <- trco + } else { + print(paste("The function TranspCoeffByVegType returned NA or does not sum to greater than 0 for this run for type", do_vegs[["adjustType"]][k])) + tasks$create <- 0L + } + } + } + } #the monthly ppt-shifts are extracted, but written to the weathersetup input file only at the end of the create section 'copy and make climate scenarios from datafiles', because they are multiplied with any climate change factors ppt_scShift <- rep(1, times=12) @@ -1767,23 +1704,23 @@ do_OneSite <- function(i_sim, i_labels, i_SWRunInformation, i_sw_input_soillayer if(do.wind <- datafile.windspeedAtHeightAboveGround != SoilWat.windspeedAtHeightAboveGround) wind <- adjust.WindspeedHeight(uz=wind, height=datafile.windspeedAtHeightAboveGround) - if(sum(sw_input_cloud_use[-1]) > 0 | do.wind){ + if (any(sw_input_cloud_use) | do.wind) { #sky cover - if(sum(sw_input_cloud_use[grepl(pattern="SkyC", x=names(sw_input_cloud_use))]) > 0) { + if (any(sw_input_cloud_use[grepl("SkyC", names(sw_input_cloud_use))])) { sky <- with(i_sw_input_cloud, data.frame(SkyC_1, SkyC_2, SkyC_3, SkyC_4, SkyC_5, SkyC_6, SkyC_7, SkyC_8, SkyC_9, SkyC_10, SkyC_11, SkyC_12)) swCloud_SkyCover(swRunScenariosData[[1]]) <- round(as.double(sky), 0) } #wind speed - if(sum(sw_input_cloud_use[grepl(pattern="wind", x=names(sw_input_cloud_use))]) > 0 | do.wind) { + if (any(sw_input_cloud_use[grepl("wind", names(sw_input_cloud_use))]) | do.wind) { swCloud_WindSpeed(swRunScenariosData[[1]]) <- round(as.double(wind), 2) } #relative humidity - if(sum(sw_input_cloud_use[grepl(pattern="RH", x=names(sw_input_cloud_use))]) > 0) { + if (any(sw_input_cloud_use[grepl("RH", names(sw_input_cloud_use))])) { rh <- with(i_sw_input_cloud, data.frame(RH_1, RH_2, RH_3, RH_4, RH_5, RH_6, RH_7, RH_8, RH_9, RH_10, RH_11, RH_12)) swCloud_Humidity(swRunScenariosData[[1]]) <- round(as.double(rh), 0) } #snow density - if(sum(sw_input_cloud_use[grepl(pattern="snowd", x=names(sw_input_cloud_use))]) > 0) { + if (any(sw_input_cloud_use[grepl("snowd", names(sw_input_cloud_use))])) { snowd <- with(i_sw_input_cloud, data.frame(snowd_1, snowd_2, snowd_3, snowd_4, snowd_5, snowd_6, snowd_7, snowd_8, snowd_9, snowd_10, snowd_11, snowd_12)) if(i_SWRunInformation$Y_WGS84 < 0 && i_sw_input_cloud$SnowD_Hemisphere == "N" || i_SWRunInformation$Y_WGS84 > 0 && i_sw_input_cloud$SnowD_Hemisphere == "S"){ #adjust for hemisphere only if location and data are opposite snowd <- c(snowd[7:12], snowd[1:6]) @@ -1794,107 +1731,96 @@ do_OneSite <- function(i_sim, i_labels, i_SWRunInformation, i_sw_input_soillayer #add vegetation information from datafile to prodin if(print.debug) print("Start of prodin") - if(sum(sw_input_prod_use[-1]) > 0){ + + if (any(sw_input_prod_use)) { #composition - if(sum(use_comp <- unlist(sw_input_prod_use[grepl(pattern="Composition", x=names(sw_input_prod_use))])) > 0) { - comp.datfile <- with(i_sw_input_prod, data.frame(Composition_GrassFraction, Composition_ShrubFraction, Composition_TreeFraction, Composition_ForbFraction, Composition_BareGround)) - comp.datfile[is.na(comp.datfile)]<-0 - swRunScenariosData[[1]]@prod@Composition[1:5]<- as.numeric(comp.datfile[1:length(use_comp)]) + use_it <- sw_input_prod_use[grepl("Composition", names(sw_input_prod_use))] + if (any(use_it)) { + temp <- with(i_sw_input_prod, cbind( + Composition_GrassFraction, Composition_ShrubFraction, Composition_TreeFraction, + Composition_ForbFraction, Composition_BareGround)) + temp[is.finite(temp) | !use_it] <- 0 # if one is is requested, then put others to 0 + swRunScenariosData[[1]]@prod@Composition <- temp } #albedo - if(sum(use_albedo <- unlist(sw_input_prod_use[grepl(pattern="Albedo", x=names(sw_input_prod_use))])) > 0) { - albedo.datfile <- with(i_sw_input_prod, data.frame(Grass_Albedo, Shrub_Albedo, Tree_Albedo, Forb_Albedo, BareGround_Albedo)) - swProd_Albedo(swRunScenariosData[[1]])[use_albedo] <- albedo.datfile[use_albedo] + use_it <- sw_input_prod_use[grepl("Albedo", names(sw_input_prod_use))] + if (any(use_it)) { + temp <- with(i_sw_input_prod, cbind(Grass_Albedo, Shrub_Albedo, Tree_Albedo, + Forb_Albedo, BareGround_Albedo)) + temp[is.finite(temp)] <- 0 + swProd_Albedo(swRunScenariosData[[1]])[use_it] <- temp[use_it] } #constant canopy height - if(sum(use_height <- unlist(sw_input_prod_use[grepl(pattern="CanopyHeight_Constant", x=names(sw_input_prod_use))])) > 0) { - height.datfile <- with(i_sw_input_prod, data.frame(Grass_CanopyHeight_Constant_cm, Shrub_CanopyHeight_Constant_cm, Tree_CanopyHeight_Constant_cm,Forb_CanopyHeight_Constant_cm)) - height.datfile[is.na(height.datfile)]<-0 - swRunScenariosData[[1]]@prod@CanopyHeight[5,] <- as.numeric(height.datfile[1:length(use_height)]) + use_it <- sw_input_prod_use[grepl("CanopyHeight_Constant", names(sw_input_prod_use))] + if (any(use_it)) { + temp <- with(i_sw_input_prod, cbind(Grass_CanopyHeight_Constant_cm, + Shrub_CanopyHeight_Constant_cm, Tree_CanopyHeight_Constant_cm, + Forb_CanopyHeight_Constant_cm)) + temp[is.finite(temp)] <- 0 + swRunScenariosData[[1]]@prod@CanopyHeight[5, use_it] <- height.datfile[use_it] } #flag for hydraulic redistribution - if(sum(use_HD <- unlist(sw_input_prod_use[grepl(pattern="HydRed", x=names(sw_input_prod_use))])) > 0) { - HD.datfile <- with(i_sw_input_prod, data.frame(Grass_HydRed_OnOff, Shrub_HydRed_OnOff, Tree_HydRed_OnOff, Forb_HydRed_OnOff)) - swProd_HydrRedstro_use(swRunScenariosData[[1]])[use_HD] <- as.logical(HD.datfile[use_HD]) - } - #biomass components TODO Check This - biomassComponents <- function(FunctGroup){ - if( sum(litt <- sw_input_prod_use[grepl(pattern=paste(FunctGroup, "_Litter", sep=""), x=names(sw_input_prod_use))]) + - sum(biom <- sw_input_prod_use[grepl(pattern=paste(FunctGroup, "_Biomass", sep=""), x=names(sw_input_prod_use))]) + - sum(live <- sw_input_prod_use[grepl(pattern=paste(FunctGroup, "_FractionLive", sep=""), x=names(sw_input_prod_use))]) + - sum(laiconv <- sw_input_prod_use[grepl(pattern=paste(FunctGroup, "_LAIconv", sep=""), x=names(sw_input_prod_use))]) > 0) { - - for (m in st_mo){ - input_in_veg<-with(i_sw_input_prod,eval(parse(text=paste0("swRunScenariosData[[1]]@prod@MonthlyProductionValues_",tolower(FunctGroup),"[m,]")))) - mo.dat <- with(i_sw_input_prod, c( ifelse(litt[m], eval(parse(text=paste(FunctGroup, "_Litter_m", m, sep=""))),input_in_veg[1]), - ifelse(biom[m], eval(parse(text=paste(FunctGroup, "_Biomass_m", m, sep=""))), input_in_veg[2]), - ifelse(live[m], eval(parse(text=paste(FunctGroup, "_FractionLive_m", m, sep=""))), input_in_veg[3]), - ifelse(laiconv[m], eval(parse(text=paste(FunctGroup, "_LAIconv_m", m, sep=""))), input_in_veg[4]))) - if(FunctGroup=="Grass") - swRunScenariosData[[1]]@prod@MonthlyProductionValues_grass[m,] <- mo.dat - if(FunctGroup=="Shrub") - swRunScenariosData[[1]]@prod@MonthlyProductionValues_shrub[m,] <- mo.dat - if(FunctGroup=="Tree") - swRunScenariosData[[1]]@prod@MonthlyProductionValues_tree[m,] <- mo.dat - if(FunctGroup=="Forb") - swRunScenariosData[[1]]@prod@MonthlyProductionValues_forb[m,] <- mo.dat - } - } - if(FunctGroup=="Grass") - return(swProd_MonProd_grass(swRunScenariosData[[1]])) - if(FunctGroup=="Shrub") - return(swProd_MonProd_shrub(swRunScenariosData[[1]])) - if(FunctGroup=="Tree") - return(swProd_MonProd_tree(swRunScenariosData[[1]])) - if(FunctGroup=="Forb") - return(swProd_MonProd_forb(swRunScenariosData[[1]])) + use_it <- sw_input_prod_use[grepl("HydRed", names(sw_input_prod_use))] + if (any(use_it)) { + temp <- with(i_sw_input_prod, cbind(Grass_HydRed_OnOff, Shrub_HydRed_OnOff, + Tree_HydRed_OnOff, Forb_HydRed_OnOff)) + temp[is.finite(temp)] <- 0 + swProd_HydrRedstro_use(swRunScenariosData[[1]])[use_it] <- temp[use_it] } - swProd_MonProd_grass(swRunScenariosData[[1]]) <- biomassComponents(FunctGroup="Grass") - swProd_MonProd_shrub(swRunScenariosData[[1]]) <- biomassComponents(FunctGroup="Shrub") - swProd_MonProd_tree(swRunScenariosData[[1]]) <- biomassComponents(FunctGroup="Tree") - swProd_MonProd_forb(swRunScenariosData[[1]]) <- biomassComponents(FunctGroup="Forb") + + swProd_MonProd_grass(swRunScenariosData[[1]]) <- update_biomass( + funct_veg = "Grass", use = sw_input_prod_use, prod_input = i_sw_input_prod, + prod_default = swRunScenariosData[[1]]@prod) + swProd_MonProd_shrub(swRunScenariosData[[1]]) <- update_biomass( + funct_veg = "Shrub", use = sw_input_prod_use, prod_input = i_sw_input_prod, + prod_default = swRunScenariosData[[1]]@prod) + swProd_MonProd_tree(swRunScenariosData[[1]]) <- update_biomass( + funct_veg = "Tree", use = sw_input_prod_use, prod_input = i_sw_input_prod, + prod_default = swRunScenariosData[[1]]@prod) + swProd_MonProd_forb(swRunScenariosData[[1]]) <- update_biomass( + funct_veg = "Forb", use = sw_input_prod_use, prod_input = i_sw_input_prod, + prod_default = swRunScenariosData[[1]]@prod) } #Moved adjust to southern Hemi #add site information to siteparamin if(print.debug) print("Start of siteparamin") - if(sum(sw_input_site_use[-1]) > 0){ - site_swc_use <- as.logical(c(sw_input_site_use$SWC_min,sw_input_site_use$SWC_init,sw_input_site_use$SWC_wet)) - if(any(site_swc_use)){ - swSite_SWClimits(swRunScenariosData[[1]])[temp] <- c(i_sw_input_site$SWC_min,i_sw_input_site$SWC_init,i_sw_input_site$SWC_wet)[temp] - } - site_modelflag_use <- as.logical(c(sw_input_site_use$SWC_YearlyReset,sw_input_site_use$SWC_Deepdrain)) - if(any(site_modelflag_use)){ - swSite_ModelFlags(swRunScenariosData[[1]])[site_modelflag_use] <- c(i_sw_input_site$SWC_YearlyReset,i_sw_input_site$SWC_Deepdrain)[site_modelflag_use] - } - site_modelcoef_use <- as.logical(c(sw_input_site_use$PET_multiplier,sw_input_site_use$RunoffPercent_fromPondedWater)) - if(any(site_modelcoef_use)){ - swSite_ModelCoefficients(swRunScenariosData[[1]])[site_modelcoef_use] <- c(i_sw_input_site$PET_multiplier, i_sw_input_site$RunoffPercent_fromPondedWater)[site_modelcoef_use] - } - #replace if in treatment file - if(sw_input_site_use$Param_UnsaturatedPercolation){ - swSite_DrainageCoefficient(swRunScenariosData[[1]]) <- i_sw_input_site$Param_UnsaturatedPercolation - } - if(sw_input_site_use$Latitude) { - swSite_IntrinsicSiteParams(swRunScenariosData[[1]])[1] <- i_sw_input_site$Latitude - } - if(sw_input_site_use$Altitude) { - swSite_IntrinsicSiteParams(swRunScenariosData[[1]])[2] <- i_sw_input_site$Altitude - } - if(sw_input_site_use$Slope){ - swSite_IntrinsicSiteParams(swRunScenariosData[[1]])[3] <- i_sw_input_site$Slope - } - if(sw_input_site_use$Aspect){ - swSite_IntrinsicSiteParams(swRunScenariosData[[1]])[4] <- i_sw_input_site$Aspect - } - #Moved sw_input_site_use$SoilTempC_atLowerBoundary - if(sw_input_site_use$SoilTemp_Flag){ - swSite_SoilTemperatureFlag(swRunScenariosData[[1]]) <- i_sw_input_site$SoilTemp_Flag - } - rm(site_swc_use,site_modelflag_use,site_modelcoef_use) - } - swSite_IntrinsicSiteParams(swRunScenariosData[[1]])[1] <- i_SWRunInformation$Y_WGS84 * pi / 180 - if(is.finite(i_SWRunInformation$ELEV_m)) swSite_IntrinsicSiteParams(swRunScenariosData[[1]])[2] <- i_SWRunInformation$ELEV_m + if (any(sw_input_site_use)) { + flags <- c("SWC_min", "SWC_init", "SWC_wet") + site_use <- sw_input_site_use[flags] + if (any(site_use)) + swSite_SWClimits(swRunScenariosData[[1]])[site_use] <- + as.numeric(i_sw_input_site[flags][site_use]) + + flags <- c("SWC_YearlyReset", "SWC_Deepdrain") + site_use <- sw_input_site_use[flags] + if (any(site_use)) + swSite_ModelFlags(swRunScenariosData[[1]])[site_use] <- + as.numeric(i_sw_input_site[flags][site_use]) + + flags <- c("PET_multiplier", "RunoffPercent_fromPondedWater") + site_use <- sw_input_site_use[flags] + if (any(site_use)) + swSite_ModelCoefficients(swRunScenariosData[[1]])[site_use] <- + as.numeric(i_sw_input_site[flags][site_use]) + + if (sw_input_site_use["Param_UnsaturatedPercolation"]) { + swSite_DrainageCoefficient(swRunScenariosData[[1]]) <- i_sw_input_site$Param_UnsaturatedPercolation + } + + flags <- c("Latitude", "Altitude", "Slope", "Aspect") + site_use <- sw_input_site_use[flags] + if (any(site_use)) + swSite_IntrinsicSiteParams(swRunScenariosData[[1]])[site_use] <- + as.numeric(i_sw_input_site[flags][site_use]) + + if (sw_input_site_use["SoilTemp_Flag"]) { + swSite_SoilTemperatureFlag(swRunScenariosData[[1]]) <- i_sw_input_site$SoilTemp_Flag + } + } + swSite_IntrinsicSiteParams(swRunScenariosData[[1]])[1] <- i_SWRunInformation$Y_WGS84 * pi / 180 + if (is.finite(i_SWRunInformation$ELEV_m)) + swSite_IntrinsicSiteParams(swRunScenariosData[[1]])[2] <- i_SWRunInformation$ELEV_m #add soil information to soilsin if (print.debug) @@ -1907,13 +1833,13 @@ do_OneSite <- function(i_sim, i_labels, i_SWRunInformation, i_sw_input_soillayer dimnames(soil_swdat)[[2]] <- soil_cols done.Imperm_L1 <- FALSE - if (sw_input_soils_use$Imperm_L1 == 1 && any(create_treatments == "soilsin")) { + if (sw_input_soils_use["Imperm_L1"] && any(create_treatments == "soilsin")) { soil_swdat[1, "imperm"] <- i_sw_input_soils$Imperm_L1 done.Imperm_L1 <- TRUE } - use_transpregion <- as.numeric(sw_input_soils_use[paste0("TranspRegion_L", ld)]) - if (sum(sw_input_soils_use[-1]) + {if (done.Imperm_L1) -1 else 0} - sum(use_transpregion) > 0) { + use_transpregion <- sw_input_soils_use[paste0("TranspRegion_L", ld)] + if (sum(sw_input_soils_use) + {if (done.Imperm_L1) -1 else 0} - sum(use_transpregion) > 0) { # Calculate soil layer structure, because any(create_treatments=="soilsin") and soilsin may have a different soil layer structure than the datafiles temp <- as.numeric(na.omit(unlist(i_sw_input_soillayers[paste0("depth_L", seq_len(SoilLayer_MaxNo))]))) @@ -1950,9 +1876,9 @@ do_OneSite <- function(i_sim, i_labels, i_SWRunInformation, i_sw_input_soillayer icol <- icol[ld] if (length(temp) > 0) { - luse <- list(use = which(as.logical(sw_input_soils_use[icol])), + luse <- list(use = which(sw_input_soils_use[icol]), other = intersect( - which(!as.logical(sw_input_soils_use[icol])), + which(!sw_input_soils_use[icol]), seq_len(dim(soil_swdat)[1]))) for (k in 1:2) if (any(luse[[k]])) { temp <- if (k == 1L) { @@ -2064,11 +1990,11 @@ do_OneSite <- function(i_sim, i_labels, i_SWRunInformation, i_sw_input_soillayer } #add weather setup information to weatherin - if(sw_input_weather_use$SnowFlag) + if (sw_input_weather_use["SnowFlag"]) swWeather_UseSnow(swRunScenariosData[[1]]) <- as.logical(i_sw_input_weather$SnowFlag) - if(sw_input_weather_use$SnowDrift_Percent) + if (sw_input_weather_use["SnowDrift_Percent"]) swWeather_pct_SnowDrift(swRunScenariosData[[1]]) <- i_sw_input_weather$SnowDrift_Percent - if(sw_input_weather_use$RunOffOnPerSnowmelt_Percent) + if (sw_input_weather_use["RunOffOnPerSnowmelt_Percent"]) swWeather_pct_SnowRunoff(swRunScenariosData[[1]]) <- i_sw_input_weather$RunOffOnPerSnowmelt_Percent swWeather_FirstYearHistorical(swRunScenariosData[[1]]) <- simstartyr @@ -2144,22 +2070,22 @@ do_OneSite <- function(i_sim, i_labels, i_SWRunInformation, i_sw_input_soillayer SiteClimate_Ambient <- sw_SiteClimate_Ambient(weatherList=i_sw_weatherList[[1]], year.start=min(simTime$useyrs), year.end=max(simTime$useyrs), do.C4vars=do.C4vars, simTime2=simTime2) } } - if(!getScenarioWeatherDataFromDatabase) { - #get climate change information - use_pptValscen <- sw_input_climscen_values_use[, pptVal.colnames <- paste("PPTmm_m", st_mo, "_sc", formatC(sc-1, width=2, format="d", flag="0"), sep="")] - use_tempValMinScen <- sw_input_climscen_values_use[, tempValMin.colnames <- paste("TempC_min_m", st_mo, "_sc", formatC(sc-1, width=2, format="d", flag="0"), sep="")] - use_tempValMaxScen <- sw_input_climscen_values_use[, tempValMax.colnames <- paste("TempC_max_m", st_mo, "_sc", formatC(sc-1, width=2, format="d", flag="0"), sep="")] + if (!getScenarioWeatherDataFromDatabase) { + #get climate change information + cols_climscen_val <- lapply(c("PPTmm_m", "TempC_min_m", "TempC_max_m"), function(flag) + paste0(flag, st_mo, "_sc", formatC(sc - 1, width = 2, format = "d", flag = "0"))) + use_climscen_val <- any(unlist(sw_input_climscen_values_use[unlist(cols_climscen_val)])) - use_pptscen <- sw_input_climscen_use[, ppt.colnames <- paste("PPTfactor_m", st_mo, "_sc", formatC(sc-1, width=2, format="d", flag="0"), sep="")] - use_tempMinScen <- sw_input_climscen_use[, tempMin.colnames <- paste("deltaTempC_min_m", st_mo, "_sc", formatC(sc-1, width=2, format="d", flag="0"), sep="")] - use_tempMaxScen <- sw_input_climscen_use[, tempMax.colnames <- paste("deltaTempC_max_m", st_mo, "_sc", formatC(sc-1, width=2, format="d", flag="0"), sep="")] + cols_climscen_delta <- lapply(c("PPTfactor_m", "deltaTempC_min_m", "deltaTempC_max_m"), function(flag) + paste0(flag, st_mo, "_sc", formatC(sc - 1, width = 2, format = "d", flag = "0"))) + use_climscen_delta <- any(unlist(sw_input_climscen_use[unlist(cols_climscen_delta)])) - if( sum(use_pptValscen) + sum(use_tempValMinScen) + sum(use_tempValMaxScen) > 0){ + if (any(use_climscen_val)) { #convert climate change values to factors #read values from datafile - pptVal_sc <- unlist(i_sw_input_climscen_values[, pptVal.colnames]) - tVal_min_sc <- unlist(i_sw_input_climscen_values[, tempValMin.colnames]) - tVal_max_sc <- unlist(i_sw_input_climscen_values[, tempValMax.colnames]) + pptVal_sc <- unlist(i_sw_input_climscen_values[, cols_climscen_val[[1]]]) + tVal_min_sc <- unlist(i_sw_input_climscen_values[, cols_climscen_val[[2]]]) + tVal_max_sc <- unlist(i_sw_input_climscen_values[, cols_climscen_val[[3]]]) #calculate change factors ppt_sc <- pptVal_sc / (10 * SiteClimate_Ambient$meanMonthlyPPTcm) if(sum(abs(tVal_max_sc - tVal_min_sc)) > tol){ @@ -2168,22 +2094,21 @@ do_OneSite <- function(i_sim, i_labels, i_SWRunInformation, i_sw_input_soillayer } else { #no information for tmin, tmax by GCM -> tmin=tmax=tmean t_min_sc <- t_max_sc <- tVal_min_sc - SiteClimate_Ambient$meanMonthlyTempC } - } else if( sum(use_pptscen) + sum(use_tempMinScen) + sum(use_tempMaxScen) > 0){ + } else if (any(use_climscen_delta)) { #read climate change factors from datafile - ppt_sc <- unlist(i_sw_input_climscen[, ppt.colnames]) - t_min_sc <- unlist(i_sw_input_climscen[, tempMin.colnames]) - t_max_sc <- unlist(i_sw_input_climscen[, tempMax.colnames]) + ppt_sc <- unlist(i_sw_input_climscen[, cols_climscen_delta[[1]]]) + t_min_sc <- unlist(i_sw_input_climscen[, cols_climscen_delta[[2]]]) + t_max_sc <- unlist(i_sw_input_climscen[, cols_climscen_delta[[3]]]) } else { - ppt_sc <- rep(1, times=12) - t_min_sc <- rep(0, times=12) - t_max_sc <- rep(0, times=12) + ppt_sc <- rep(1, times = 12) + t_min_sc <- t_max_sc <- rep(0, times=12) } #guarantee that all entries are finite: this may not be the case for instance if any(meanMonthlyClimate$meanMonthlyPPTcm == 0) ppt_sc <- temp_ppt_sc <- ifelse(is.finite(ppt_sc), ppt_sc, 1) t_min_sc <- ifelse(is.finite(t_min_sc), t_min_sc, 0) t_max_sc <- ifelse(is.finite(t_max_sc), t_max_sc, 0) - if(sc > 1){ + if (sc > 1) { if(any(create_treatments=="ClimateScenario_Temp_PerturbationInMeanSeasonalityBothOrNone") && !grepl("Both", i_sw_input_treatments$ClimateScenario_Temp_PerturbationInMeanSeasonalityBothOrNone, ignore.case=T)){ if(grepl("Mean", i_sw_input_treatments$ClimateScenario_Temp_PerturbationInMeanSeasonalityBothOrNone, ignore.case=T)) { t_min_sc <- rep(mean(t_min_sc), times=12) @@ -2194,8 +2119,7 @@ do_OneSite <- function(i_sim, i_labels, i_SWRunInformation, i_sw_input_soillayer t_max_sc <- t_max_sc - mean(t_max_sc) } if(grepl("None", i_sw_input_treatments$ClimateScenario_Temp_PerturbationInMeanSeasonalityBothOrNone, ignore.case=T)) { - t_min_sc <- rep(0, times=12) - t_max_sc <- rep(0, times=12) + t_min_sc <- t_max_sc <- rep(0, times=12) } } if(any(create_treatments=="ClimateScenario_PPT_PerturbationInMeanSeasonalityBothOrNone") && !grepl("Both", i_sw_input_treatments$ClimateScenario_PPT_PerturbationInMeanSeasonalityBothOrNone, ignore.case=T)){ @@ -2206,12 +2130,13 @@ do_OneSite <- function(i_sim, i_labels, i_SWRunInformation, i_sw_input_soillayer } } - ppt_old <- (temp <- swWeather_MonScalingParams(swRunScenariosData[[sc]]))[,1] + temp <- swWeather_MonScalingParams(swRunScenariosData[[sc]]) + ppt_old <- temp[, 1] t_max_old <- temp[,2] t_min_old <- temp[,3] #write information into weatherin - if(sum(use_pptValscen) + sum(use_tempValMinScen) + sum(use_tempValMaxScen) + sum(use_pptscen) + sum(use_tempMinScen) + sum(use_tempMaxScen) > 0){ + if (any(use_climscen_val, use_climscen_delta)) { ppt_f <- ppt_sc t_min_f <- t_min_sc t_max_f <- t_max_sc @@ -2226,9 +2151,9 @@ do_OneSite <- function(i_sim, i_labels, i_SWRunInformation, i_sw_input_soillayer rownames(MonthlyScalingParams)<-c("January","February","March","April","May","June","July","August","September","October","November","December") swWeather_MonScalingParams(swRunScenariosData[[sc]])[, 1:3] <- MonthlyScalingParams - ClimatePerturbationsVals[sc,1:12] <- MonthlyScalingParams[,1] - ClimatePerturbationsVals[sc,13:24] <- MonthlyScalingParams[,2] - ClimatePerturbationsVals[sc,25:36] <- MonthlyScalingParams[,3] + ClimatePerturbationsVals[sc,1:12] <- MonthlyScalingParams[, 1] + ClimatePerturbationsVals[sc,13:24] <- MonthlyScalingParams[, 2] + ClimatePerturbationsVals[sc,25:36] <- MonthlyScalingParams[, 3] #Update climate data with climate scenario information if(do.GetClimateMeans){ @@ -2356,14 +2281,14 @@ do_OneSite <- function(i_sim, i_labels, i_SWRunInformation, i_sw_input_soillayer #out.temp <- data.frame(i_sim, i_label, soilTUpper, soilTlower) #write.csv(out.temp, file=paste(dir.out.temp, .Platform$file.sep, flag.icounter, "_", "SoilTempC_atLowerBoundary.csv", sep=""), quote=FALSE, row.names=FALSE) } - if(sw_input_site_use$SoilTempC_atUpperBoundary) { - soilTUpper <- ifelse(exists("soilTUpper"), soilTUpper, i_sw_input_site$SoilTempC_atUpperBoundary) + if (sw_input_site_use["SoilTempC_atUpperBoundary"]) { + soilTUpper <- if (exists("soilTUpper")) soilTUpper else i_sw_input_site$SoilTempC_atUpperBoundary } - if(sw_input_site_use$SoilTempC_atLowerBoundary){ - soilTlower <- ifelse(exists("soilTlower"), soilTlower, i_sw_input_site$SoilTempC_atLowerBoundary) + if (sw_input_site_use["SoilTempC_atLowerBoundary"]) { + soilTlower <- if (exists("soilTlower")) soilTlower else i_sw_input_site$SoilTempC_atLowerBoundary swSite_SoilTemperatureConsts(swRunScenariosData[[sc]])[8] <- soilTlower } - if(pcalcs$EstimateInitialSoilTemperatureForEachSoilLayer){ + if (pcalcs$EstimateInitialSoilTemperatureForEachSoilLayer) { init.soilTprofile <- EstimateInitialSoilTemperatureForEachSoilLayer(layers_depth=layers_depth, lower.Tdepth=as.numeric(swRunScenariosData[[sc]]@site@SoilTemperatureConstants[10]), soilTupper=soilTUpper, soilTlower=soilTlower) #lower.Tdepth needs to be adjusted if it changes in soilparam.in #temporaly save data #TODO get this working #out.temp <- data.frame(i_sim, i_label, t(c(init.soilTprofile, rep(NA, times=SoilLayer_MaxNo-length(init.soilTprofile))))) @@ -2371,8 +2296,8 @@ do_OneSite <- function(i_sim, i_labels, i_SWRunInformation, i_sw_input_soillayer } #adjust init soil temperatures to climatic conditions - if(any(use_soil_temp <- as.logical(sw_input_soils_use[paste("SoilTemp_L", ld, sep="")]))){ - + use_soil_temp <- sw_input_soils_use[paste0("SoilTemp_L", ld)] + if (any(use_soil_temp)) { temp <- 1:nrow(swSoils_Layers(swRunScenariosData[[sc]])) if(exists("init.soilTprofile")) { swSoils_Layers(swRunScenariosData[[sc]])[,12][use_soil_temp] <- init.soilTprofile @@ -3761,9 +3686,9 @@ do_OneSite <- function(i_sim, i_labels, i_SWRunInformation, i_sw_input_soillayer impermeability <- 0.9 #impermeable layer #Required soil layers - soildat <- swSoils_Layers(swRunScenariosData[[sc]])[, c("depth_cm", "sand_frac", "clay_frac", "impermeability_frac")] + soildat <- swSoils_Layers(swRunScenariosData[[sc]])[, c("depth_cm", "sand", "clay", "imperm")] #50cm soil depth or impermeable layer (whichever is shallower; Soil Survey Staff 2014: p.31) - imp_depth <- which(soildat[, "impermeability_frac"] >= impermeability) + imp_depth <- which(soildat[, "imperm"] >= impermeability) imp_depth <- min(imp_depth, max(soildat[, "depth_cm"])) #Interpret maximum soil depth as possible impermeable layer Fifty_depth <- min(50, imp_depth) @@ -3846,7 +3771,7 @@ do_OneSite <- function(i_sim, i_labels, i_SWRunInformation, i_sw_input_soillayer paste(layers_depth, collapse = ", "), "}")) swp_dy_nrsc <- get_SWPmatric_aggL(vwc_dy_nrsc, texture = texture, - sand = soildat[, "sand_frac"], clay = soildat[, "clay_frac"]) + sand = soildat[, "sand"], clay = soildat[, "clay"]) } else { swp_dy_nrsc <- swpmatric.dy.all @@ -5505,21 +5430,49 @@ do_OneSite <- function(i_sim, i_labels, i_SWRunInformation, i_sw_input_soillayer } #end if do aggregate - if(tasks$aggregate > 0L && (length(SQL) > 0 || length(SQLcurrent) > 0)){ - tasks$aggregate <- 2L - if(length(SQLcurrent) > 0) write(SQLcurrent, dbTempFileCurrent, append=TRUE) - if(length(SQL) > 0) write(SQL, dbTempFile, append=TRUE) - } + if (tasks$aggregate > 0L && (length(SQL) > 0 || length(SQLcurrent) > 0)) { + tasks$aggregate <- 2L + if (length(SQLcurrent) > 0) + cat(SQLcurrent, file = dbTempFileCurrent, append = TRUE, sep = "\n") + if (length(SQL) > 0) + cat(SQL, file = dbTempFile, append = TRUE, sep = "\n") + } - if(all(unlist(tasks) != 0)){ - #ETA estimation - dt <- difftime(Sys.time(), time.sys, units="secs") - times <- read.csv(file=file.path(dir.out, timerfile), header=FALSE, colClasses=c("NULL", "numeric")) - if(!be.quiet) print(paste(i_sim, ":", i_label, "done in", round(dt, 2), units(dt), ":", round(nrow(times)/runsN_total*100, 2), "% complete, ETA =", Sys.time()+ceiling((runsN_total-(nrow(times)-1))/workersN)*mean(unlist(c(times, dt)), na.rm=TRUE) )) - write.table(data.frame(i=i_sim,dt=round(dt, 2)), file=file.path(dir.out, timerfile), append=TRUE, sep=",", dec=".", col.names=FALSE,row.names=FALSE) - } else { - print(paste(i_sim, ":", i_label, " unsuccessful:", paste(names(tasks), "=", tasks, collapse=", "))) - } + if (all(unlist(tasks) != 0)) { + #ETA estimation + dt <- round(difftime(Sys.time(), time.sys, units = "secs"), 2) + ftemp <- file.path(dir.out, timerfile) + cat(paste(i_sim, dt, sep = ","), file = ftemp, append = TRUE, sep = "\n") + times <- swsf_read_csv(file = ftemp, header = FALSE, + nrows = runsN_total + 1, colClasses = c("NULL", "numeric"))[[1]] + + if (!be.quiet) { + n <- length(times) - 1 + temp <- paste0(i_sim, ": ", i_label, " done in ", dt, " ", units(dt), ": ", + round(n / runsN_total * 100, 2), "% complete") + + if (eta.estimate) { + deta <- round(ceiling((runsN_total - n) / workersN) * + sapply(list(mean, sd), function(f) f(times, na.rm = TRUE))) + pi95 <- deta[2] * sqrt(1 + 1 / n) * qt(0.975, n) # 95% prediction interval + pi95 <- if (is.na(pi95)) "NA" else if (pi95 > 3600) { + paste(round(pi95 / 3600), "h") + } else if (pi95 > 60) { + paste(round(pi95 / 60), "min") + } else { + paste(round(pi95), "s") + } + temp <- paste0(temp, ", ETA (mean ± 95%-PI) = ", + Sys.time() + deta[1], " ± ", pi95) + } + + print(temp) + } + + } else { + print(paste0(i_sim, ": ", i_label, " unsuccessful: ", + paste(names(tasks), "=", tasks, collapse = ", "))) + } on.exit() @@ -5565,62 +5518,7 @@ if(actionWithSoilWat && runsN_todo > 0){ #Used for weather from files filebasename <- basename(swFiles_WeatherPrefix(swDataFromFiles)) - #objects to export (sorted alphabetically) - list.export <- c("accountNSHemispheres_agg", "accountNSHemispheres_veg", "AdjMonthlyBioMass", - "adjust.soilDepth", "adjust.WindspeedHeight", "adjustLayersDepth", "aon", - "be.quiet", "bin.prcpfreeDurations", "bin.prcpSizes", "circ.mean", "circ.range", - "circ.sd", "climate.conditions", "cloudin", - "continueAfterAbort", "cor2", "counter.digitsN", "create_experimentals", - "create_filename_for_Maurer2002_NorthAmerica", "create_treatments", "cut0Inf", - "daily_lyr_agg", "daily_no", - "datafile.windspeedAtHeightAboveGround", "dbOverallColumns", - "dbWeatherDataFile", "debug.dump.objects", "DegreeDayBase", - "Depth_TopLayers", "dir.create2", "dir.ex.daymet", "dir.ex.maurer2002", - "dir.out.temp", "dir.out", "dir.prj", "dir.sw.in.tr", "dir.sw.runs", - "dirname.sw.runs.weather", "do_OneSite", "do.GetClimateMeans", - "done_prior", "endDoyAfterDuration", "endyr", "estabin", - "establishment.delay", "establishment.duration", "establishment.swp.surface", - "EstimateInitialSoilTemperatureForEachSoilLayer", "exec_c_prefix", - "ExpInput_Seperator", "expN", - "ExtractGriddedDailyWeatherFromDayMet_NorthAmerica_swWeather", - "ExtractGriddedDailyWeatherFromMaurer2002_NorthAmerica", - "filebasename.WeatherDataYear", "filebasename", "fill_empty", "finite01", - "germination.duration", "germination.swp.surface", "get_DayMet_cellID", - "get_DayMet_NorthAmerica", "get.LookupFromTable", "get.month", - "getCurrentWeatherDataFromDatabase", "getLayersWidth", - "getScenarioWeatherDataFromDatabase", "getStartYear", - "growing.season.threshold.tempC", "increment_soiltemperature_deltaX_cm", - "it_exp", "it_Pid", "it_site", "lmax", "local_weatherDirName", - "makeInputForExperimentalDesign", "max.duration", "name.OutputDB", - "no.species_regeneration", "ouput_aggregated_ts", "output_aggregate_daily", - "parallel_backend", "parallel_runs", "param.species_regeneration", "pcalcs", - "PotentialNaturalVegetation_CompositionShrubsC3C4_Paruelo1996", "print.debug", - "prodin", "runIDs_sites", "runsN_master", "runsN_sites", "runsN_todo", - "runsN_total", "saveRsoilwatInput", "saveRsoilwatOutput", "scenario_No", - "season.end", "season.start", "setAggSoilLayerForAggDailyResponses", - "setBottomLayer", "setDeepestTopLayer", "setLayerSequence", "setTopLayer", - "shrub.fraction.limit", "simstartyr", "simTime_ForEachUsedTimeUnit_North", - "simTime_ForEachUsedTimeUnit_South", "simTime", "simTiming_ForEachUsedTimeUnit", - "simTiming", "simulation_timescales", "siteparamin", "SoilLayer_MaxNo", - "soilsin", "SoilWat.windspeedAtHeightAboveGround", "st_mo", - "startDoyOfDuration", "startyr", "sw_aet", "sw_dailyC4_TempVar", "sw_deepdrain", - "sw_evapsurface", "sw_evsoil", "sw_hd", "sw_inf_soil", "sw_input_climscen_use", - "sw_input_climscen_values_use", "sw_input_cloud_use", - "sw_input_experimentals_use", "sw_input_experimentals", "sw_input_prod_use", - "sw_input_site_use", "sw_input_soils_use", "sw_input_weather_use", - "sw_interception", "sw_percolation", "sw_pet", "sw_precip", "sw_runoff", - "sw_SiteClimate_Ambient", "sw_snow", "sw_soiltemp", "sw_swcbulk", - "sw_swpmatric", "sw_temp", "sw_transp", "sw_vwcbulk", "sw_vwcmatric", - "sw.inputs", "sw.outputs", "swcsetupin", "swDataFromFiles", "swFilesIn", - "swOutSetupIn", "SWPcrit_MPa", "SWPtoVWC", "tempError", "timerfile", - "Tmax_crit_C", "Tmin_crit_C", "tol", "toln", "tr_cloud", "tr_files", - "tr_input_climPPT", "tr_input_climTemp", "tr_input_EvapCoeff", - "tr_input_shiftedPPT", "tr_input_SnowD", "tr_input_TranspCoeff_Code", - "tr_input_TranspCoeff", "tr_input_TranspRegions", "tr_prod", "tr_site", - "tr_soil", "tr_VegetationComposition", "tr_weather", - "transferExpDesignToInput", "TranspCoeffByVegType", "VWCtoSWP", "weatherin", - "mpi_work", "workersN", "yearsin") - # List of objects required by do_OneSite which are not in rSWSF + # List of objects to export which are required by do_OneSite and are not in rSWSF (sorted alphabetically) list.export <- c("accountNSHemispheres_agg", "accountNSHemispheres_veg", "adjust.soilDepth", "aon", "be.quiet", "bin.prcpfreeDurations", "bin.prcpSizes", "climate.conditions", "cloudin", "continueAfterAbort", "counter.digitsN", @@ -5631,7 +5529,7 @@ if(actionWithSoilWat && runsN_todo > 0){ "dir.prj", "dir.sw.in.tr", "dir.sw.runs", "dirname.sw.runs.weather", "do_OneSite", "do.GetClimateMeans", "done_prior", "endyr", "estabin", "establishment.delay", "establishment.duration", "establishment.swp.surface", - "exec_c_prefix", "ExpInput_Seperator", "expN", "filebasename", + "eta.estimate", "exec_c_prefix", "ExpInput_Seperator", "expN", "filebasename", "filebasename.WeatherDataYear", "germination.duration", "germination.swp.surface", "get.month", "getCurrentWeatherDataFromDatabase", "getScenarioWeatherDataFromDatabase", "growing.season.threshold.tempC", "increment_soiltemperature_deltaX_cm", @@ -5761,9 +5659,10 @@ tryCatch({ } else if (tag == 4L) { #The slave had a problem with Soilwat record Slave number and the Run number. print("Problem with run:", complete, "on slave:", save_id, "at", Sys.time()) - write.csv(data.frame(Slave = slave_id, Run = complete), - file = file.path(dir.out, "ProblemRuns.csv"), - append = TRUE, row.names = FALSE, col.names = TRUE) + ftemp <- file.path(dir.out, "ProblemRuns.csv") + if (!file.exists(ftemp)) + cat("Slave,Run", file = ftemp, sep = "\n") + cat(paste(slave_id, complete, sep = ","), file = ftemp, append = TRUE, sep = "\n") } }, interrupt=function(interrupt) { @@ -5965,7 +5864,8 @@ if (any(actions == "concatenate")) { # Clean up and report if (OK) { - write(file.path(dir.out.temp, theFileList[j]), file = file.path(dir.out.temp,concatFile), append = TRUE) + cat(file.path(dir.out.temp, theFileList[j]), + file = file.path(dir.out.temp,concatFile), append = TRUE, sep = "\n") if (deleteTmpSQLFiles) try(file.remove(file.path(dir.out.temp, theFileList[j])), silent = TRUE) } @@ -6339,10 +6239,14 @@ if(!be.quiet && do.ensembles) print(paste("SWSF calculates ensembles: ended afte #--------------------------------------------------------------------------------------------------# #------------------------OVERALL TIMING -delta.overall <- difftime(Sys.time(), t.overall, units="secs") -if(!be.quiet) print(paste("SWSF: ended after", round(delta.overall, 2), "s")) +delta.overall <- difftime(Sys.time(), t.overall, units = "secs") +if (!be.quiet) + print(paste("SWSF: ended after", round(delta.overall, 2), "s")) -write.timer <- function(label, time_sec="", number=""){ write.table(t(c(label, time_sec, number)), file=file.path(dir.out, timerfile2), append=TRUE, sep=",", dec=".", col.names=FALSE, row.names=FALSE) } +write.timer <- function(label, time_sec = "", number = "") { + cat(paste(label, time_sec, number, sep = ","), file = file.path(dir.out, timerfile2), + append = TRUE, sep = "\n") +} write.timer("Time_Total", time_sec=delta.overall) write.timer("Time_Check", time_sec=delta.check) diff --git a/2_SWSF_p5of5_Functions_v51.R b/2_SWSF_p5of5_Functions_v51.R index 3be5aad3..bcec91b9 100644 --- a/2_SWSF_p5of5_Functions_v51.R +++ b/2_SWSF_p5of5_Functions_v51.R @@ -4,8 +4,6 @@ slot <- methods::slot - - #------ Constants output_timescales_maxNo <- 4L SoilLayer_MaxNo <- 20L @@ -16,6 +14,42 @@ tol <- sqrt(.Machine$double.eps) toln <- sqrt(.Machine$double.neg.eps) #------ Funtions +swsf_read_csv <- compiler::cmpfun(function(file, ...) { + if (requireNamespace("iotools", quietly = TRUE)) { + # faster than utils::read.csv + temp <- try(iotools::read.csv.raw(file = file, ...), silent = TRUE) + if (inherits(temp, "try-error")) { + read.csv(file = file, ...) + } else { + names(temp) <- gsub("\"", "", names(temp)) + temp + } + + } else { + read.csv(file = file, ...) + } +}) + +swsf_read_inputfile <- compiler::cmpfun(function(file, header_rows = 1) { + sw_use <- tryCatch(swsf_read_csv(file, nrows = header_rows), + error = function(e) print(paste("Failed to read file:", shQuote(basename(file)), "with", e))) + sw <- swsf_read_csv(file, skip = header_rows) + names(sw) <- names(sw_use) + sw_use <- c(FALSE, as.logical(as.numeric(sw_use[, -1]))) + sw_use[is.na(sw_use)] <- FALSE + names(sw_use) <- names(sw) + + list(use = sw_use, data = sw) +}) + +reconstitute_inputfile <- compiler::cmpfun(function(sw_use, data) { + temp <- as.data.frame(matrix(as.integer(sw_use), nrow = 1L)) + colnames(temp) <- names(sw_use) + temp[1, 1] <- "UseInformationToCreateSoilWatRuns" + rbind(temp, data) +}) + + getStartYear <- compiler::cmpfun(function(simstartyr) simstartyr + 1) set_PRAGMAs <- compiler::cmpfun(function(con, settings) { @@ -1484,11 +1518,11 @@ get.LookupFromTable <- compiler::cmpfun(function(pattern, trtype, tr_input, sw_i # add data to datafiles and set the use flags sw_input[, seq1_icols_in] <- res[, icols_res] - sw_input_use[seq1_icols_in] <- 1L + sw_input_use[seq1_icols_in] <- TRUE if (length(seq2_icols_in) > 0) { sw_input[, seq2_icols_in] <- NA - sw_input_use[seq2_icols_in] <- 0L + sw_input_use[seq2_icols_in] <- FALSE } list(sw_input_use = sw_input_use, @@ -1506,7 +1540,7 @@ fill_empty <- compiler::cmpfun(function(data, pattern, fill, tol = tol) { iempty <- is.na(data$sw_input[, ic]) | abs(data$sw_input[, ic]) < tol if (any(iempty)) { data$sw_input[iempty, ic] <- fill - data$sw_input_use[icols[k, "sw_input_use"]] <- 1L + data$sw_input_use[icols[k, "sw_input_use"]] <- TRUE } } @@ -1601,7 +1635,7 @@ EstimateInitialSoilTemperatureForEachSoilLayer <- compiler::cmpfun(function(laye #--put information from experimental design into appropriate input variables; create_treatments and the _use files were already adjusted for the experimental design when files were read in/created transferExpDesignToInput <- compiler::cmpfun(function(x, i_exp, df_exp, df_exp_use) { - temp <- match(names(df_exp)[df_exp_use == 1], names(x), nomatch = 0) + temp <- match(names(df_exp)[df_exp_use], names(x), nomatch = 0) ctemp <- temp[!(temp == 0)] if (length(ctemp) > 0) { cexp <- match(names(x)[ctemp], names(df_exp), nomatch = 0) @@ -1892,7 +1926,7 @@ get_Runoff_yr <- compiler::cmpfun(function(sc, x, st) { }) cor2 <- compiler::cmpfun(function(y) { - res <- try(cor(y[,1], y[,2]), silent = TRUE) + res <- try(cor(y[, 1], y[, 2]), silent = TRUE) if (inherits(res, "try-error")) NA else res }) @@ -2545,6 +2579,24 @@ GriddedDailyWeatherFromNCEPCFSR_Global <- compiler::cmpfun(function(site_ids, da }) +update_biomass <- compiler::cmpfun(function(funct_veg = c("Grass", "Shrub", "Tree", "Forb"), use, prod_input, prod_default) { + funct_veg <- match.arg(funct_veg) + + comps <- c("_Litter", "_Biomass", "_FractionLive", "_LAIconv") + veg_ids = lapply(comps, function(x) + grep(paste0(funct_veg, x), names(use))) + veg_incl = lapply(vegs_ids, function(x) use[x]) + + temp <- slot(prod_default, paste0("MonthlyProductionValues_", tolower(funct_veg))) + if (any(unlist(veg_incl))) { + for (k in seq_along(comps)) if (any(veg_incl[[k]])) + temp[veg_incl[[k]], k] <- prod_input[, veg_ids[[k]][veg_incl[[k]]]] + } + + temp +}) + + ######################## #------ GISSM functions # Schlaepfer, D.R., Lauenroth, W.K. & Bradford, J.B. (2014). Modeling regeneration responses of big sagebrush (Artemisia tridentata) to abiotic conditions. Ecol Model, 286, 66-77. diff --git a/DESCRIPTION b/DESCRIPTION index 38697586..eb69bcd5 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: SWSF Title: SOILWAT Simulation Framework -Version: 1.6.1 -Date: 2016-09-23 +Version: 1.6.2 +Date: 2016-09-24 Authors@R: c(person("Daniel", "Schlaepfer", email = "daniel.schlaepfer@unibas.ch", role = c("aut", "cre")), person("Caitlin", "Andrews", role = "ctb", email = "not@available.com")) person("Ryan", "Murphy", role = "ctb", email = "not@available.com")) @@ -20,6 +20,7 @@ Suggests: rgdal (>= 1.1.10), rgeos (>= 0.3.19), raster (>= 2.5.8), sp (>= 1.2.3), ncdf4 (>= 1.15), maps (>= 3.1.0), RCurl (>= 1.95.4.8), fastmatch (>= 1.0.4), + iotools (>= 0.1-12), Hmisc (>= 3.17.4), qmap (>= 1.0.4), DaymetR (>= 0.1)