diff --git a/R/reclass_raster.R b/R/reclass_raster.R index 947ae14d..0dfb0df6 100644 --- a/R/reclass_raster.R +++ b/R/reclass_raster.R @@ -16,6 +16,9 @@ #' @param new_val character. The name of the column with the new raster values, #' which need to be integer values. In case of floating point values, consider #' multiplying the values e.g. by 1000 to keep three decimals. +#' @param reclass_value integer. Value to reclassify the entire raster. +#' Default is FALSE. Note that in case reclass_value and new_val +#' is provided the raster is reclassified on reclass_value. #' @param raster_layer Full path to the input raster .tif layer. #' @param recl_layer character. Full path of the output .tif layer, i.e., the #' reclassified raster file. @@ -43,7 +46,7 @@ #' @importFrom terra rast #' @export #' -#' @author Marlene Schürz +#' @author Marlene Schürz, Thomas Tomiczek #' #' @references #' https://grass.osgeo.org/grass82/manuals/r.reclass.html @@ -56,9 +59,7 @@ #' download_test_data(my_directory) #' #' # Read the stream order for each sub-catchment as a data.table -#' my_dt <- read_geopackage(gpkg= paste0(my_directory, -#' "/hydrography90m_test_data", -#' "/order_vect_59.gpkg"), +#' my_dt <- read_geopackage(paste0(my_directory, "/hydrography90m_test_data/order_vect_59.gpkg"), #' import_as = "data.table") #' #' @@ -79,10 +80,15 @@ #' raster_layer = stream_raster, #' recl_layer = recl_raster) #' +# Reclassify the stream network with the value 1 across the network +#' str_ord_rast <- reclass_raster(data = str_ord$stream, +#' reclass_value = 1, +#' rast_val = "stream", +#' raster_layer = stream_raster, +#' recl_layer = recl_raster) - -reclass_raster <- function(data, rast_val, new_val, raster_layer, - recl_layer, +reclass_raster <- function(data, rast_val, new_val = FALSE, raster_layer, + recl_layer, reclass_value = FALSE, no_data = -9999, type = "Int32", compression = "low", bigtiff = TRUE, read = FALSE, quiet = TRUE) { @@ -92,33 +98,44 @@ reclass_raster <- function(data, rast_val, new_val, raster_layer, if (missing(data)) stop("data: Input data.frame is missing.") - # Check if input data is of type data.frame, - # data.table or tibble + # Check if input data is of type data.frame, data.table or tibble if (!is(data, "data.frame")) stop("data: Has to be of class 'data.frame'.") - # Check if rast_val and new_val is defined + # Check if rast_val is defined if (missing(rast_val)) stop("rast_val: Column name of current raster value is missing.") - if (missing(new_val)) - stop("new_val: Column name of new raster value is missing.") - # Check if rast_val/new_val column names exist + # Check if rast_val column names exist if (is.null(data[[rast_val]])) stop(paste0("rast_val: Column name '", rast_val, - "' does not exist.")) - if (is.null(data[[new_val]])) - stop(paste0("new_val: Column name '", new_val, - "' does not exist.")) + "' does not exist.")) - # Check if values of the rast_val/new_val columns are numeric + # Check if values of the rast_val columns are numeric if (!is.integer(data[[rast_val]])) stop( paste0("rast_val: Values of column ", rast_val, - " have to be integers.")) - if (!is.integer(data[[new_val]])) - stop(paste0("new_val: Values of column ", new_val, - " have to be integers.")) + " have to be integers.")) + + # Check if new_val column names exist when no reclass_values is given + if (isFALSE(new_val) && isFALSE(reclass_value)) + stop(paste0("new_val: Column name '", new_val, + "' does not exist.")) + + # Check if values of the new_val columns are numeric when no reclass_values is given + if (isFALSE(reclass_value)) { + if (!is.integer(data[[new_val]])) { + stop(paste0("reclass_value:", reclass_value, " must be integers.")) + } + + } + + # Check if reclass_value is an numeric or integer value + if (!isFALSE(reclass_value)) { + if (!(is.numeric(reclass_value) || is.integer(reclass_value))) { + stop(paste0("reclass_value:", reclass_value, " must be integers.")) + } + } # Check if raster_layer is defined if (missing(raster_layer)) @@ -131,8 +148,8 @@ reclass_raster <- function(data, rast_val, new_val, raster_layer, # Check if raster_layer ends and recl_layer with .tif if (!endsWith(raster_layer, ".tif")) stop("raster_layer: Input raster is not a .tif file.") - if (!endsWith(recl_layer, ".tif")) - stop("recl_layer: Output raster file path needs to end with .tif.") + if (!endsWith(recl_layer, ".tif")) + stop("recl_layer: Output raster file path needs to end with .tif.") # Check if recl_layer is defined if (missing(recl_layer)) @@ -145,7 +162,6 @@ reclass_raster <- function(data, rast_val, new_val, raster_layer, stop("type: Has to be 'Byte', 'Int16', 'UInt16', 'Int32', 'UInt32', 'CInt16', or 'CInt32' ") - # Check and translate compression into the compression type and the # compression level which is applied to the tiff file when writing it. if(compression == "none") { @@ -170,7 +186,7 @@ reclass_raster <- function(data, rast_val, new_val, raster_layer, } if (!is.logical(read)) - stop("read: Has to be TRUE or FALSE.") + stop("read: Has to be TRUE or FALSE.") # Check if quiet is logical if (!is.logical(quiet)) @@ -179,12 +195,100 @@ reclass_raster <- function(data, rast_val, new_val, raster_layer, # Make bash scripts executable make_sh_exec() - # The r.reclass function of GRASS GIS requires a text file - # including the old and the new value with an = between - # (e.g. 1 = 20) - rules <- data.table(old = data[[rast_val]], - equal = "=", - new = data[[new_val]]) + # To check the raster values and data values + # load raster and convert to data.frame + rast_dat <- rast(raster_layer) + rast_dat <- as.data.frame(rast_dat) + rast_dat <- as.data.frame(unique(rast_dat[[1]])) + colnames(rast_dat) <- "val" + + # Check if rast_val is missing raster values + if (length(data[[rast_val]]) < length(rast_dat[[1]])) { + print("Reclassification is missing raster values: Warning NA's are introduced!") + } + # Check and handle if raster values are provided in rast_val that are not in the raster tif file. + if (length(data[[rast_val]]) != length(rast_dat[[1]])) { + # Index all raster values which are not in the input data table + indx_miss_raster <- which(rast_dat[[1]] %in% data[[rast_val]]) + # Get missing raster values + miss_raster <- rast_dat[-c(indx_miss_raster),] + print(paste0("These values of the raster were not found in the data table:", + paste(miss_raster, collapse = ", "))) + # Write all values found in raster tif file and input data table as data frame + same_val1 <- as.data.frame(rast_dat[c(indx_miss_raster),]) + # Set name for raster values + colnames(same_val1) <- "val" + # Index all input data table values which are not in the raster tif file + indx_miss_rast_val <- which(data[[rast_val]] %in% rast_dat[[1]]) + # Get missing input data values + miss_rast_val <- data[-c(indx_miss_rast_val),] + # Get only missing raster input data values to throw out message + missing_rast_values <- data[-c(indx_miss_rast_val),1] + print(paste0("These raster values of the data table were not found in the raster:", + paste(missing_rast_values, collapse = ", "))) + # Write all values found in input data table and raster file as data frame + same_val2 <- data[c(indx_miss_raster),] + # Combine all values found in raster and input data table + same_val <- as.data.frame(cbind(same_val1, same_val2)) + # Write missing raster values as data frame + miss_raster <- as.data.frame(miss_raster) + # Set name for raster values + colnames(miss_raster) <- "val" + # Give missing values NA + miss_raster$rast_val <- NA + miss_raster$new_val <- NA + # Set to the same names to combine with all values table + colnames(miss_raster) <- c("val", rast_val, new_val) + # Combine all values table with missing values table and use as input data + data <- rbind(same_val, miss_raster) + + # List both rast_val and new_val columns to check if they are of equal length + dat <- list(data[[rast_val]], list(data[[new_val]])) + + # In case new_val is bigger than rast_val length of rast_val will be used to reclassify + data <- setNames(do.call(cbind.data.frame, + lapply(lapply(dat, unlist), + `length<-`, max(lengths(dat)))), paste0(c(rast_val, new_val))) + if (isFALSE(reclass_value)) { + # The r.reclass function of GRASS GIS requires a text file + # including the old and the new value with an = between + # (e.g. 1 = 20) + rules <- data.table::data.table(old = data[[rast_val]], + equal = "=", + new = data[[new_val]]) + } + } + + if (!isFALSE(reclass_value)) { + + # use reclass_value for reclassification + data$reclass <- reclass_value + data$reclass <- as.integer(data$reclass) + # + # + # + # + # + # rand_string <- stri_rand_strings(n = 1, length = 8, pattern = "[A-Za-z0-9]") + # rules_path <- paste0(tempdir(), "/reclass_rules_", rand_string, ".txt") + # fwrite(data, rules_path, sep = " ", col.names = TRUE) + + # The r.reclass function of GRASS GIS requires a text file + # including the old and the new value with an = between + # (e.g. 1 = 20) + rules <- data.table::data.table(old = data[[rast_val]], + equal = "=", + new = data[["reclass"]]) + } + # else { + # + # # The r.reclass function of GRASS GIS requires a text file + # # including the old and the new value with an = between + # # (e.g. 1 = 20) + # rules <- data.table::data.table(old = data[[rast_val]], + # equal = "=", + # new = data[[new_val]]) + # } # Create random string to attach to the file name of the temporary # rules .txt file rand_string <- stri_rand_strings(n = 1, length = 8, pattern = "[A-Za-z0-9]") @@ -193,14 +297,19 @@ reclass_raster <- function(data, rast_val, new_val, raster_layer, # Write rules as a .txt file to the temporary folder fwrite(rules, rules_path, sep = " ", col.names = FALSE) + # Remove all temporary data files + rm(indx_miss_raster, miss_raster, same_val1, indx_miss_rast_val, + miss_rast_val, missing_rast_values, same_val2, same_val) + if (sys_os == "linux" || sys_os == "osx") { - # Open GRASS GIS session - # Call external GRASS GIS command r.reclass - processx::run(system.file("sh", "reclass_raster.sh", - package = "hydrographr"), - args = c(raster_layer, rules_path, recl_layer, - no_data, type, compression_type, compression_level, bigtiff), - echo = !quiet) + + # Open GRASS GIS session + # Call external GRASS GIS command r.reclass + processx::run(system.file("sh", "reclass_raster.sh", + package = "hydrographr"), + args = c(raster_layer, rules_path, recl_layer, + no_data, type, compression_type, compression_level, bigtiff), + echo = !quiet) } else { # Check if WSL and Ubuntu are installed @@ -211,16 +320,16 @@ reclass_raster <- function(data, rast_val, new_val, raster_layer, wsl_rules_path <- fix_path(rules_path) wsl_sh_file <- fix_path( system.file("sh", "reclass_raster.sh", - package = "hydrographr")) + package = "hydrographr")) # Open GRASS GIS session on WSL # Call external GRASS GIS command r.reclass processx::run(system.file("bat", "reclass_raster.bat", - package = "hydrographr"), - args = c(wsl_raster_layer, wsl_rules_path, wsl_recl_layer, - no_data, type, compression_type, compression_level, bigtiff, - wsl_sh_file), - echo = !quiet) + package = "hydrographr"), + args = c(wsl_raster_layer, wsl_rules_path, wsl_recl_layer, + no_data, type, compression_type, compression_level, bigtiff, + wsl_sh_file), + echo = !quiet) } # Remove temporary rules file diff --git a/R/utils.R b/R/utils.R index 7013b72d..b7e822d7 100644 --- a/R/utils.R +++ b/R/utils.R @@ -11,7 +11,7 @@ get_os <- function() { sysinf <- Sys.info() if (!is.null(sysinf)) { - os <- sysinf["sysname"] + os <- sysinf[["sysname"]] if (os == "Darwin") { os <- "osx" } @@ -85,6 +85,6 @@ fix_path <- function(path) { stri_replace_all_fixed(., "\\", "/") %>% stri_replace_first_fixed(., drive, mnt) %>% stri_replace_first_fixed(., "Program Files (x86)", "PROGRA~2") %>% - stri_replace_first_fixed(., "Program Files", "PROGRA~1") + stri_replace_first_fixed(., "Program Files", "PROGRA~1") } diff --git a/hydrographr.Rproj b/hydrographr.Rproj index 488ba327..f0d61875 100644 --- a/hydrographr.Rproj +++ b/hydrographr.Rproj @@ -1,21 +1,21 @@ -Version: 1.0 - -RestoreWorkspace: Default -SaveWorkspace: Default -AlwaysSaveHistory: Default - -EnableCodeIndexing: Yes -UseSpacesForTab: Yes -NumSpacesForTab: 2 -Encoding: UTF-8 - -RnwWeave: Sweave -LaTeX: pdfLaTeX - -AutoAppendNewline: Yes -StripTrailingWhitespace: Yes - -BuildType: Package -PackageUseDevtools: Yes -PackageInstallArgs: --no-multiarch --with-keep.source -PackageRoxygenize: rd,collate,namespace,vignette +Version: 1.0 + +RestoreWorkspace: Default +SaveWorkspace: Default +AlwaysSaveHistory: Default + +EnableCodeIndexing: Yes +UseSpacesForTab: Yes +NumSpacesForTab: 2 +Encoding: UTF-8 + +RnwWeave: Sweave +LaTeX: pdfLaTeX + +AutoAppendNewline: Yes +StripTrailingWhitespace: Yes + +BuildType: Package +PackageUseDevtools: Yes +PackageInstallArgs: --no-multiarch --with-keep.source +PackageRoxygenize: rd,collate,namespace,vignette diff --git a/man/reclass_raster.Rd b/man/reclass_raster.Rd index 58d1c51e..35f48014 100644 --- a/man/reclass_raster.Rd +++ b/man/reclass_raster.Rd @@ -7,9 +7,12 @@ reclass_raster( data, rast_val, - new_val, + new_val = FALSE, raster_layer, recl_layer, + reclass_value = FALSE, + remaining = NULL, + remaining_value = -9999, no_data = -9999, type = "Int32", compression = "low", @@ -34,6 +37,16 @@ multiplying the values e.g. by 1000 to keep three decimals.} \item{recl_layer}{character. Full path of the output .tif layer, i.e., the reclassified raster file.} +\item{reclass_value}{integer. Value to reclassify the entire raster. +Default is FALSE. Note that in case reclass_value and new_val +is provided the raster is reclassified on reclass_value. +#' @param remaining character. How to treat raster values not listed in the reclassification table: + `"same"` retains their original values (equivalent to `* = *` in GRASS), + `"value"` assigns a fixed value (`remaining_value`), and + `NULL` (default) does nothing.} + +\item{remaining_value}{numeric. Value to assign if `remaining = "value"`. Default is -9999.} + \item{no_data}{numeric. The no_data value of the new .tif layer. Default is -9999.} @@ -74,10 +87,8 @@ my_directory <- tempdir() download_test_data(my_directory) # Read the stream order for each sub-catchment as a data.table -my_dt <- read_geopackage(gpkg= paste0(my_directory, - "/hydrography90m_test_data", - "/order_vect_59.gpkg"), - import_as = "data.table") +# my_dt <- read_geopackage(paste0(my_directory, "/order_vect_59.gpkg"), +type = "net", as_dt = T) # Select the stream segment ID and and the Strahler stream order @@ -97,10 +108,15 @@ str_ord_rast <- reclass_raster(data = str_ord, raster_layer = stream_raster, recl_layer = recl_raster) +str_ord_rast <- reclass_raster(data = str_ord$stream, + reclass_value = 1, + rast_val = "stream", + raster_layer = stream_raster, + recl_layer = recl_raster) } \references{ https://grass.osgeo.org/grass82/manuals/r.reclass.html } \author{ -Marlene Schürz +Marlene Schürz, Thomas Tomiczek, Afroditi Grigoropoulou } diff --git a/vignettes/example_other_stream_networks.Rmd b/vignettes/example_other_stream_networks.Rmd index 49e3e1ad..4b44dc3c 100644 --- a/vignettes/example_other_stream_networks.Rmd +++ b/vignettes/example_other_stream_networks.Rmd @@ -1,11 +1,20 @@ --- title: "Usage of hydrographr with other stream networks" +<<<<<<< HEAD +subtitle: "Examples showing how to use the hydrographr package with the NHDPlus dataset " +author: "" +date: "`r Sys.Date()`" +output: rmarkdown::html_vignette +vignette: > + %\VignetteIndexEntry{Other stream networks} +======= subtitle: "Examples showing how to use the hydrographr package with the NHDPlusHR dataset" author: "" date: "`r Sys.Date()`" output: html_vignette vignette: > %\VignetteIndexEntry{Usage of hydrographr with other stream networks} +>>>>>>> 51c1714db6631f1198315012d2bd5ab6cf66e1aa %\VignetteEncoding{UTF-8} %\VignetteEngine{knitr::rmarkdown} bibliography: bib/references_nhdp.bib @@ -41,8 +50,11 @@ library(leafem) library(htmlwidgets) library(mapview) library(here) +<<<<<<< HEAD +======= library(data.table) library(igraph) +>>>>>>> 51c1714db6631f1198315012d2bd5ab6cf66e1aa ``` Load required libraries: @@ -55,8 +67,11 @@ library(dplyr) library(archive) library(terra) library(leaflet) +<<<<<<< HEAD +======= library(data.table) library(igraph) +>>>>>>> 51c1714db6631f1198315012d2bd5ab6cf66e1aa ``` Define working directory: @@ -73,16 +88,27 @@ dir.create("data") ``` ```{r, eval = FALSE, include=FALSE} +<<<<<<< HEAD +wdir <- paste0("/home/mueblacker/Documents/Glowabio/Code/hydrographr/vignettes/data_hawaii") +``` + +### Download NHDPlus data +======= wdir <- paste0( "/home/mueblacker/Documents/Glowabio/Code/hydrographr/vignettes/data_hawaii") ``` ### Download NHDPlusHR data +>>>>>>> 51c1714db6631f1198315012d2bd5ab6cf66e1aa ```{r, eval=F} # Define the working directory for the region Hawaii wdir <- "my/working/directory/data_hawaii" +<<<<<<< HEAD +wdir <- "/home/afroditi/Documents/PhD/hydrographr_case_study/nhdplus" +======= +>>>>>>> 51c1714db6631f1198315012d2bd5ab6cf66e1aa data_dir <- paste0(wdir, "/data") rast_dir <- paste0(data_dir, "/raster") vect_dir <- paste0(data_dir, "/vector") @@ -97,15 +123,23 @@ dir.create(sp_dir) dir.create(out_dir) ``` +<<<<<<< HEAD +Data of the NHDPlus dataset can be downloaded using the R package +======= Data of the NHDPlusHR dataset can be downloaded using the R package +>>>>>>> 51c1714db6631f1198315012d2bd5ab6cf66e1aa `nhdplusTools` or from this USGS website: ```{r download, eval=F} # Download NHDPlus raster data for Hawaii download_nhdplusv2(outdir = rast_dir, +<<<<<<< HEAD + url = paste0("https://prd-tnm.s3.amazonaws.com/StagedProducts/", +======= url = paste0( "https://prd-tnm.s3.amazonaws.com/StagedProducts/", +>>>>>>> 51c1714db6631f1198315012d2bd5ab6cf66e1aa "Hydrography/NHDPlusHR/Beta/GDB/", "NHDPLUS_H_2006_HU4_RASTER.7z"), progress = TRUE) @@ -121,10 +155,16 @@ archive_extract(archive = paste0(rast_dir, "/NHDPLUS_H_2006_HU4_RASTER.7z"), # Download NHDPlus vector data for Hawaii download_nhdplusv2(outdir = vect_dir, +<<<<<<< HEAD + url = paste0("https://prd-tnm.s3.amazonaws.com/StagedProducts/", + "Hydrography/NHDPlusHR/Beta/GDB/", + "NHDPLUS_H_2006_HU4_GDB.zip"), +======= url = paste0( "https://prd-tnm.s3.amazonaws.com/StagedProducts/", "Hydrography/NHDPlusHR/Beta/GDB/", "NHDPLUS_H_2006_HU4_GDB.zip"), +>>>>>>> 51c1714db6631f1198315012d2bd5ab6cf66e1aa progress = TRUE) # Unzip -- can also be done manually @@ -192,12 +232,20 @@ m ![](../man/figures/hawaii_map.png) +<<<<<<< HEAD +By inspecting the NHDPlus data, we observe that the coordinate reference +system is not WGS84, but the EPSG:26904. + +```{r, eval=F} +rast(paste0(rast_dir, "/NHDPLUS_H_2006_HU4_RASTER/HRNHDPlusRasters2006/cat.tif")) +======= By inspecting the NHDPlusHR data, we observe that the coordinate reference system is not WGS84, but the EPSG:26904. ```{r, eval=F} rast(paste0(rast_dir, "/NHDPLUS_H_2006_HU4_RASTER/HRNHDPlusRasters2006/cat.tif")) +>>>>>>> 51c1714db6631f1198315012d2bd5ab6cf66e1aa ``` Therefore we will need to convert the species points' coordinates to @@ -287,6 +335,15 @@ spdata_ids ```{r, eval = FALSE} # Don't run: +<<<<<<< HEAD +# If we wanted to change the no data value of a layer, we could do it using the function: +set_no_data(data_dir = paste0(rast_dir, "/NHDPLUS_H_2006_HU4_RASTER/HRNHDPlusRasters2006"), var_layer = "cat.tif", + no_data = -32768) + +# Remember to set it back to the original value! + + +======= # If we wanted to change the no data value of a layer, we could do it using the # function: set_no_data(data_dir = paste0(rast_dir, "/NHDPLUS_H_2006_HU4_RASTER/HRNHDPlusRasters2006"), @@ -294,6 +351,7 @@ set_no_data(data_dir = paste0(rast_dir, "/NHDPLUS_H_2006_HU4_RASTER/HRNHDPlusRas no_data = -32768) # Remember to set it back to the original value! +>>>>>>> 51c1714db6631f1198315012d2bd5ab6cf66e1aa ``` ### Extract zonal statistics @@ -310,14 +368,22 @@ compatible with the coordinate system of CHELSA Bioclim files. ```{r, eval=F} # Download bio12 in the rast_dir +<<<<<<< HEAD +download.file("https://os.zhdk.cloud.switch.ch/envicloud/chelsa/chelsa_V2/GLOBAL/climatologies/1981-2010/bio/CHELSA_bio12_1981-2010_V.2.1.tif", +======= download.file( "https://os.zhdk.cloud.switch.ch/envicloud/chelsa/chelsa_V2/GLOBAL/climatologies/1981-2010/bio/CHELSA_bio12_1981-2010_V.2.1.tif", +>>>>>>> 51c1714db6631f1198315012d2bd5ab6cf66e1aa destfile = paste0(rast_dir, "/bio12_1981-2010.tif"), mode = "wb") # Crop bio12 to bounding box +<<<<<<< HEAD +# We first need to convert the bounding box coordinates to WGS84, that is the coordinate system of the CHELSA Bioclim layers. +======= # We first need to convert the bounding box coordinates to WGS84, that is the # coordinate system of the CHELSA Bioclim layers. +>>>>>>> 51c1714db6631f1198315012d2bd5ab6cf66e1aa bbox_wgs84 <- project(ext(bbox), from = "EPSG:26904", to = "EPSG:4326") # Now crop crop_to_extent(raster_layer = paste0(rast_dir, "/bio12_1981-2010.tif"), @@ -332,6 +398,17 @@ cat_rast <- rast(paste0(rast_dir, "/NHDPLUS_H_2006_HU4_RASTER/HRNHDPlusRasters20 # Reproject to WGS84 and write out the file with a new name cat_rast <- project(cat_rast, +<<<<<<< HEAD + y = "epsg:4326", + filename = paste0(rast_dir, "/NHDPLUS_H_2006_HU4_RASTER/HRNHDPlusRasters2006/cat_wgs84.tif")) + +# Extract the zonal statistics of bio12 based on the NHDPlus basins +zonal_stats <- extract_zonal_stat(data_dir = out_dir, + subc_id = "all", + subc_layer = paste0(rast_dir, "/NHDPLUS_H_2006_HU4_RASTER/HRNHDPlusRasters2006/cat_wgs84.tif"), + var_layer = "bio12_crop.tif", + quiet = F) +======= y = "epsg:4326", filename = paste0(rast_dir, "/NHDPLUS_H_2006_HU4_RASTER/HRNHDPlusRasters2006/cat_wgs84.tif")) @@ -341,6 +418,7 @@ zonal_stats <- extract_zonal_stat(data_dir = out_dir, subc_layer = paste0(rast_dir, "/NHDPLUS_H_2006_HU4_RASTER/HRNHDPlusRasters2006/cat_wgs84.tif"), var_layer = "bio12_crop.tif", quiet = F) +>>>>>>> 51c1714db6631f1198315012d2bd5ab6cf66e1aa # Same as before, we need to filter out the no data values to get a better # grasp of the resulting table @@ -395,11 +473,23 @@ gdb_dir <- paste0(vect_dir, "/NHDPLUS_H_2006_HU4_GDB") # Write out a .gpkg file based on the .gdb files get_nhdplushr(gdb_dir, file.path(gdb_dir, "nhdplus_2006-HU4.gpkg")) +<<<<<<< HEAD +======= +>>>>>>> 51c1714db6631f1198315012d2bd5ab6cf66e1aa # Import flow line attribute table as a data.table flow_line <- read_geopackage(gpkg = paste0(gdb_dir, "/nhdplus_2006-HU4.gpkg"), layer_name = "NHDFlowline") +<<<<<<< HEAD +# Import flow line as a graph +flow_line_graph <- read_geopackage(gpkg = paste0(gdb_dir, "/nhdplus_2006-HU4.gpkg"), + layer_name = "NHDFlowline", + import_as = "graph", name = "Permanent_Identifier") + +# Import the catchment layer as a spatial vector +catchment_sp <- read_geopackage(gpkg = paste0(gdb_dir, "/nhdplus_2006-HU4.gpkg"), +======= # To import flow line as a graph, first import the attribute table as a # data.table and then transform to a graph. Note that the first two columns @@ -413,6 +503,7 @@ flow_line_graph <- graph_from_data_frame(flow_line, directed = TRUE) # Import the catchment layer as a spatial vector catchment_sp <- read_geopackage( gpkg = paste0(gdb_dir, "/nhdplus_2006-HU4.gpkg"), +>>>>>>> 51c1714db6631f1198315012d2bd5ab6cf66e1aa layer_name = "NHDPlusCatchment", name = "REACHCODE", import_as = "SpatVect") @@ -421,13 +512,21 @@ catchment_sp <- read_geopackage( ### Snap points to the network +<<<<<<< HEAD +```{r} +======= ```{r, eval = F} +>>>>>>> 51c1714db6631f1198315012d2bd5ab6cf66e1aa # Define full path to the stream network raster layer stream_rast <- paste0(rast_dir, "/NHDPLUS_H_2006_HU4_RASTER/HRNHDPlusRasters2006/swnet.tif") # Define full path to the flow accumulation raster layer +<<<<<<< HEAD +flow_rast <- paste0(rast_dir, "/NHDPLUS_H_2006_HU4_RASTER/HRNHDPlusRasters2006/fac.tif") +======= flow_rast <- paste0(rast_dir, "/NHDPLUS_H_2006_HU4_RASTER/HRNHDPlusRasters2006/fac.tif") +>>>>>>> 51c1714db6631f1198315012d2bd5ab6cf66e1aa # Define thresholds for the flow accumulation of the stream segment, where # the point location should be snapped to