Various flownet changes/fixes

wburke24 · wburke24 · commit d4ee85c5f1ea · 2021-05-03T10:16:41.000-07:00
diff --git a/R/CreateFlownet.R b/R/CreateFlownet.R
@@ -5,7 +5,7 @@
 #' @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 asp_list List of aspatial structure and inputs. Used internally
+#' @param asp_rules List of aspatial structure and inputs. Also can be path to rules file - must be used along with template input.
 #' @param road_width >0, defaults to 1.
 #' @inheritParams RHESSysPreprocess
 #' @author Will Burke
@@ -15,7 +15,7 @@ CreateFlownet = function(flownet_name,
                          readin = NULL,
                          type = "raster",
                          typepars = NULL,
-                         asp_list = NULL,
+                         asp_rules = NULL,
                          streams = NULL,
                          overwrite = FALSE,
                          roads = NULL,
@@ -24,6 +24,7 @@ CreateFlownet = function(flownet_name,
                          roofs = NULL,
                          parallel = TRUE,
                          make_stream = 4,
+                         skip_hillslope_check = FALSE,
                          wrapper = FALSE) {
 
   # ------------------------------ Read and check inputs ------------------------------
@@ -40,7 +41,9 @@ CreateFlownet = function(flownet_name,
 
   if (!wrapper & is.character(readin)) { #if run outside of rhessyspreprocess.R, and if readin is character. readin is the template (and path)
     template_list = template_read(readin)
-    map_info = template_list[[5]]
+    template_clean = template_list[[1]] # template in list form
+    var_names = template_list[[2]] # names of template vars
+    map_info = template_list[[5]] # tables of maps and their inputs/names in the template
     cfmaps = rbind(map_info,c("cell_length","none"), c("streams","none"), c("roads","none"), c("impervious","none"),c("roofs","none"))
   } else if (wrapper | (!wrapper & is.matrix(readin))) { # map info is passsed directly from world gen - either in wrapper or outside of wrapper and readin is matrix
     cfmaps = readin
@@ -113,6 +116,18 @@ CreateFlownet = function(flownet_name,
   raw_impervious_data = NULL
   if (!is.null(impervious)) {raw_impervious_data =  map_list[[cfmaps[cfmaps[,1] == "impervious",2]]]}
 
+  # read aspatial rules if needed
+  if (!is.null(asp_rules) && is.character(asp_rules)) {
+    asp_map = template_clean[[which(var_names == "asp_rule")]][3] # get rule map/value
+    if (suppressWarnings(any(is.na(as.numeric(asp_map))))) { # if it's a map
+      asp_map = gsub(".tif|.tiff","",asp_map)
+      asp_mapdata = as.data.frame(readmap)[asp_map]
+    } else if (suppressWarnings(all(!is.na(as.numeric(asp_map))))) { # if is a single number
+      asp_mapdata = as.numeric(asp_map)
+    }
+    asp_rules = aspatial_patches(asprules = asp_rules, asp_mapdata = asp_mapdata)
+  }
+
   # ------------------------------ Make flownet list ------------------------------
   cat("Building flowtable")
   CF1 = make_flow_list(
@@ -127,18 +142,19 @@ CreateFlownet = function(flownet_name,
     road_width = road_width,
     cell_length = cell_length,
     parallel = parallel,
-    make_stream = make_stream)
+    make_stream = make_stream,
+    skip_hillslope_check = skip_hillslope_check)
 
   # ------------------------------ Multiscale routing/aspatial patches ------------------------------
-  if (!is.null(asp_list)) {
+  if (!is.null(asp_rules)) {
     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])
       #cfmaps = rbind(cfmaps, c("asp_rule", "asp_rule"))
       cfmaps[cfmaps[,1] == "asp_rule", 2] = "asp_rule"
     }
 
-    CF1 = multiscale_flow(CF1 = CF1, map_list = map_list, cfmaps = cfmaps, asp_list = asp_list)
+    CF1 = multiscale_flow(CF1 = CF1, map_list = map_list, cfmaps = cfmaps, asp_list = asp_rules)
   }
 
   # ---------- Flownet list to flow table file ----------
diff --git a/R/GIS_read.R b/R/GIS_read.R
@@ -99,8 +99,14 @@ GIS_read = function(maps_in, type, typepars, map_info = NULL, seq_patch_IDs = FA
       if (length(file) == 0) { # if there were no matches
         stop(paste("No file named:",name,"at path:",typepars))
       }
+      if (length(file) > 1) { # ignore .xml files
+        file = file[!grepl(".xml$",file)]
+      }
       if (length(file) > 1) { # if multiple files, use tif preferentially
         file = file[grep(".tif$",file)]
+        if (length(file) > 1) {
+          file = file[grep(".tiff$",file)]
+        }
       }
       if (length(file) > 1) { # if STILL multiple files, can only be one file for each name in maps_in
         stop(paste("multiple files containing name:",name,"check directory:",typepars))
diff --git a/R/RHESSysPreprocess.R b/R/RHESSysPreprocess.R
@@ -35,6 +35,8 @@
 #' @param make_stream The maximum distance (cell lengths) away from an existing stream that a patch can be automatically coerced to be a stream.
 #' Setting to TRUE will include patches at any distance. This is needed for hillslope parallelization, as all hillslopes must have an outlet stream patch.
 #' Default is 4.
+#' @param skip_hillslope_check TRUE/FALSE to skip the recursive check for segmented hillslopes. Segmented hillslopes will break the routing, but the
+#' recursive check can trigger various recursion protections when hillslopes are large.
 #' @param wrapper internal argument to track if being run as all-in-one
 #' @seealso \code{\link[raster]{raster}}, \code{\link[RHESSysIOinR]{run_rhessys}}
 #' @author Will Burke
@@ -59,6 +61,7 @@ RHESSysPreprocess = function(template,
                              fire_grid_out = FALSE,
                              parallel = TRUE,
                              make_stream = 4,
+                             skip_hillslope_check = FALSE,
                              wrapper = TRUE) {
 
   # ---------- Check Inputs ----------
diff --git a/R/make_flow_list.R b/R/make_flow_list.R
@@ -20,17 +20,18 @@
 #' @importFrom utils setTxtProgressBar
 
 make_flow_list <- function(raw_patch_data,
-                                raw_patch_elevation_data,
-                                raw_hill_data,
-                                raw_basin_data,
-                                raw_zone_data,
-                                raw_slope_data,
-                                raw_stream_data,
-                                raw_road_data,
-                                cell_length,
-                                road_width = NULL,
-                                parallel,
-                                make_stream) {
+                           raw_patch_elevation_data,
+                           raw_hill_data,
+                           raw_basin_data,
+                           raw_zone_data,
+                           raw_slope_data,
+                           raw_stream_data,
+                           raw_road_data,
+                           cell_length,
+                           road_width = NULL,
+                           parallel,
+                           make_stream,
+                           skip_hillslope_check) {
 
   # -------------------- Error checking and NULL handling --------------------
   if (cell_length <= 0) {stop("Cell length is <=0")}
@@ -258,8 +259,7 @@ make_flow_list <- function(raw_patch_data,
     }
 
 
-
-    find_connected = function(patch_borders, start, history = NULL, found = NULL, missing, hill) {
+    find_connected = function(patch_borders, start, history = NULL, found = NULL, missing, hill, itr_ct = NULL, itr_max = 10000) {
       # add start to found, only for first itr tho
       if (is.null(history)) {
         found = start
@@ -268,6 +268,7 @@ make_flow_list <- function(raw_patch_data,
       history = c(history, start)
       # find the neighbors
       neighbors = names(patch_borders[[start]])[-1][names(patch_borders[[start]])[-1] %in% flw_struct[flw_struct$Hill == hill, "Number"]]
+      neighbors = as.numeric(neighbors)
       # update the found vector to include those neighbors
       found = c(found, neighbors)
       # remove the found from missing
@@ -280,62 +281,81 @@ make_flow_list <- function(raw_patch_data,
       if (all(found %in% history)) {
         return(missing)
       }
+
+      if (!is.null(itr_ct)) {
+        itr_ct = itr_ct + 1
+        if (itr_ct == itr_max) {
+          cat("Reached iteration limit of ",itr_max)
+          return(list(missing, found, history, start))
+        }
+      }
+
       # if we haven't found everything, keep iterating through found patches not in the history
-      next_patch = as.numeric(found[!found %in% history][1])
+      opts = as.numeric(found[!found %in% history])
+      next_patch = opts[round(runif(1,min = 1, max = length(opts)))]
+      #next_patch = as.numeric(found[!found %in% history][1])
+
       find_connected(
         patch_borders = patch_borders,
         start = next_patch,
         history = history,
         found = found,
         missing = missing,
-        hill = hill
+        hill = hill,
+        itr_ct = itr_ct,
+        itr_max = itr_max
       )
     }
 
-    # check for segmented hillslopes - this is going to be so slow but oh well
-    cat("Checking for segmented hillslopes")
-    segmented_hills = NULL
-    segmented_patch_ct = NULL
-    segmented_patch_list = list()
-
-    pb = txtProgressBar(min = 0,max = length(unique(flw_struct$Hill)),style = 3)
-    ct = 0
+    if (!skip_hillslope_check) {
+      # check for segmented hillslopes - this is going to be so slow but oh well
+      cat("Checking for segmented hillslopes")
+      segmented_hills = NULL
+      segmented_patch_ct = NULL
+      segmented_patch_list = list()
 
-    for (i in unique(flw_struct$Hill)) {
-      ct = ct + 1
-      setTxtProgressBar(pb,ct)
-      # start from outlet
-      min_patch = flw_struct[flw_struct$Hill == i,][min_hill_patch[min_hill_patch$hillID == i, "patch_ind"],]
-      all_patches = flw_struct[flw_struct$Hill == i, "Number"]
+      pb = txtProgressBar(min = 0,max = length(unique(flw_struct$Hill)),style = 3)
+      ct = 0
 
-      missing = find_connected(
-        patch_borders = patch_borders,
-        start = min_patch$Number,
-        history = NULL,
-        found = NULL,
-        missing = all_patches,
-        hill = i
-      )
+      # for testing:
+      i = unique(flw_struct$Hill)[4]
 
-      if (!is.null(missing)) {
-
-        segmented_hills = c(segmented_hills, i)
-        segmented_patch_ct = c(segmented_patch_ct, length(missing))
-        segmented_df = data.frame("Hillslope ID" = segmented_hills, "Segmented patch count" = segmented_patch_ct)
-        segmented_patch_list[[as.character(i)]] = missing
-
-        warning(
-          "One or more hillslopes are disconnected with segmented segmented patches and do not reach the outlet patch.\n",
-          "This will lead to problems in your simulaion. See table below for hilslope IDs and segmented patch counts.\n",
-          "Regenerate hillslopes with a higher threshold for accumulated upslope area to correct this.\n"
+      for (i in unique(flw_struct$Hill)) {
+        ct = ct + 1
+        setTxtProgressBar(pb,ct)
+        # start from outlet
+        min_patch = flw_struct[flw_struct$Hill == i,][min_hill_patch[min_hill_patch$hillID == i, "patch_ind"],]
+        all_patches = flw_struct[flw_struct$Hill == i, "Number"]
+
+        missing = find_connected(
+          patch_borders = patch_borders,
+          start = min_patch$Number,
+          history = NULL,
+          found = NULL,
+          missing = all_patches,
+          hill = i,
+          itr_ct = 0,
+          itr_max = length(all_patches)
         )
-        print(segmented_df)
-        #error("")
 
-      }
+        if (!is.null(missing)) {
+          segmented_hills = c(segmented_hills, i)
+          segmented_patch_ct = c(segmented_patch_ct, length(missing))
+          segmented_df = data.frame("Hillslope ID" = segmented_hills, "Segmented patch count" = segmented_patch_ct)
+          segmented_patch_list[[as.character(i)]] = missing
+          warning(
+            "One or more hillslopes are disconnected with segmented segmented patches and do not reach the outlet patch.\n",
+            "This will lead to problems in your simulaion. See table below for hilslope IDs and segmented patch counts.\n",
+            "Regenerate hillslopes with a higher threshold for accumulated upslope area to correct this.\n"
+          )
+          print(segmented_df)
+          #error("")
+
+        }
 
-    } # end for segmented hills
-    close(pb)
+      } # end for segmented hills
+      close(pb)
+    }
 
 
   } # end parallel if
diff --git a/R/multiscale_flow.R b/R/multiscale_flow.R
