Skip to content

Commit 6f9c173

Browse files
wburke24wburke24
committed
Added update_world.R and spatial_input_gen.R
Added update_world.R to modify existing worldfiles, and spatial_input_gen.R to generate rhessys input maps with grass7. More documentation and support for both still to come. More changes to vignette
1 parent 1ec42e5 commit 6f9c173

4 files changed

Lines changed: 358 additions & 1 deletion

File tree

‎R/spatial_input_gen.R‎

Lines changed: 151 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,151 @@
1+
#' spatial_input_gen
2+
#'
3+
#' Generate the spatial inputs for RHESSys and RHESSysPreprocessing
4+
#' @param Name Name of your project
5+
#' @param basin The area defining your basin/watershed of interest. Assumes it's a raster
6+
#' @author Will Burke
7+
#' @export
8+
9+
10+
spatial_input_gen = function(name,
11+
DEM,
12+
GRASS_path,
13+
threshold,
14+
basin = NULL,
15+
patch = NULL,
16+
define_watershed = NULL,
17+
easting = NULL,
18+
northing = NULL,
19+
gisDbase = NULL,
20+
location = NULL,
21+
mapset = NULL) {
22+
23+
24+
# check and convert file paths for args
25+
26+
27+
# check GRASS location and version - requires grass 7
28+
29+
30+
out_maps = NULL
31+
32+
# set initial grass environment
33+
init = rgrass7::initGRASS(gisBase = GRASS_path,
34+
home = tempdir(),
35+
override = TRUE)
36+
# add basin map to new location
37+
rgrass7::execGRASS(cmd = "r.in.gdal", input = basin, output = "basin", location = "temp_loc")
38+
# set grass environment to new (correctly projected) location
39+
init = rgrass7::initGRASS(gisBase = GRASS_path,
40+
home = tempdir(),
41+
location = "temp_loc",
42+
mapset = "PERMANENT",
43+
override = TRUE)
44+
45+
# check if same proj - should be quiet but idk
46+
check = rgrass7::execGRASS(cmd = "r.in.gdal",flags = c("j","quiet"), input = DEM)
47+
# if projections don't match
48+
if (check == 1) {
49+
# add DEM to its own location
50+
rgrass7::execGRASS(cmd = "r.in.gdal", input = DEM, output = "DEM", location = "dem_loc")
51+
# project DEM from that location to current loc (temp_loc)
52+
rgrass7::execGRASS(cmd = "r.proj", location = "dem_loc", mapset = "PERMANENT", input = "DEM")
53+
}
54+
55+
# --- checking env and stuff ---
56+
# gmeta()
57+
# rgrass7::execGRASS(cmd = "g.mapset", flags = "p")
58+
# rgrass7::execGRASS(cmd = "g.mapset", flags = "l")
59+
# rgrass7::execGRASS(cmd = "g.list", flags = "p", type = "raster")
60+
61+
# check if basin exists and use it as region if it does, if not use raster
62+
rgrass7::execGRASS(cmd = "g.region", raster = "basin")
63+
#rgrass7::execGRASS(cmd = "g.region", flags = "p") # print it
64+
65+
# get slope and aspect maps
66+
rgrass7::execGRASS(cmd = "r.slope.aspect", elevation = "DEM", slope = "slope", aspect = "aspect")
67+
out_maps = c(out_maps, "slope", "aspect")
68+
69+
# get horizon maps, in degrees
70+
rgrass7::execGRASS(cmd = "r.horizon", flags = "d", elevation = "DEM", direction = 0, output = "east")
71+
rgrass7::execGRASS(cmd = "r.horizon", flags = "d", elevation = "DEM", direction = 180, output = "west")
72+
out_maps = c(out_maps, "east_000", "west_180")
73+
74+
rgrass7::execGRASS(cmd = "r.mapcalc", expression = 'east_horizon = sin(east_000)')
75+
rgrass7::execGRASS(cmd = "r.mapcalc", expression = 'west_horizon = sin(west_180)')
76+
out_maps = c(out_maps, "east_horizon", "west_horizon")
77+
78+
#SET TO BE OPTIONAL
79+
#rgrass7::execGRASS(cmd = "r.mapcalc", expression = 'ehr.100 = east_horizon*100')
80+
#rgrass7::execGRASS(cmd = "r.mapcalc", expression = 'whr.100 = west_horizon*100')
81+
#out_maps = c(out_maps, "ehr.100", "whr.100")
82+
83+
# watershed analysis
84+
subbasin_name = paste0("subbasin.t",as.character(threshold))
85+
stream_name = paste0("stream.t",as.character(threshold))
86+
hillslope_name = paste0("hillslope.t",as.character(threshold))
87+
88+
rgrass7::execGRASS(cmd = "r.watershed", elevation = "DEM", threshold = threshold, accumulation = "accumulation",
89+
drainage = "drainage", basin = subbasin_name, stream = stream_name, half_basin = hillslope_name)
90+
out_maps = c(out_maps, subbasin_name, stream_name, hillslope_name)
91+
92+
# ----- patch map -----
93+
if (patch == "topo") {
94+
# Topographical
95+
hinfo = rgrass7::execGRASS(cmd = "r.info", map = hillslope_name, intern = TRUE)
96+
minmax_ch = unlist(strsplit(grep("Range of data: \\s+ min = \\d\\s+max = \\d",hinfo, value = TRUE),split = "\\s+"))
97+
hmax = as.numeric(minmax_ch[which(minmax_ch == "max") + 2])
98+
99+
#rgrass7::execGRASS(cmd = "r.info", map = "DEM")
100+
rgrass7::execGRASS(cmd = "r.mapcalc", expression = paste0('p.topo = (DEM * ', hmax, ') + ', hillslope_name))
101+
rgrass7::execGRASS(cmd = "r.clump", input = "p.topo", output = "p.topo.cl")
102+
out_maps = c(out_maps, "p.topo.cl")
103+
104+
} else if (patch == "grid") {
105+
# Grid
106+
rgrass7::execGRASS(cmd = "r.random.cells", output = "grid", distance = 0)
107+
rgrass7::execGRASS(cmd = "r.mapcalc", expression = paste0('p.grid = (grid * ', hmax, ') + ', hillslope_name))
108+
rgrass7::execGRASS(cmd = "r.clump", input = "p.grid", output = "p.grid.cl")
109+
out_maps = c(out_maps, "p.grid.cl")
110+
111+
} else if (patch == "variable") {
112+
113+
# THIS IS MOSTLY HARDCODED AND PROBABLY SHOULD BE ADDED TO A R-BASED FUNCTION THAT CAN DO IT WITH DYNAMIC RESOLUTION INPUTS
114+
# Variable resolution
115+
rgrass7::execGRASS(cmd = "g.region", raster = "DEM", res = 90)
116+
rgrass7::execGRASS(cmd = "r.mapcalc", expression = 'DEM90m = DEM')
117+
rgrass7::execGRASS(cmd = "g.region", raster = "DEM") #(reset the region back to 30m)
118+
rgrass7::execGRASS(cmd = "r.topidx", input = "DEM", output = "lna")
119+
rgrass7::execGRASS(cmd = "r.fillnulls", input = "lna", output = "lna.fill")
120+
rgrass7::execGRASS(cmd = "r.mapcalc", expression = 'lna.int = round(lna.fill)')
121+
rgrass7::execGRASS(cmd = "r.mapcalc", expression = 'riparian = lna.int > 7.5')
122+
rgrass7::execGRASS(cmd = "r.mapcalc", expression = paste0('p.rip30.up90 = (riparian == 1) * (DEM *', hmax,') + (riparian < 1)* (DEM90m * 633) + ', hillslope_name))
123+
rgrass7::execGRASS(cmd = "r.clump", input = "p.rip30.up90", output = "p.30_90res.cl")
124+
out_maps = c(out_maps, "p.30_90res.cl")
125+
126+
} else if (patch == "aspatial_wetness") {
127+
# ALSO DIDNT TEST THIS ONE
128+
rgrass7::execGRASS(cmd = "r.mapcalc", expression = paste0('p.lna = (lna.int * ', hmax,') + ', hillslope_name))
129+
out_maps = c(out_maps, "p.lna")
130+
131+
}
132+
133+
# Roads - zero map
134+
rgrass7::execGRASS(cmd = "r.mapcalc", expression = 'zero = ("DEM" > 0)*0')
135+
136+
if (define_watershed) {
137+
# ADDING AN INTERACTIVE GUI MAP HERE WOULD ALLOW FOR PRETTY SIMPLE BASIN DELINEATION
138+
rgrass7::execGRASS(cmd = "r.water.outlet", input = "drainage", output = "basin",coordinates = c(easting, northing))
139+
rgrass7::execGRASS(cmd = "d.rast", "basin")
140+
rgrass7::execGRASS(cmd = "g.copy", raster = "basin", 'MASK')
141+
}
142+
143+
# Write rasters
144+
out_names = paste0(name,"_",out_maps,".tif")
145+
for (i in 1:length(out_maps)) {
146+
rgrass7::execGRASS(cmd = "r.out.gdal", input = out_maps[i], output = out_names[i], format = "GTiff", type = "Float64")
147+
}
148+
149+
150+
}
151+

‎R/update_world.R‎

Lines changed: 96 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,96 @@
1+
#' update_world
2+
#'
3+
#' Read and modify an existing worldfile. Optimized somewhat (for R) but awk may still be faster.
4+
#' @param worldfile Name and path of worldfile to be created.
5+
#' @param vars State variable(s) from the worldfile to be modified
6+
#' @param values Values to replace for `vars`, or coefficients to multiply by. N args must equal that of vars or be 1 (will be used for all vars).
7+
#' or can be a list of values where number of list elements = n vars
8+
#' @param level_names The level for which a given state variable will have it's value changed or replaced. Can be a single character string or vector
9+
#' of character strings. Valid arguments are: world, basin, hillslope, zone, patch, canopy_strata.
10+
#' @param level_IDs The IDs of the correstponding `level_names` for which a given state variable will have it's value changed or replaced. Set to NULL or "all"
11+
#' to change a value across all of a given level (e.g. for all patches). This can be a vector same length as names, or a list, with IDs for each named level in names
12+
#' @param out_file Destination file to write new worldfile
13+
#' @param overwrite TRUE/FALSE if input worldfile should be overwritten
14+
#' @author Will Burke
15+
#' @export
16+
17+
update_world = function(worldfile, out_file, vars, values,level_names = NULL, level_IDs = NULL, overwrite = FALSE) {
18+
19+
# NOTES
20+
# - vdouble check all the list options work and check var lengths/numbers when needed
21+
22+
# ---------- Check Aguments ----------
23+
if (is.null(out_file) & !overwrite) {stop(noquote("No destination file set by 'out_file' and 'overwrite' is FALSE"))}
24+
if ( file.exists(out_file) & overwrite == FALSE) {stop(noquote(paste0("File '",out_file,"' already exists and 'overwrite' argument is FALSE")))}
25+
if (length(vars) != length(values) & length(values) != 1 & !is.list(values))
26+
27+
if (!is.list(values)) {
28+
if (!is.character(values)) {
29+
values = as.character(values)
30+
}
31+
}
32+
33+
# for (i in length(vars)) {
34+
# if (level_names[i] == "all") {level_names[i] = c("world", "basin", "hillslope", "zone", "patch", "canopy_strata")}
35+
# }
36+
#
37+
# if (level_IDs == "all") {}
38+
39+
# ---------- Parse Worldfile ----------
40+
options(scipen = 999)
41+
# parsing the values as characters to retain the exact value/precision
42+
read_world = readLines(worldfile)
43+
world = strsplit(trimws(read_world), "\t+")
44+
world = data.frame(matrix(unlist(world), nrow=length(world), byrow=T),stringsAsFactors = FALSE)
45+
names(world) = c("values","vars")
46+
47+
# ---------- Find Levels----------
48+
index_all = which(world[,2] == "world_ID" | world[,2] == "basin_ID" | world[,2] == "hillslope_ID" |
49+
world[,2] == "zone_ID" | world[,2] == "patch_ID" | world[,2] == "canopy_strata_ID")
50+
index_names = gsub("_ID", "", x = world$vars[index_all])
51+
index_max = c(index_all[2:length(index_all)]-1, length(read_world))
52+
world$level = unname(unlist(mapply(rep,index_names, (index_max - index_all) + 1 )))
53+
world$ID = unname(unlist(mapply(rep, world$values[index_all], (index_max - index_all) + 1 )))
54+
55+
# ---------- Find and Replace Vars ----------
56+
for (i in length(vars)){
57+
58+
if (!is.null(level_names[i])) {
59+
find_index = world$level == level_names[i]
60+
} else {
61+
find_index = rep(TRUE,length(world$level))
62+
}
63+
if (!is.null(level_IDs[i])) {
64+
find_index = find_index & world$ID == level_IDs[i]
65+
}
66+
replace_index = which(world$vars == vars[i] & find_index)
67+
68+
# if unique values for every instance of var to be replaces were given, do nothing, otherwise repeat to get enough replacement values
69+
current_value = world$values[replace_index]
70+
if (length(values[i]) != length(replace_index)) {
71+
new_value = rep(values[i], length(replace_index)/length(values[i]))
72+
} else {
73+
new_value = values[i]
74+
}
75+
76+
if (any(startsWith(new_value,"*"))) {
77+
new_value = as.numeric(trimws(substr(new_value[startsWith(new_value,"*")],2,nchar(new_value[startsWith(new_value,"*")]))))
78+
new_value = new_value * as.numeric(current_value)
79+
} else {
80+
new_value = as.numeric(trimws(new_value))
81+
}
82+
83+
# generic sub/gsub
84+
#sub("\\d+[[:blank:]]",XXX"\t", read_world)
85+
86+
read_world[replace_index] = unname(mapply(sub,paste0(current_value,"[[:blank:]]"),paste0(new_value,"\t"),read_world[replace_index]))
87+
88+
}
89+
90+
# ---------- Write file ----------
91+
writeLines(text = read_world,out_file)
92+
93+
print(noquote(paste("Successfully rote updated worldfile to",out_file)))
94+
95+
96+
}

‎Script_Run_RHESSysPreprocess.R‎

Lines changed: 110 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,110 @@
1+
# run_RHESSysPreprocess
2+
# Will Burke 3/5/18
3+
# Last updated 4/20/19
4+
5+
# Instructions
6+
# ------------
7+
# This is an example script showing how the RHESSysPreprocess.R function should be run.
8+
# 1) Install or source the RHESSysPreprocessing package
9+
# 2) Copy this script, and edit where indicated.
10+
# 3) Run the RHESSysPreprocess.R function at the bottom.
11+
# 4) The funciton will produce:
12+
# - worldfile
13+
# - flowtable
14+
# - optionally: header
15+
# (metadata output is being worked on currently)
16+
17+
# Load Package
18+
# ------------
19+
# Assuming you installed the RHESSysPreprocessing package, you will need to load it
20+
library(RHESSysPreprocessing)
21+
22+
# Filepaths
23+
# ----------------
24+
# This script uses relative filepaths. This means that it will look folders and files relative to your current working directory.
25+
# If needed, set your current working directory to the folder of your project:
26+
setwd("~/Documents/MyProject")
27+
# This script also uses the "~", which is a shorthand method of navigating to your "home" user directory - typically the folder named for your username.
28+
29+
# Spatial Data
30+
# ------------
31+
# You will need to select your method of geospatial data input. This is the means by which the spatial data referenced in your template.
32+
33+
# Currently there are a two methods:
34+
# 1) Raster - spatial data in any raster format supported by R GDAL will be read in from a folder.
35+
# 2) < NO LONGER BEING ACTIVELY SUPPORTED> GRASS GIS > - GRASS 6 or 7, spatial data will be imported from the specified GRASS location and mapset.
36+
37+
# NOTES:
38+
# - Due to a variety of factors, spatial data import via the raster method is both more robust and faster.
39+
# - Regardless of import method, good practice for spatial data should be followed.
40+
# - Input data should have the same projections, extents, and cell sizes. THIS MAY RESULT IN ERRORS IF NOT FOLLOWED.
41+
42+
# Raster
43+
# ------
44+
# To import spatial data from a folder of rasters:
45+
# 1) Set type to "raster"
46+
type = "raster"
47+
# 2) Set typepars to the path of the folder containing your rasters
48+
typepars = "spatial_data"
49+
50+
51+
# Template
52+
# --------
53+
# The worldfile template is the key document that outlines how your worldfile will be built.
54+
# The template variable should point to the name and location of your template.
55+
template = "/templates/example.template"
56+
57+
58+
# Name
59+
# ----
60+
# Set the name and path for all function outputs.
61+
# Suffixes of .world, .flow, and .meta will be appended to the worldfile, flowtable, and metadata files respectively.
62+
name = "/output/my_watershed"
63+
64+
# Overwrite
65+
# ---------
66+
# TRUE/FALSE if an existing worldfile and flowtable be overwritten.
67+
overwrite = FALSE
68+
69+
# Streams
70+
# -------
71+
# Streams map to be used in creation of the flowtable - this is just the name of the map, to be found via the method indicated with "type"
72+
streams = "my_watershed_streams"
73+
74+
# Optional Flowtable Spatial Data
75+
# -------------------------------
76+
# These maps are optional inputs in flowtable creation
77+
# roads = "roads_map"
78+
# impervious = "impervious_map"
79+
# roofs = "roofs_map"
80+
81+
# Header
82+
# ------
83+
# TRUE/FALSE to produce a header file. Header file will be have same name(and location) set by "name", with the ".hdr" suffix.
84+
header = FALSE
85+
86+
# Parallelization
87+
# ---------------
88+
# Current (Dec 2018 and on) develop branch RHESSys is hillslope paralleized and requires a flowtable that is compatible.
89+
# Parallelization must be set to TRUE (this is a default)
90+
# parallel = TRUE
91+
92+
# The "make_stream" argument defines the distance from an existing stream that the outlet of a hillslope can be set to be a
93+
# stream.Since all hilslopes must have stream outlets, if a hillslope outlet exists outside of the distance threshold set by
94+
# "make_stream",an error will occur and indicae the problem hillslope/outlet patch of that hillslope. This typically occurs
95+
# as an artifact of how watershed analysis is done, and hillslopes are created, which sometimes results in fragmented or very
96+
# small/skinny hillslopes, far away from streams.
97+
# "make_stream" can be set to any positive value, or TRUE to always set hillslope outlets to streams.
98+
# Default is 4, which is meant to roughly account for the errors/aritifacts that might occur from GIS, without including any
99+
# extreme outlying hillslopes
100+
# make_stream = 4
101+
102+
# Finally, run the function. Depending on size, it may take a minute or two.
103+
RHESSysPreprocess(
104+
template = template,
105+
name = name,
106+
type = type,
107+
typepars = typepars,
108+
streams = streams,
109+
overwrite = overwrite,
110+
header = header)

vignettes/RHESSysPreprocess.Rmd renamed to vignettes/Run_RHESSysPreprocess.Rmd

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,5 @@
11
---
2-
title: "RHESSysPreprocess"
2+
title: "Run_RHESSysPreprocess"
33
author: "William Burke"
44
date: "`lir Sys.Date()`"
55
output: rmarkdown::html_vignette

0 commit comments

Comments
 (0)