Added update_template.R, msr bug fixes

Author: wburke24

R/CreateFlownet.R

@@ -2,9 +2,9 @@
 #'
 #' Creates the flow networkd table used by RHESSys
 #' @param flownet_name The name of the flow network file to be created.  Will be coerced to have ".flow" extension if not already present.
-#' @param readin readin indicates which maps to be used. If CreateFlowmet.R is run it's own, this should point to the template. Otherwise,
-#' if run inside of RHESSysPreprocess, readin will use the map data from world_gen.R, Streams map, and other optional maps, still need to
-#' be specified.
+#' @param readin readin indicates the maps to be used. If CreateFlowmet.R is run it's own, this should point to the template which references the maps(or values)
+#' used for the various levels and statevars. Otherwise, if run inside of RHESSysPreprocess, readin will use the map data from world_gen.R,
+#' Streams map, and other optional maps, still need to be specified.
 #' @param type Input file type to be used. Default is raster. "Raster" type will use rasters
 #' in GeoTiff or equivalent format (see Raster package), with file names  matching those indicated in the template.
 #' ASCII is supported, but 0's cannot be used as values for data. "GRASS" will attempt to autodetect the version of
@@ -69,13 +69,13 @@ CreateFlownet = function(flownet_name,
   # Check for streams map, menu allows input of stream map
   if (is.null(streams) & (cfmaps[cfmaps[,1] == "streams",2] == "none" | is.na(cfmaps[cfmaps[,1] == "streams",2]))) {
     t = menu(c("Specify map","Abort function"),
-             title = "Missing stream map. Specify one now, or abort function and edit cf_maps file?")
+             title = "Missing stream map. Specify one now, or abort function and edit cf_maps file/readin input?")
     if (t == 2) {stop("Function aborted")}
     if (t == 1) {
       streams = readline("Stream map:")
     }
   }
-  # only add stream map to cfmaps if it's not there already
+  # add stream map to cfmaps if it's not there already
   if ((cfmaps[cfmaps[,1] == "streams",2] == "none" | is.na(cfmaps[cfmaps[,1] == "streams",2]))) {
     cfmaps[cfmaps[,1] == "streams",2] = streams
   }
@@ -88,32 +88,57 @@ CreateFlownet = function(flownet_name,
   if (!is.null(impervious)) {cfmaps[cfmaps[,1] == "impervious",2] = impervious}
   if (!is.null(roofs)) {cfmaps[cfmaps[,1] == "roofs",2] = roofs}
 
-  maps_in = unique(cfmaps[cfmaps[,2] != "none" & cfmaps[,1] != "cell_length",2])
+  # remove tif and tiff extensions for simplicity
+  if ( any(endsWith(cfmaps[,2],".tif") | endsWith(cfmaps[,2],".tiff")) ) {
+    cfmaps[,2] = gsub(".tif$","",cfmaps[,2])
+    cfmaps[,2] = gsub(".tiff$","",cfmaps[,2])
+  }
+
+  # check inpputs are maps or values
+  notamap = cfmaps[suppressWarnings( which(!is.na(as.numeric(cfmaps[,2])))),1]
+  maps_in = unique(cfmaps[cfmaps[,2] != "none" & !cfmaps[,1] %in% notamap,2])
 
   # ------------------------------ Use GIS_read to get maps ------------------------------
   readmap = GIS_read(maps_in, type, typepars, map_info = cfmaps)
+
+  # going to get rid of this
   map_ar_clean = as.array(readmap)
   dimnames(map_ar_clean)[[3]] = colnames(readmap@data)
 
-  raw_patch_data = map_ar_clean[, , cfmaps[cfmaps[, 1] == "patch", 2]]
-  raw_patch_elevation_data = map_ar_clean[, , unique(cfmaps[cfmaps[, 1] == "z", 2])]
-  raw_hill_data = map_ar_clean[, , cfmaps[cfmaps[, 1] == "hillslope", 2]]
-  raw_basin_data = map_ar_clean[, , cfmaps[cfmaps[, 1] == "basin", 2]]
-  raw_zone_data = map_ar_clean[, , cfmaps[cfmaps[, 1] == "zone", 2]]
-  raw_slope_data = map_ar_clean[, , unique(cfmaps[cfmaps[, 1] == "slope", 2])]
-  raw_stream_data = map_ar_clean[, , cfmaps[cfmaps[, 1] == "streams", 2]]
+  # this is better
+  map_list = lapply(readmap@data, matrix, nrow = readmap@grid@cells.dim[1], ncol = readmap@grid@cells.dim[2])
+
+  raw_patch_data = map_list[[cfmaps[cfmaps[, 1] == "patch", 2]]]
+  raw_hill_data = map_list[[cfmaps[cfmaps[, 1] == "hillslope", 2]]]
+  raw_basin_data = map_list[[cfmaps[cfmaps[, 1] == "basin", 2]]]
+  raw_zone_data = map_list[[cfmaps[cfmaps[, 1] == "zone", 2]]]
+  raw_stream_data = map_list[[cfmaps[cfmaps[, 1] == "streams", 2]]]
+
+  if ("slope" %in% notamap) {
+    raw_slope_data = raw_patch_data
+    raw_slope_data[!is.na(raw_slope_data)] = as.numeric(cfmaps[cfmaps[,1] == "slope",2])
+  } else {
+    raw_slope_data = map_list[[unique(cfmaps[cfmaps[, 1] == "slope", 2])]]
+  }
+  if ("z" %in% notamap) {
+    raw_patch_elevation_data = raw_patch_data
+    raw_patch_elevation_data[!is.na(raw_patch_elevation_data)] = as.numeric(cfmaps[cfmaps[,1] == "z",2])
+  } else {
+    raw_patch_elevation_data = map_list[[unique(cfmaps[cfmaps[, 1] == "z", 2])]]
+  }
+
   cell_length = readmap@grid@cellsize[1]
 
   # Roads
   raw_road_data = NULL
-  if (!is.null(roads)) {raw_road_data = map_ar_clean[, ,cfmaps[cfmaps[,1] == "roads",2]]}
+  if (!is.null(roads)) {raw_road_data =  map_list[[cfmaps[cfmaps[,1] == "roads",2]]]}
 
   # Roofs and impervious is not yet implemented - placeholders for now -----
   if (!is.null(roofs) | !is.null(impervious)) {print("Roofs and impervious are not yet working",quote = FALSE)}
   raw_roof_data = NULL
-  if (!is.null(roofs)) {raw_roof_data = map_ar_clean[, ,cfmaps[cfmaps[,1] == "roofs",2]]}
+  if (!is.null(roofs)) {raw_roof_data =  map_list[[cfmaps[cfmaps[,1] == "roofs",2]]]}
   raw_impervious_data = NULL
-  if (!is.null(impervious)) {raw_impervious_data = map_ar_clean[, ,cfmaps[cfmaps[,1] == "impervious",2]]}
+  if (!is.null(impervious)) {raw_impervious_data =  map_list[[cfmaps[cfmaps[,1] == "impervious",2]]]}
 
   # ----- SMOOTH FLAG STILL NEEDS TO BE LOOKED AT AGAIN -----
   smooth_flag = FALSE
@@ -138,7 +163,12 @@ CreateFlownet = function(flownet_name,
 
   # ------------------------------ Multiscale routing/aspatial patches ------------------------------
   if (!is.null(asp_list)) {
-    CF1 = multiscale_flow(CF1 = CF1, map_ar_clean = map_ar_clean, cfmaps = cfmaps, asp_list = asp_list)
+    if ("asp_rule" %in% notamap) {
+      map_list[["asp_rule"]] = raw_basin_data
+      map_list[["asp_rule"]][!is.na(map_list[["asp_rule"]])] = as.numeric(cfmaps[cfmaps[,1] == "asp_rule",2])
+    }
+
+    CF1 = multiscale_flow(CF1 = CF1, map_list = map_list, cfmaps = cfmaps, asp_list = asp_list)
   }
 
   # ---------- Flownet list to flow table file ----------

R/multiscale_flow.R

@@ -2,21 +2,24 @@
 # function that transforms flowtable into multiscale flow table
 # Will Burke 1/16/19
 
-multiscale_flow = function(CF1, map_ar_clean, cfmaps, asp_list){
+multiscale_flow = function(CF1, map_list, cfmaps, asp_list){
 
   # Inputs: exisiting flow table (list), map arrays, map list, aspatial/multiscale rule list
 
   # ----- Variable setup -----
-  asp_map = map_ar_clean[, ,cfmaps[cfmaps[,1] == "asp_rule",2]] # matrix of aspatial rules
+  asp_map = map_list[["asp_rule"]] # matrix of aspatial rules
+
   patch_ID = unlist(lapply(CF1, "[[",9)) # patch IDs from cf1
   numbers = unlist(lapply(CF1, "[[",1)) # flow list numbers
-  raw_patch_data = map_ar_clean[, , cfmaps[cfmaps[, 1] == "patch", 2]] # get patch matrix inside the function
+  raw_patch_data = map_list[[cfmaps[cfmaps[, 1] == "patch", 2]]] # get patch matrix inside the function
   rulevars = asp_list # get rules - state variable overrides
   CF2 = list() # empty list for new flow list
 
   # ----- iterate through (spatial) patches -----
-  for (p in raw_patch_data[!is.na(raw_patch_data)]) {
+  #for (p in raw_patch_data[!is.na(raw_patch_data)]) {
+   for (p in patch_ID) {
     id = unique(asp_map[which(raw_patch_data == p)]) # get unique rule ID for patch p
+    id = paste0("rule_",id)
 
     if (length(id) > 1) {stop(paste("multiple aspatial rules for patch",p))} # if multiple rules for a single patch
     asp_count = ncol(rulevars[[id]]$patch_level_vars[1,]) - 1 # get number of aspatial patches for current patch
@@ -60,12 +63,13 @@ multiscale_flow = function(CF1, map_ar_clean, cfmaps, asp_list){
       CF2[[length(CF2)]]$Slopes = new_slope
       CF2[[length(CF2)]]$Boarder = new_boarder
 
-      for(i in 1:length(CF2[[length(CF2)]])){
-        if(is.na(CF2[[length(CF2)]])[[i]]) {stop("shouldn't have NAs")}
-        if(length(CF2[[length(CF2)]][[i]])==0) {stop("shouldn't have numeric(0)'s")}
+      if (length(CF2) != 1 & unique(sapply(CF2, "[[", "PatchFamilyID")) != 1) {
+        for (i in 1:length(CF2[[length(CF2)]])) {
+          if(is.na(CF2[[length(CF2)]])[[i]]) {stop("shouldn't have NAs")}
+          if(length(CF2[[length(CF2)]][[i]])==0) {stop("shouldn't have numeric(0)'s")}
+        }
       }
 
-
     } # end aspatial patch loop
   } # end spatial patch loop
 

R/patch_data_analysis.R

@@ -396,7 +396,7 @@ patch_data_analysis <- function(raw_patch_data,
 
   # -------------------- Single Patch World --------------------
   # If there's only 1 patch, build list manually and return
-  if (length(raw_patch_data) == 1) {
+  if (length(raw_patch_data) == 1 | length(flw_struct$Number) == 1) {
     lst[[1]]$Slopes = 0
     lst[[1]]$TotalG = flw_struct$Area
     return(lst)

R/update_template.R

@@ -0,0 +1,97 @@
+#' update_template
+#' 
+#' Read and modify an existing template
+#' @param template Name and path of templatefile to be created.
+#' @param out_file Destination file to write new template
+#' @param vars State variable(s) to be modified
+#' @param operator Operator to be (optionally) replaced. If  not included the existing operator (e.g. 'value') will be left.
+#' @param values Values to replace for `vars`. N args must equal that of vars or be 1 (will be used for all vars).
+#' @param level_names Not implamented yet - for rare cases where you only want to replace 1 occurance of a state var (e.g. 'z' at zone but not patch ele)
+#' @param overwrite TRUE/FALSE if input file should be overwriten
+#' @author Will Burke
+#' @export
+
+update_template = function(template, out_file, vars, values, operator = NULL, level_names = NULL, overwrite = FALSE) {
+  
+  # NOTES
+  # - vdouble check all the list options work and check var lengths/numbers when needed
+  
+  # ---------- Check Aguments ----------
+  if (is.null(out_file) & !overwrite) {stop(noquote("No destination file set by 'out_file' and 'overwrite' is FALSE"))}
+  if ( file.exists(out_file) & overwrite == FALSE) {stop(noquote(paste0("File '",out_file,"' already exists and 'overwrite' argument is FALSE")))}
+  if (length(vars) != length(values) & length(values) != 1) { stop(noquote(paste0("Vars and values have different lengths"))) }
+    
+
+  if (!is.character(values)) {
+    values = as.character(values)
+  }
+  
+  
+  # for (i in length(vars)) {
+  #   if (level_names[i] == "all") {level_names[i] = c("world", "basin", "hillslope", "zone", "patch", "canopy_strata")}
+  # }
+  
+  # ---------- Parse templatefile ----------
+  options(scipen = 999)
+  # parsing the values as characters to retain the exact value/precision
+  read_template = readLines(template)
+  template =  strsplit(trimws(read_template), "\\s+")
+  
+  #template = data.frame(matrix(unlist(template), nrow=length(template), byrow=T),stringsAsFactors = FALSE)
+  #names(template) = c("values","vars")
+  
+  # ---------- Find Levels----------
+  template_vars = sapply(template,"[[",1)
+  index_all = which(template_vars == "_world" | template_vars == "_basin" | template_vars == "_hillslope" | 
+                      template_vars == "_zone" | template_vars == "_patch" | template_vars == "_canopy_strata")
+  index_names = gsub("^_", "", x = template_vars[index_all])
+  index_max = c(index_all[2:length(index_all)]-1, length(read_template))
+  template_levels = unname(unlist(mapply(rep,index_names, (index_max - index_all) + 1 )))
+  
+  # ---------- Find and Replace Vars ----------
+  for (i in length(vars)){
+    
+    if (!is.null(level_names[i])) {
+      find_index = template_levels == level_names[i]
+    } else {
+      find_index = rep(TRUE,length(template_levels))
+    }
+
+    replace_index = which(template_vars == vars[i] & find_index)
+    
+    # if unique values for every instance of var to be replaces were given, do nothing, otherwise repeat to get enough replacement values
+    if (length(replace_index) > 1) {stop(noquote("Too many things to replace havet coded this yet"))}
+    
+    if (template_levels[replace_index] == "canopy_strata") {
+      current_value = template[[replace_index]][c(length(template[[replace_index]]) - 1, length(template[[replace_index]]))]
+    } else {
+      current_value = template[[replace_index]][length(template[[replace_index]])]
+    }
+    
+    if (length(values[i]) != length(replace_index)) {
+      new_value = rep(values[i], length(replace_index)/length(values[i]))
+    } else {
+      new_value = values[i]
+    }
+    
+    if (any(startsWith(new_value,"*"))) {
+      new_value = as.numeric(trimws(substr(new_value[startsWith(new_value,"*")],2,nchar(new_value[startsWith(new_value,"*")]))))
+      new_value = new_value * as.numeric(current_value)
+    } else {
+      new_value = as.numeric(trimws(new_value))
+    }
+    
+    # generic sub/gsub
+    #sub(paste0("[[:blank:]]",current_value,"$"), paste0("\t",new_value), read_template[replace_index])
+    
+    read_template[replace_index] = unname(mapply(sub,paste0("[[:blank:]]",current_value,"$"), paste0("\t",new_value),read_template[replace_index]))
+    
+  }
+  
+  # ---------- Write file ----------
+  writeLines(text = read_template,out_file)
+  
+  print(noquote(paste("Successfully wrote updated template to",out_file)))
+  
+  
+}

R/world_gen.R

@@ -295,15 +295,20 @@ world_gen = function(template,
         if (asp_check) {
           patches = unique(levels[levels[,4] == z & levels[,3] == h & levels[,2] == b, 5])
           asp_ct = sapply(rulevars, FUN = function(x) ncol(x[[1]]) - 1)
-          total_patches = sum(asp_ct[levels[levels[,5] == patches & levels[,4] == z & levels[,3] == h & levels[,2] == b, 7]])
+          if (length(patches) == 1 & length(asp_ct) == 1){
+            total_patches = length(patches) * asp_ct
+          } else {
+            total_patches = sum(asp_ct[levels[levels[,5] == patches & levels[,4] == z & levels[,3] == h & levels[,2] == b, 7]])
+          }
 
           writeChar(paste("\t\t\t",total_patches,"\t\t\t","num_patches\n",sep = ""),con = wcon,eos = NULL)
 
           # ----- Patches (spatial) -----
           for (p in patches) {
             ruleid = unique(levels[(levels[,5] == p & levels[,4] == z & levels[,3] == h & levels[,2] == b),7])
+            ruleid = paste0("rule_",ruleid)
             if (length(ruleid) != 1) {stop("something's wrong with the ruleid")}
-            asp_index = 1:(length(rulevars[[ruleid]]$patch_level_vars[1,]) - 1)
+            asp_index = 1:(length(rulevars[[(ruleid)]]$patch_level_vars[1,]) - 1)
 
             # ----- Patches (non-spatial) -----
             for (asp in asp_index) {
@@ -443,6 +448,19 @@ world_gen = function(template,
                  c("cell_length",cell_len),
                  c("streams","none"), c("roads","none"), c("impervious","none"),c("roofs","none"))
 
+  if (any(!map_info[,1] == "slope")) {
+    slope = NULL
+    for (i in 1:length(which(var_names == "slope"))) {
+      slope = c(slope, template_clean[[which(var_names == "slope")[i]]][3])
+    }
+    cfmaps = rbind(cfmaps, c("slope", mean(as.numeric(slope))))
+  }
+  if (asp_check) {
+    if (!"asp_rule" %in% cfmaps[,1]) {
+      cfmaps = rbind(cfmaps, c("asp_rule", mean(as.numeric(template_clean[[which(var_names == "asp_rule")]][3]))))
+    }
+  }
+
   world_gen_out = list(cfmaps,rulevars)
 
   return(world_gen_out)