diff --git a/R/GIS_read.R b/R/GIS_read.R index 0598b12..a59aa45 100644 --- a/R/GIS_read.R +++ b/R/GIS_read.R @@ -144,9 +144,10 @@ GIS_read = function(maps_in, type = "raster", typepars, map_info = NULL, seq_pat NAs_in_wrld[[map_info[map_info[, 1] == "streams", 2]]] = NULL } if (any(NAs_in_wrld > 0) ) { - cat("One or more maps have NAs within the bounds of the world map, see maps and counts of NAs below:\n") + cat("Warning: One or more maps have NAs within the bounds of the world map, see maps and counts of NAs below:\n") print(NAs_in_wrld[NAs_in_wrld > 0]) - stop("See above and check your input maps.") + #09282022: LML for some layers, it's okay to have NA values + #stop("See above and check your input maps.") } } diff --git a/R/RHESSysPreprocess.R b/R/RHESSysPreprocess.R index f7fd51f..68f9436 100644 --- a/R/RHESSysPreprocess.R +++ b/R/RHESSysPreprocess.R @@ -46,6 +46,7 @@ RHESSysPreprocess = function(template, streams = NULL, overwrite = FALSE, roads = NULL, + road_width = NULL, impervious = NULL, roofs = NULL, header = FALSE, @@ -78,12 +79,12 @@ RHESSysPreprocess = function(template, basename = substr(basename, 0, nchar(basename) - 5) } name_clean = file.path(dirname(name), basename) - worldfile = name_clean - flownet_name = name_clean + worldfile = paste(name_clean, ".world", sep = "") + flownet_name = paste(name_clean,".flow", sep = "") if (!dir.exists(dirname(name))) { # check if output dir exists, menu to create t = utils::menu( - c("Yes", "No [Exit]"), + c("Yes", "No [Skip]"), title = paste("Ouput directory path:",dirname(name),"is not valid. Create folder(s)?")) if (t == 1) { dir.create(dirname(name), recursive = TRUE) @@ -103,24 +104,26 @@ RHESSysPreprocess = function(template, # ---------- Run world_gen ---------- cat("Begin world_gen.R\n") - + t = 1 if (file.exists(worldfile) & overwrite == FALSE) { # check for worldfile overwrite - t = utils::menu(c("Yes", "No [Exit]"), title = noquote(paste( + t = utils::menu(c("Yes", "No [Skip]"), title = noquote(paste( "Worldfile", worldfile, "already exists. Overwrite?" ))) if (t == 2) { - stop("RHESSysPreprocess.R exited without completing") + #stop("RHESSysPreprocess.R exited without completing") + cat("Use existing world file\n") } } - - world_gen_out = world_gen(template = template, + if (t == 1) { + world_gen_out = world_gen(template = template, worldfile = worldfile, type = type, typepars = typepars, - overwrite = overwrite, + overwrite = TRUE, header = header, unique_strata_ID = unique_strata_ID, asprules = asprules) + } #readin = world_gen_out[[1]] #asp_rules = world_gen_out[[2]] @@ -128,29 +131,33 @@ RHESSysPreprocess = function(template, # ---------- Run create_flownet ---------- cat("Begin create_flownet.R") + t = 1 if (file.exists(flownet_name) & overwrite == FALSE) { # check for flownet overwrite - t = utils::menu(c("Yes", "No [Exit]"), title = noquote(paste( + t = utils::menu(c("Yes", "No [Skip]"), title = noquote(paste( "Flowtable", flownet_name, "already exists. Overwrite?" ))) if (t == 2) { - stop("RHESSysPreprocess.R exited without completing") + #stop("RHESSysPreprocess.R exited without completing") + cat("Use existing flownet file\n") } } - - create_flownet(flownet_name = flownet_name, + if (t == 1) { + create_flownet(flownet_name = flownet_name, template = template, type = type, typepars = typepars, asprules = asprules, streams = streams, - overwrite = overwrite, + overwrite = TRUE, roads = roads, + road_width = road_width, impervious = impervious, roofs = roofs, wrapper = wrapper, parallel = parallel, make_stream = make_stream, skip_hillslope_check = skip_hillslope_check) + } # ---------- Run build_meta ---------- # if (meta) { diff --git a/R/create_flownet.R b/R/create_flownet.R index 02412d3..3f48fd7 100644 --- a/R/create_flownet.R +++ b/R/create_flownet.R @@ -92,12 +92,14 @@ create_flownet = function(flownet_name, } 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 @@ -120,7 +122,6 @@ create_flownet = function(flownet_name, if (!is.null(asprules)) { asp_map = template_clean[[which(var_names == "asp_rule")]][3] # get rule map/value patch_map = map_df[[cfmaps[cfmaps[,1] == "patch",2]]] # set for use later - overwrite if using mode - if (suppressWarnings(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] @@ -177,7 +178,7 @@ create_flownet = function(flownet_name, } # ------------------------------ Make flownet list ------------------------------ - cat("Building flowtable") + cat("Building flowtable\n") CF1 = make_flow_list( raw_patch_data = raw_patch_data, raw_patch_elevation_data = raw_patch_elevation_data, diff --git a/R/find_stream.R b/R/find_stream.R index d8ecbee..9a2781f 100644 --- a/R/find_stream.R +++ b/R/find_stream.R @@ -1,40 +1,42 @@ -# Searches through nodes, if a node has a road, it finds the nearest stream in the nodes descendants. -# If there are no streams it uses the final basin (sum of gamma ==0). Once the stream is found, the -# landtype is changed to 2, and the roadtype variable gets the node number. -find_stream<-function(flw,road_width){ - len_flw<-length(flw) - for (i in 1:len_flw){ - if ((flw[[i]]$Roadtype==1)&(sum(flw[[i]]$Gamma_i)!=0)){ - - # returns a list of all nodes that can be reached from start_node - len<-length(flw) - clist_bool<-rep(FALSE,len) - clist_bool<-find_children(flw,i,clist_bool) - if (sum(clist_bool)==0){ #if there are no downstream nodes, return 0 - desc = 0 - } else { # otherwise, return a list of all downstream nodes. - desc<-seq.int(1,len)[clist_bool] - } - - strm_chld<-c() #list of descendants with a stream or bottom - strm_dist<-c() # 2d distance to each stream - for (j in desc) { - if ((flw[[j]]$Landtype==1)|(sum(flw[[j]]$Gamma_i)==0)) { #has a stream or is bottom - strm_chld<-c(strm_chld,j) - # find 2d distance between two nodes - strm_dist<-c(strm_dist,sqrt((flw[[i]]$Centroidx-flw[[j]]$Centroidx)^2+(flw[[i]]$Centroidy-flw[[j]]$Centroidy)^2)) - } - } - strm_chld<-strm_chld[order(strm_dist)] # order list by 2d distance - flw[[i]]$Roadtype<-strm_chld[[1]] #pick closest one - } - } - for (i in 1:len_flw){ - if (flw[[i]]$Roadtype!=0){ - strm<-flw[[i]]$Roadtype - flw[[i]]$Landtype<-2 - flw[[i]]$Roadtype<-c(flw[[strm]]$PatchID,flw[[strm]]$ZoneID,flw[[strm]]$HillID,road_width) - } - } - return(flw) -} +# Searches through nodes, if a node has a road, it finds the nearest stream in the nodes descendants. +# If there are no streams it uses the final basin (sum of gamma ==0). Once the stream is found, the +# landtype is changed to 2, and the roadtype variable gets the node number. +find_stream<-function(flw,road_width){ + len_flw<-length(flw) + for (i in 1:len_flw){ + if ((flw[[i]]$Roadtype==1)&(sum(flw[[i]]$Gamma_i)!=0)){ + + # returns a list of all nodes that can be reached from start_node + len<-length(flw) + clist_bool<-rep(FALSE,len) + clist_bool<-find_children(flw,i,clist_bool) + if (sum(clist_bool)==0){ #if there are no downstream nodes, return 0 + desc = 0 + } else { # otherwise, return a list of all downstream nodes. + desc<-seq.int(1,len)[clist_bool] + } + + strm_chld<-c() #list of descendants with a stream or bottom + strm_dist<-c() # 2d distance to each stream + for (j in desc) { + if ((flw[[j]]$Landtype==1)|(sum(flw[[j]]$Gamma_i)==0)) { #has a stream or is bottom + strm_chld<-c(strm_chld,j) + # find 2d distance between two nodes + strm_dist<-c(strm_dist,sqrt((flw[[i]]$Centroidx-flw[[j]]$Centroidx)^2+(flw[[i]]$Centroidy-flw[[j]]$Centroidy)^2)) + } + } + if(length(strm_dist) > 0) { + strm_chld<-strm_chld[order(strm_dist)] # order list by 2d distance + flw[[i]]$Roadtype<-strm_chld[[1]] #pick closest one + } + } + } + for (i in 1:len_flw){ + if (flw[[i]]$Roadtype!=0){ + strm<-flw[[i]]$Roadtype + flw[[i]]$Landtype<-2 + flw[[i]]$Roadtype<-c(flw[[strm]]$PatchID,flw[[strm]]$ZoneID,flw[[strm]]$HillID,road_width) + } + } + return(flw) +} diff --git a/R/make_flow_list.R b/R/make_flow_list.R index 57c43ce..6756235 100644 --- a/R/make_flow_list.R +++ b/R/make_flow_list.R @@ -47,7 +47,7 @@ make_flow_list <- function(raw_patch_data, options(scipen = 999) # IF DEBUGGING YOU WILL HAVE ERRORS ON LARGE BASINS WITHOUT THIS - comes from numeric to character conversion # -------------------- Build unique patch IDs -------------------- - cat("Generating unique patch IDs") + cat("Generating unique patch IDs\n") id_data = data.frame("Basin" = as.vector(raw_basin_data[!is.na(raw_basin_data)]), "Hill" = as.vector(raw_hill_data[!is.na(raw_basin_data)]), @@ -121,6 +121,8 @@ make_flow_list <- function(raw_patch_data, flw_struct <- data.frame(patches,xy_data,patch_mean_elev,patch_mean_slope,patch_landtype,patch_roadtype) colnames(flw_struct) <- c("Number","Centroidx","Centroidy","Area","Centroidz","Mean_Slope","Landtype","Roadtype") + #12212022LML if patch is both road and stream, set it as stream + flw_struct$Roadtype[flw_struct$Roadtype == 1 & flw_struct$Landtype == 1] <- 0 flw_struct <- cbind(flw_struct,id_data_unique[-5]) # -------------------- find neighbors cell by cell (previously find border row) -------------------- @@ -128,7 +130,7 @@ make_flow_list <- function(raw_patch_data, # patch_borders is an array, initailly 0 but will be filled with the number of times patch i # touches patch j. the diagonal will be the number of times patch i touches anything. # -------------------- D8 neighbor search and border count -------------------- - cat("Finding patch neighbors") + cat("Finding patch neighbors\n") # new - list instead of matrix patch_borders = list(list("Total" = 0))[rep(1,length(patches))] @@ -225,7 +227,7 @@ make_flow_list <- function(raw_patch_data, } if (any(hill_no_outlets[,2] == 0) & length(flw_struct[,1]) != 1) { # if there are any hillslopes without streams - cat("Correcting for hillslopes missing stream outlets") + cat("Correcting for hillslopes missing stream outlets\n") streams = flw_struct[flw_struct$Landtype == 1,] # make var of streams for (i in hill_no_outlets[hill_no_outlets[,2] == 0,1]) { hill_patches = flw_struct[flw_struct$Hill == i,] # get patches from problem hillslope @@ -250,16 +252,18 @@ make_flow_list <- function(raw_patch_data, } if (!is.null(no_stream_fix)) { cat(paste(length(no_stream_fix),"hillslopes had their lowest elevation patches set to streams, at a max distance from existing streams of", - max(print_dist_fix),"cell lengths.")) + max(print_dist_fix),"cell lengths.\n")) } if (!is.null(no_stream)) { # output hillslopes that weren't corrected - cat("The following outlet patches were outside of the",make_stream,"cell length distance threshold set by the 'make_stream' argument.") + cat("The following outlet patches were outside of the",make_stream,"cell length distance threshold set by the 'make_stream' argument.\n") cat("Dist2Stream" = print_dist,flw_struct[no_stream,]) stop(noquote("The above hillslopes must have stream patch outlets, either increase the value of the 'make_stream' argument, or fix via GIS.")) } find_connected = function(patch_borders, start, history = NULL, found = NULL, missing, hill, itr_ct = NULL, itr_max = 10000) { + while(length(missing) > 0) { + # add start to found, only for first itr tho if (is.null(history)) { found = start @@ -269,6 +273,9 @@ make_flow_list <- function(raw_patch_data, # 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) + + #cat("neighbors:",neighbors, "missing(size):", length(missing),"\n") + # update the found vector to include those neighbors found = c(found, neighbors) # remove the found from missing @@ -284,8 +291,13 @@ make_flow_list <- function(raw_patch_data, if (!is.null(itr_ct)) { itr_ct = itr_ct + 1 + + #cat("itr_ct:",itr_ct," length_found:", length(found), + # " length_history:",length(history), "length_miss:",length(missing), + # "\n") + if (itr_ct == itr_max) { - cat("Reached iteration limit of ",itr_max) + cat("Reached iteration limit of ",itr_max, "\n") return(list(missing, found, history, start)) } } @@ -293,23 +305,25 @@ make_flow_list <- function(raw_patch_data, # if we haven't found everything, keep iterating through found patches not in the history 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, - itr_ct = itr_ct, - itr_max = itr_max - ) + start = next_patch + #next_patch = as.numeric(found[!found %in% history][1]) + } #while + #find_connected( + # patch_borders = patch_borders, + # start = next_patch, + # history = history, + # found = found, + # missing = missing, + # hill = hill, + # itr_ct = itr_ct, + # itr_max = itr_max + #) } if (!skip_hillslope_check) { # check for segmented hillslopes - this is going to be so slow but oh well - cat("Checking for segmented hillslopes") + cat("Checking for segmented hillslopes\n") segmented_hills = NULL segmented_patch_ct = NULL segmented_patch_list = list() @@ -318,7 +332,7 @@ make_flow_list <- function(raw_patch_data, ct = 0 # for testing: - i = unique(flw_struct$Hill)[4] + #i = unique(flw_struct$Hill)[4] for (i in unique(flw_struct$Hill)) { ct = ct + 1 @@ -337,18 +351,20 @@ make_flow_list <- function(raw_patch_data, itr_ct = 0, itr_max = length(all_patches) ) - 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 + options(warn=1) 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" ) + #cat("\nHill:",i," num_total_patches:",length(all_patches),"\n") print(segmented_df) + #error("") } @@ -364,7 +380,7 @@ make_flow_list <- function(raw_patch_data, # Build list for output. Turn border count into probabilities and lists of neighbors lst <- list() - cat("Buildling flowtable list") + cat("Buildling flowtable list\n") pb = txtProgressBar(min = 0,max = length(flw_struct$Number),style = 3) for (i in 1:length(flw_struct$Number)) { @@ -404,12 +420,19 @@ make_flow_list <- function(raw_patch_data, if (gamma_tot != 0) { # if there's a downslope neighbor tp_gamma <- tp_gamma/gamma_tot # normalize gamma by total (% or proportion) - tp_TotalG <- (gamma_tot/perim_sum)*flw_struct$Area[i] # gamma_tot/perim_sum = sum of slopes * area = volume + #09212022LML the horizental flow (from each reservoir layer) should be calculated as layer depth X Ksat X perim X slope (i.e dy/dx) + #the calculation of transimissivity profile should be revised accordingly in RHESSys + #tp_TotalG <- (gamma_tot/perim_sum)*flw_struct$Area[i] # gamma_tot/perim_sum = sum of slopes * area = volume + tp_TotalG <- gamma_tot # sum perim * slope of all neighbors } else { if (is.null(slope_i) | length(slope_i) == 0) { # if there is only one patch slope_i will be null, - tp_TotalG <- flw_struct$Area[i] + #09212022LML + #tp_TotalG <- flw_struct$Area[i] + tp_TotalG <- 0.0 } else {#if all upslope, take slope from closest neighbor in height - tp_TotalG <- -max(slope_i)*flw_struct$Area[i] + #09212022LML + #tp_TotalG <- -max(slope_i)*flw_struct$Area[i] + tp_TotalG <- -max(slope_i)*perim_sum } }