From 6e0acf02c17c0fa0e9b4c4769bb7ae49ca7b63f1 Mon Sep 17 00:00:00 2001 From: liuminggee <99829902+liuminggee@users.noreply.github.com> Date: Wed, 21 Sep 2022 11:26:08 -0700 Subject: [PATCH 1/7] Update make_flow_list.R --- R/make_flow_list.R | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/R/make_flow_list.R b/R/make_flow_list.R index 57c43ce..f27e5ef 100644 --- a/R/make_flow_list.R +++ b/R/make_flow_list.R @@ -404,12 +404,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 } } From b47dd59f6e7385261f2df8b7c7632181111c9607 Mon Sep 17 00:00:00 2001 From: Mingliang Liu Date: Tue, 27 Sep 2022 00:22:29 -0700 Subject: [PATCH 2/7] replace recursive find_connected() to iteration --- R/RHESSysPreprocess.R | 16 ++++++------- R/create_flownet.R | 5 ++-- R/make_flow_list.R | 53 ++++++++++++++++++++++++++----------------- 3 files changed, 43 insertions(+), 31 deletions(-) diff --git a/R/RHESSysPreprocess.R b/R/RHESSysPreprocess.R index f7fd51f..a52e776 100644 --- a/R/RHESSysPreprocess.R +++ b/R/RHESSysPreprocess.R @@ -113,14 +113,14 @@ RHESSysPreprocess = function(template, } } - world_gen_out = world_gen(template = template, - worldfile = worldfile, - type = type, - typepars = typepars, - overwrite = overwrite, - header = header, - unique_strata_ID = unique_strata_ID, - asprules = asprules) + #world_gen_out = world_gen(template = template, + # worldfile = worldfile, + # type = type, + # typepars = typepars, + # overwrite = overwrite, + # header = header, + # unique_strata_ID = unique_strata_ID, + # asprules = asprules) #readin = world_gen_out[[1]] #asp_rules = world_gen_out[[2]] 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/make_flow_list.R b/R/make_flow_list.R index f27e5ef..2c7f150 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)]), @@ -128,7 +128,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 +225,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 +250,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 +271,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 +289,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 +303,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 +330,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,7 +349,6 @@ 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)) @@ -364,7 +375,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)) { From 6e37d89b438ddcdb76576a9a88039a476d214927 Mon Sep 17 00:00:00 2001 From: Mingliang Liu Date: Wed, 28 Sep 2022 16:51:18 -0700 Subject: [PATCH 3/7] test --- R/RHESSysPreprocess.R | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/R/RHESSysPreprocess.R b/R/RHESSysPreprocess.R index a52e776..f7fd51f 100644 --- a/R/RHESSysPreprocess.R +++ b/R/RHESSysPreprocess.R @@ -113,14 +113,14 @@ RHESSysPreprocess = function(template, } } - #world_gen_out = world_gen(template = template, - # worldfile = worldfile, - # type = type, - # typepars = typepars, - # overwrite = overwrite, - # header = header, - # unique_strata_ID = unique_strata_ID, - # asprules = asprules) + world_gen_out = world_gen(template = template, + worldfile = worldfile, + type = type, + typepars = typepars, + overwrite = overwrite, + header = header, + unique_strata_ID = unique_strata_ID, + asprules = asprules) #readin = world_gen_out[[1]] #asp_rules = world_gen_out[[2]] From d1084e1dd0754814bf8d58c6e710606db107b01f Mon Sep 17 00:00:00 2001 From: Mingliang Liu Date: Wed, 28 Sep 2022 17:07:18 -0700 Subject: [PATCH 4/7] add road_width argument into RHESSysPreprocess function --- R/RHESSysPreprocess.R | 2 ++ 1 file changed, 2 insertions(+) diff --git a/R/RHESSysPreprocess.R b/R/RHESSysPreprocess.R index f7fd51f..54e01db 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, @@ -145,6 +146,7 @@ RHESSysPreprocess = function(template, streams = streams, overwrite = overwrite, roads = roads, + road_width = road_width, impervious = impervious, roofs = roofs, wrapper = wrapper, From 058d0ec4d6ca3ff92f3b45ce1780275b2bacdc30 Mon Sep 17 00:00:00 2001 From: Mingliang Liu Date: Wed, 28 Sep 2022 23:08:03 -0700 Subject: [PATCH 5/7] fix a bug in filling pits (find stream) --- R/GIS_read.R | 5 +-- R/RHESSysPreprocess.R | 27 ++++++++------ R/find_stream.R | 82 ++++++++++++++++++++++--------------------- 3 files changed, 61 insertions(+), 53 deletions(-) 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 54e01db..2ed725f 100644 --- a/R/RHESSysPreprocess.R +++ b/R/RHESSysPreprocess.R @@ -79,8 +79,8 @@ 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( @@ -104,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( "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]] @@ -129,22 +131,24 @@ 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( "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, @@ -153,6 +157,7 @@ RHESSysPreprocess = function(template, parallel = parallel, make_stream = make_stream, skip_hillslope_check = skip_hillslope_check) + } # ---------- Run build_meta ---------- # if (meta) { 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) +} From 33f44977134f946421a83027290ceb16a0f38011 Mon Sep 17 00:00:00 2001 From: Mingliang Liu Date: Thu, 13 Oct 2022 15:51:14 -0700 Subject: [PATCH 6/7] print right place for warnning message --- R/RHESSysPreprocess.R | 6 +++--- R/make_flow_list.R | 3 +++ 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/R/RHESSysPreprocess.R b/R/RHESSysPreprocess.R index 2ed725f..68f9436 100644 --- a/R/RHESSysPreprocess.R +++ b/R/RHESSysPreprocess.R @@ -84,7 +84,7 @@ RHESSysPreprocess = function(template, 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) @@ -106,7 +106,7 @@ RHESSysPreprocess = function(template, 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) { @@ -133,7 +133,7 @@ RHESSysPreprocess = function(template, 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) { diff --git a/R/make_flow_list.R b/R/make_flow_list.R index 2c7f150..097ba5d 100644 --- a/R/make_flow_list.R +++ b/R/make_flow_list.R @@ -354,12 +354,15 @@ make_flow_list <- function(raw_patch_data, 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("") } From dc6055a41079047d3b485e5448529dd4e18495d3 Mon Sep 17 00:00:00 2001 From: Mingliang Liu Date: Wed, 21 Dec 2022 16:56:36 -0800 Subject: [PATCH 7/7] fix the bug when the tif raster shows a patch is a stream AND a road (set stream instead) --- R/make_flow_list.R | 2 ++ 1 file changed, 2 insertions(+) diff --git a/R/make_flow_list.R b/R/make_flow_list.R index 097ba5d..6756235 100644 --- a/R/make_flow_list.R +++ b/R/make_flow_list.R @@ -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) --------------------