Initial commit to gitlab

R/SED_climate_from_dataframe_parallel.R

+#' @title Script to create the Landclim SED climate files from a climate data frame
+#' @description  This script takes a climate data frame in the form:
+#'               x (lon), y (lat), year, temperature (jan - dec), precipitiation (jan - dec)
+#'               and creates the necessary spatially explicit data files (SED) for landclim
+#' @param
+#' @param
+#' @return NULL
+#' @keywords CoFoLaMo, Dischma, spin-up, climate
+#' @author Austin Haffenden
+#' @example /don't run {}
+#'
+create_SED_file_parallel <- function(clim_file, latitude, ref_elev, file_prefix) {
+
+  years <- unique(clim_file$year)
+
+  library(foreach)
+  library(doParallel)
+
+  # Calculate the number of cores
+  #no_cores <- detectCores() -2
+  no_cores <- 2
+
+  # Initiate cluster
+  cl <- makeCluster(no_cores, type = "FORK", outfile = "")
+  registerDoParallel(cl)
+
+  foreach(i = 1:length(years)) %dopar% {
+    gc()
+    curr_file <- file.path(lc_clim_path, paste0(paste(file_prefix, i-1, sep = "_"), ".txt"))
+
+    con <- file(curr_file, open = "wt")
+    writeLines("#header	info#", con)
+    writeLines(paste(latitude, "#average latitude of the sites -	used	for	cacluating	drought	index	in	model	bugmann#"), con)
+    writeLines(paste(ref_elev, "#reference elevation#"), con)
+    writeLines("#	mean	monthly	temperature	?C	total	monthly	precipitation	mm	#", con)
+    writeLines("#	x y year jan	feb	marz	april	mai	june	july	aug	sept	oct	nov	dec	jan	feb	marz	april	mai	june	july	aug	sept	oct	nov	dec	#", con)
+
+    temp_clim <- subset(clim_file, year == years[i])
+
+    #     for(j in 1:nrow(climate)){
+    #
+    #       writeLines(paste(temp_clim[j, ], collapse = " "), con)
+    #     }
+    #
+    write(t(as.matrix(temp_clim)), file = con, ncolumns = 27, append = TRUE)
+    close(con)
+
+
+  }
+
+  stopCluster(cl)
+}

R/SED_climate_from_list.R

+#' @title Script to create and check the Landclim SED climate files from a
+#' climate list
+#' @description  This script takes a list of climate data frames with
+#'               each element in the form:
+#'               x (lon), y (lat), year, temperature (jan - dec),
+#'               precipitiation (jan - dec) and creates the necessary spatially
+#'               explicit data files (SED) for landclim
+#' @param
+#' @param
+#' @return NULL
+#' @keywords CoFoLaMo, Dischma, spin-up, climate
+#' @author Austin Haffenden
+#' @example /don't run {}
+#'
+SED_climate_from_list <- function(clim_list, 
+                                  latitude, 
+                                  ref_elev, 
+                                  file_prefix, 
+                                  clim_dir, 
+                                  start_file_number = 0) { # number you want the years to start at
+  
+  save_clim <- function(list_no, start_file_number){
+    
+    curr_clim <- clim_list[[list_no]]
+    #print(curr_clim[1,])
+    
+    #curr_year <- unique(curr_clim$year)
+    #lc_year <- curr_year - min_year + year_offset
+    lc_year <- (list_no-1) + start_file_number
+    # cat("list_no is: ", list_no, "\n")
+    # cat("start_file_number is: ", start_file_number, "\n")
+    # cat("lc_year is: ", lc_year, "\n")
+    # stop()
+
+    curr_file <- file.path(clim_dir, paste0(paste(file_prefix, lc_year, sep = "_"), ".txt"))
+    cat(curr_file, "\n")
+    con <- file(curr_file, open = "wt")
+
+    writeLines("#header	info#", con)
+    writeLines(paste(latitude, "#average latitude of the sites -	used	for	cacluating	drought	index	in	model	bugmann#"), con)
+    writeLines(paste(ref_elev, "#reference elevation#"), con)
+    writeLines("#	mean	monthly	temperature	?C	total	monthly	precipitation	mm	#", con)
+    writeLines("#	x y year jan	feb	marz	april	mai	june	july	aug	sept	oct	nov	dec	jan	feb	marz	april	mai	june	july	aug	sept	oct	nov	dec	#", con)
+
+    write.table(as.matrix(curr_clim),
+                #file = file.path(clim_dir, paste0(paste(file_prefix, i-1, sep = "_"), ".txt") ),
+                file = con,
+                row.names = FALSE, col.names = FALSE,
+                append = TRUE)
+    close(con)
+    
+  }
+  
+  # now just for test to call function
+  # years <- unique(clim_list$year)  
+  # min_year <- min(years)
+  
+  #doParallel::registerDoParallel(cores = 3)
+  #plyr::ddply(clim_list, ~ year, save_clim, min_year) 
+              #.parallel = TRUE)
+  
+  list_levels <- seq(1, length(clim_list), 1)
+  
+  # start_no <- start_file_number
+  # return(start_no)
+  lapply(list_levels, save_clim, start_file_number)
+  
+  return(1)
+}
+
+

R/add_to_init_file.R

+#' @title
+#' @description
+#' @param site_dir
+#' @param site
+#' @param curr_init_file
+#' @param xy
+#' @param site_r
+
+add_to_init_file <- function(site_init_dir, 
+                             site, 
+                             curr_init_file, 
+                             cell_no, 
+                             site_r)  {
+
+  site_init <- read.csv(file.path(site_init_dir, paste0(site, "_tree_init.csv")), header = T)
+
+  site_init["row"] <- raster::rowFromCell(site_r, cell_no)-1
+  site_init["col"] <- raster::colFromCell(site_r, cell_no)-1
+  site_init["cell"] <- cell_no -1
+
+  temp_init <- rbind(curr_init_file, site_init)
+
+
+  return(temp_init)
+
+}

R/check_SED_files.R

+#' @title function to check climate files for consistency with the original data
+#' @description currently set up for the dischma CoFoLaMo runs so needs to be
+#' tested on other data - maybe this can be removed with testthat
+check_SED_files <- function(clim_file_prefix,
+                            clim_dir,
+                            temp_RTS,
+                            prcp_RTS,
+                            spin_up) {
+
+  
+  sink(file.path(clim_dir, 
+                 paste0(clim_file_prefix, 
+                        "RTS_and_clim_file_comparison_log.txt")))
+  # There are 914 climate files 
+  
+  # pick 91 numbers at random between 0 and 913
+  file_nos <- floor(runif(91, min=0, max=913))
+  
+  for(i in 1:length(file_nos)) {
+  
+    # read in the first climate file
+    in_file <- read.csv(file.path(clim_dir, 
+                                  paste0(clim_file_prefix, "_", file_nos[i], ".txt")), 
+                        skip = 5, sep = " ", header = FALSE)
+    
+    cat("*******************************************************", "\n")
+    cat("Current file is: ", paste0(clim_file_prefix, "_", file_nos[i], ".txt"), "\n")
+  
+    names(in_file) <- c("x", "y", "year",
+                        "J", "F", "M", "A", "M", "J", "J", "A", "S", "O", "N", "D",
+                        "J", "F", "M", "A", "M", "J", "J", "A", "S", "O", "N", "D")
+    # get the year 
+    curr_year <- unique(in_file[,3])
+    
+    # if the year is < 1930 get the year it corresponds to in the spinup
+    if(curr_year<1930) {
+      
+      curr_year <- subset(spin_up, year==curr_year)$yearOfClimateDataSet
+    }
+  
+    # from that climate file select 10% of the rows at random 
+    # there are 3306 rows so 33 rows
+    curr_samp <- dplyr::sample_n(in_file, 33)
+    
+      # for the first row
+      for(j in 1:nrow(curr_samp)) {
+        
+        # get the cell number from the xy coords
+        cell_no <- rts::cellFromXY(temp_RTS, c(curr_samp[j,]$x, curr_samp[j,]$y))
+  
+        # extract the row temp and precip data from the climate file
+        tcor_fm_file <- as.vector(curr_samp[j, 4:15])
+        prcp_fm_file <- as.vector(curr_samp[j, 16:27])
+    
+        # extract the tcor from the RTS climate file
+        tcor_fm_RTS <- as.vector(extract(temp_RTS, cell_no, as.character(curr_year)))
+        
+        # proof of failure - keep commented out
+        #tcor_fm_RTS <- as.vector(extract(temp_RTS, 100, as.character(curr_year)))
+        
+        # extract the prcp from the RTS climate file
+        prcp_fm_RTS <- as.vector(extract(prcp_RTS, cell_no, as.character(curr_year)))
+        
+        # proof of failure - keep commented out
+        #prcp_fm_RTS <- as.vector(extract(prcp_RTS, 357, as.character(curr_year)))
+        
+        # conversion of prcp to mm's is done later in origanl code
+        prcp_fm_RTS <- prcp_fm_RTS * 10
+        
+        # data is rounded differently so make the comparison
+        # with tolerance to allow for that
+        if(isTRUE(all.equal(tcor_fm_file, tcor_fm_RTS, 
+                  tolerance = 0.000001, 
+                  check.attributes = FALSE))){
+          
+          cat("X: ", curr_samp[j,]$x,
+              "Y: ", curr_samp[j,]$y,
+              "Cell #: ", cell_no, 
+              "- RTS and climate file TCOR compared: OK", "\n")
+      
+        } else {
+          cat("* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *", "\n")
+          cat("* X: ", curr_samp[j,]$x,
+              "Y: ", curr_samp[j,]$y,
+              "Cell #: ", cell_no, 
+              "- RTS and climate file TCOR compared: FAILED *", "\n")
+          cat("* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *", "\n")
+          sink()
+          stop("FAILURE! See log file in climate files directory", "\n")
+        }
+        if(isTRUE(all.equal(prcp_fm_file, prcp_fm_RTS, 
+                  tolerance = 0.000001, 
+                  check.attributes = FALSE))){
+          
+          cat("X: ", curr_samp[j,]$x,
+              "Y: ", curr_samp[j,]$y,
+              "Cell #: ", cell_no, 
+              "- RTS and climate file PRCP compared: OK", "\n")
+          
+        } else {
+          cat("* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *", "\n")
+          cat("* X: ", curr_samp[j,]$x,
+              "Y: ", curr_samp[j,]$y,
+              "Cell #: ", cell_no, 
+              "- RTS and climate file PRCP compared: FAILED*", "\n")
+          cat("* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *", "\n")
+          sink()
+          stop("FAILURE! See log file in climate files directory", "\n")
+        }
+        cat("\n")
+    }
+  }
+  sink()
+}
\ No newline at end of file

R/create_IFN_species_pars.R

+# creates init file parameters for the species not represented by randomly selecting
+# values from a species of the same leaf type
+create_IFN_species_pars <- function(spp_names, 
+                                    base_species_df, 
+                                    curr_sites, 
+                                    tree_data_dir) {
+  
+  spp_names <- sort(spp_names)
+  # if the non base species names are zero return
+  if(length(spp_names)==0) return(base_species_df)
+
+  #============================================================================
+  # collect the tree data to match spp_names to 
+  # names <- c("idp", "espar", "a", "veget", "simplif", "acci", "ori", "lib", "forme", "tige", "mortb", "sfgui",
+  #            "sfgeliv", "sfpied", "sfdorge", "sfcoeur", "c13", "ir5", "htot", "hdec", "decoupe", "q1", "q2", "q3",
+  #            "r", "lfsd", "age", "v", "w", "old_espar", "frspecies", "xl93", "yl93", "tplant", "region", "dbh",
+  #            "altitude", "fudge_density", "stems", "basal_area", "scspecies", "species", "leaves", "stand_spp",
+  #            "stand_spp_not_chosen", "biomass", "start_dbh", "start_biomass"  )
+  # 
+  # tree_data <- setNames(data.frame(matrix(ncol = length(names), nrow = 0)), names)
+  # 
+  # for(j in 1:length(curr_sites)) {
+  # 
+  #   curr_data <- read.csv(file.path(tree_data_dir, paste0(curr_sites[j], "_tree_data.csv")))
+  #   common_cols <- intersect(colnames(curr_data), colnames(tree_data))
+  # 
+  #   tryCatch({
+  #     tree_data <- rbind(curr_data[common_cols], tree_data[common_cols])  
+  #   }, warning = function(warn) {
+  #     print("Warning recieved joining rows")
+  #   })
+  #   
+  # 
+  # }
+  
+  # for each spp. that doesn't exist in the base spp.
+  # for(i in 1:length(spp_names)) {
+  # 
+  #   # Get the leaf type for the current spp
+  #   curr_spp <- spp_names[i]
+  #   #print(paste0("current species is: ", curr_spp))
+  #   
+  #   leaf_type <- unique(tree_data[tree_data$species==curr_spp,]$leaves)
+  #   leaf_type <-factor(leaf_type, levels = c("DECIDUOUS", "EVERGREEN"))
+  #   #print(paste0("leaf type is: ", leaf_type))
+  #   rand_base <- base_species_df[sample( which( base_species_df$leafHabit==leaf_type), 1),]
+  #   rand_base$name <- spp_names[i]
+  #   print(rand_base)
+  #   all_species <- rbind(base_species_df, rand_base)
+  # }
+  spp_nm_tbl <- ifnHelpers::species_names_table()
+  #return(spp_nm_tbl)
+  # for each spp. that doesn't exist in the base spp.
+  
+  new_pars <- function(in_spp) {
+    
+    #print(class(spp_nm_tbl$species))
+    curr_spp <- as.character(in_spp)
+    #print(curr_spp)
+    curr_leaf_type <- subset(spp_nm_tbl, species==curr_spp)$leaves
+    #print(curr_leaf_type)
+    
+    if(is.null(curr_leaf_type)) {
+      stop(paste0("Current species ", curr_spp, " doesn't exist in ifnHelpers::species_names_table"))
+    }
+
+    leaf_type <- curr_leaf_type
+    #print(paste0("current species is: ", curr_spp))
+    #leaf_type <- unique(tree_data[tree_data$species==curr_spp,]$leaves)
+    #leaf_type <-factor(curr_leaf_type, levels = c("DECIDUOUS", "EVERGREEN"))
+    #print(paste0("leaf type is: ", leaf_type))
+    rand_base <- base_species_df[sample( which( base_species_df$leafHabit==leaf_type), 1),]
+    rand_base$name <- curr_spp
+    #print(rand_base)
+    #all_species <- rbind(base_species_df, rand_base)
+    return(rand_base)
+  }
+  
+  all_species <- do.call(rbind, lapply(spp_names, new_pars))
+  return(all_species)
+}

R/create_SED_XML_species.R

+create_SED_XML_species <- function(sim_is_dir,
+                                   file_prefix, 
+                                   sim_xml_dir,
+                                   curr_sites,
+                                   source_tree_dir) {
+  #print("in create_SED_XML_species")
+  # read in the initial state file from the simulation setup
+  # this step should have been completed first in create_SED_files
+  tryCatch({
+    is_file <- read.csv(file.path(sim_is_dir, paste0(file_prefix, "_tree_init.csv")),
+                        header = TRUE)
+  }, warning = function(warn){
+    print(paste0("Warning reading is_file: ", warn))
+  }, error = function(err) {
+    print(paste0("Error reading is_file: ", err))
+  })
+  #return(is_file)
+  # get the unique species names from that file i.e. all the species that 
+  # we need params for 
+  spp_names <- unique(is_file$species)
+  #return(spp_names)
+  # take in the Landclim species parameter files for parameterised species 
+  # and convert it to a dataframe
+  base_species <- XML::xmlToDataFrame(file.path(sim_xml_dir,
+                                                "LandClim_31species_CenEurope.xml"))
+  
+  # remove the rows from the base species that aren't in spp_names
+  rm_bs <- !(base_species$name %in% spp_names)
+  base_spp_sub <- base_species[!rm_bs,]
+  not_base_spp <- spp_names[!(spp_names %in% base_species$name)]
+  
+  # stop if the length of the names of base_spp and not_base_spp
+  # don't equal the length of the spp_names
+  stopifnot(length(base_spp_sub$name) + length(not_base_spp) == length(spp_names))
+  
+  # create dummy parameter values for those species that don't exist in 
+  # the base_spp
+  spp_fm_base <- create_IFN_species_pars(not_base_spp, 
+                                         base_spp_sub,
+                                         curr_sites, 
+                                         source_tree_dir)
+  
+  all_spp <- dplyr::arrange(rbind(base_spp_sub, spp_fm_base), as.character(name))
+  all_spp[, "foliageWeightFactor"] <- as.numeric(all_spp[, "foliageWeightFactor"])
+  all_spp[, "foliageWeightExponent"] <- as.numeric(all_spp[, "foliageWeightExponent"])
+  
+  # add the forrester foliage pars stored in the sysdata.rda as for_pars
+  for(i in 1:nrow(all_spp)) {
+  
+    params <- IFN_forrester_foliage_pars(as.character(all_spp[i ,"name"]),
+                                         as.character(all_spp[i,"leafHabit"]))
+    all_spp[i , "foliageWeightFactor"] <- unlist(params[1])
+    all_spp[i , "foliageWeightExponent"] <- params[2]
+    
+                
+  }
+  
+  #species_2_write <- subset(spp_fm_base, name %in% spp_names)
+  species_2_write <- all_spp
+  
+  write_species_xml(species_2_write,
+                    file.path(sim_xml_dir, paste0(file_prefix, "_species.xml")))
+  
+}
\ No newline at end of file

R/create_SED_climate_file.R

+# create SED climate file
+create_SED_climate_file <- function(no_sim_years, 
+                                    file_prefix, 
+                                    sim_clim_dir, 
+                                    source_clim_dir, 
+                                    source_land_dir, 
+                                    site_r) {
+
+  
+  #================================================
+  # FILE HEADERS
+  
+  # create the climate file headers to be read to
+  
+  # landclim currently needs 10 years of climate data
+  # if the simulation is < 10 years duplicate the missing years
+  tmp_sim_years <- ifelse(no_sim_years < 10, 10, no_sim_years)
+  
+  write_SED_climate_header(file_prefix,
+                           tmp_sim_years,
+                           sim_clim_dir)
+  
+  #================================================
+  # FILE TEMP, PRECIPITATION
+  ref_lats <- c()
+  
+  for(i in 1:raster::ncell(site_r)) {
+    
+    curr_site <- site_r[i]
+    curr_data <- read_climate_file(source_clim_dir, curr_site)
+    
+    # add the latitude to the list
+    ref_lats <- c(ref_lats, curr_data$lat)
+    
+    # retrieves the plots elevation raster
+    curr_elev_r <- raster::raster(file.path(source_land_dir, curr_site,
+                                            paste0(curr_site,
+                                                   "_elevation.asc")))
+    
+    # calculates and returns the geographically adjusted temp/precip data
+    adj_data <- lapse_adjust_climate(curr_data, raster::getValues(curr_elev_r))
+    
+    # get the coords for the specific site ...
+    adj_data_plus_rowcol <- row_col_to_clim_SED(site_r, i, adj_data)
+    
+    # output the climate data to the relevent SED file
+    write_SED_climate_data(adj_data_plus_rowcol,
+                           file_prefix,
+                           tmp_sim_years,
+                           sim_clim_dir)
+    }
+  
+  #================================================
+  # FILE LATITUDE
+
+  # Add latitude to the climate files
+  write_latitude_to_climate_file(file_prefix, tmp_sim_years,
+                                 sim_clim_dir, ref_lats)
+  
+}
\ No newline at end of file

R/create_SED_files.R

+#' @title Over riding script that pulls the SED creation together
+#' @description This should be what is called - not the individual functions
+#' @param file_prefix
+#' @param no_sim_years
+#' @sim_clim_dir
+#' @source_clim_dir
+#==============================================================================
+
+create_SED_files <- function(file_prefix,
+                             no_sim_years,
+                             sim_dir, 
+                             source_dir, 
+                             curr_sites,
+                             skip_landscape = FALSE,
+                             skip_climate = FALSE,
+                             skip_initialstate = FALSE,
+                             skip_XML_input = FALSE,
+                             skip_simulation = FALSE) {
+  
+  sim_clim_dir <- file.path(sim_dir, "climate")
+  sim_land_dir <- file.path(sim_dir, "landscape")
+  sim_is_dir <- file.path(sim_dir, "initialstate")
+  sim_xml_dir <- file.path(sim_dir, "xml_input")
+  sim_simulation_dir <- file.path(sim_dir, "simulations")
+
+  source_clim_dir <- file.path(source_dir, "climate")
+  source_land_dir <- file.path(source_dir, "landscape")
+  source_is_dir <- file.path(source_dir, "initialstate")
+  source_tree_dir <- file.path(source_dir, "tree_data")
+  
+  # prepare reference raster of the correct size
+  # with site names as values
+  # only use this now not curr_sites to be consistent with 
+  # this raster
+  site_r <- create_SED_raster(as.numeric(curr_sites))
+
+  #====================================
+  # LANDSCAPE
+  if(!skip_landscape) {
+    
+    # To make the relevent slope, aspect, elevation, soil depth rasters
+    make_SED_landscape(site_r, source_land_dir, sim_land_dir)
+    
+  } # End skip landscape
+  
+  #====================================
+  # CLIMATE
+  # Create the climate files ready to be populated
+  if(!skip_climate) {
+    
+    create_SED_climate_file(no_sim_years, 
+                            file_prefix, 
+                            sim_clim_dir, 
+                            source_clim_dir, 
+                            source_land_dir,
+                            site_r)
+    
+  } # End skip climate
+  
+  #====================================
+  # INITIAL STATE
+  if(!skip_initialstate){
+ 
+    # creates the SED initialstate file and returns the species names   
+    spp_names <- create_SED_initialstate_file(site_r,
+                                              source_is_dir,
+                                              sim_is_dir,
+                                              file_prefix)
+ 
+  } # End skip initial state
+  
+  #====================================
+  # XML Input files
+  #================
+  if(!skip_XML_input) {
+    
+    
+    create_SED_XML_species(sim_is_dir,
+                           file_prefix,
+                           sim_xml_dir,
+                           curr_sites,
+                           source_tree_dir)
+
+  } # End skip XML_input
+    
+  #====================================
+  # SIMULATION FILE
+  if(!skip_simulation) {
+    
+    create_SED_sim_file(sim_simulation_dir, 
+                        sim_xml_dir,
+                        file_root)
+
+  } # End skip simulation file
+  
+} # end of function
+
+

R/create_SED_initialstate_file.R

+create_SED_initialstate_file <- function(site_r,
+                                         source_is_dir,
+                                         sim_is_dir, 
+                                         file_prefix) {
+  
+  # create an empty initialstate file
+  init_file <- setup_init_file()
+  
+  # container for the species names
+  spp_names <- c()
+  
+  for(i in 1:raster::ncell(site_r)) {
+    
+    curr_site <- site_r[i]
+    init_file <- add_to_init_file(source_is_dir, 
+                                  curr_site,
+                                  init_file, 
+                                  i, 
+                                  site_r)
+  }
+  
+  # save the initialstate file
+  dir.create(sim_is_dir, recursive = TRUE, showWarnings = FALSE)
+  write.table(init_file, file.path(sim_is_dir, paste0(file_prefix, "_tree_init.csv")),
+              col.names = TRUE, sep = ",", row.names = FALSE, quote = FALSE)
+  
+