From 6cbf9a18893f09522a13604bfdb875d21db16e75 Mon Sep 17 00:00:00 2001 From: afrogri37 Date: Mon, 4 Sep 2023 17:10:07 +0200 Subject: [PATCH 1/6] update DESCRIPTION --- DESCRIPTION | 65 ++++++++++++++++++++++++++++++++++++++++++++--------- 1 file changed, 55 insertions(+), 10 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index df6ff057..931f13af 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -2,10 +2,10 @@ Package: hydrographr Type: Package Title: Scalable Hydrographic Data Processing in R Date: 2023 -Version: 0.1.9009 +Version: 1.0.15 Authors@R: - c(person("Maria", "Üblacker", - email = "maria.ueblacker@igb-berlin.de", + c(person("Marlene", "Schürz", + email = "marlene.schuerz@igb-berlin.de", role = c("aut", "cre")), person("Afroditi", "Grigoropoulou", email = "afroditi.grigoropoulou@igb-berlin.de", @@ -26,7 +26,7 @@ Authors@R: email = "thomas.tomiczek@igb-berlin.de", role = c("aut")), person("Vanessa", "Bremerich", - email = "bremerich@igb-berlin.de", + email = "vanessa.bremerich@igb-berlin.de", role = c("aut")), person("Giuseppe", "Amatulli", email = "giuseppe.amatulli@gmail.com", @@ -34,20 +34,40 @@ Authors@R: person("Sami", "Domisch", email = "sami.domisch@igb-berlin.de", role = c("aut"))) -Maintainer: Maria Üblacker -Description: A collection of R function wrappers for GDAL and GRASS-GIS functions to efficiently work with Hydrography90m spatial data. +Maintainer: Marlene Schürz and + Afroditi Grigoropoulou +Description: Scalable hydrographic geospatial data processing tools using + open-source command-line utilities. The package provides functions to + download the Hydrography90m data (https://essd.copernicus.org/articles/14/4525/2022/), + processing, reading and extracting information, as well as assessing network + distances and network connectivity. While the functions are, by default, + tailored towards the Hydrography90m data, they can also be generalized + towards other data and purposes, such as efficient cropping and merging + of raster and vector data, point-raster extraction, raster reclassification, + and data aggregation. The package depends on the open-source software + GDAL/OGR, GRASS-GIS and the AWK programming language in the Linux + environment, allowing a seamless language integration. Since the data is + processed outside R, hydrographr allows creating scalable geo-processing + workflows. Please see the installation guide of the additional software + at https://glowabio.github.io/hydrographr/articles/hydrographr.html. + Windows users need to to first activate the Windows Subsystem for Linux + (WSL) feature. License: GPL-3 URL: https://glowabio.github.io/hydrographr/ BugReports: https://github.com/glowabio/hydrographr/issues Encoding: UTF-8 LazyData: true -RoxygenNote: 7.2.2 +RoxygenNote: 7.2.3 Imports: processx (>= 3.7.0), data.table (>= 1.14.2), tidyr (>= 1.2.1), dplyr (>= 1.0.10), stringi (>= 1.7.8), + stringr (>= 1.4.1), + rlang (>= 1.0.6), + DBI (>= 1.1.3), + RSQLite (>= 2.2.19), terra (>= 1.6-41), sf (>= 1.0-9), parallel (>= 4.2.2), @@ -56,7 +76,32 @@ Imports: future (>= 1.29.0), doFuture (>= 0.12.2), future.apply (>= 1.10.0), - DBI (>= 1.1.3), - RSQLite (>= 2.2.19), memuse (>= 4.2-2), - igraph (>= 1.3.5) + igraph (>= 1.3.5), + magrittr (>= 2.0.3), + methods (>= 4.3.0) +Collate: + 'check_tiles_filesize.R' + 'crop_to_extent.R' + 'download_test_data.R' + 'download_tiles.R' + 'download_tiles_base.R' + 'extract_from_gpkg.R' + 'extract_ids.R' + 'extract_zonal_stat.R' + 'get_catchment_graph.R' + 'get_distance.R' + 'get_distance_graph.R' + 'get_pfafstetter_basins.R' + 'get_regional_unit_id.R' + 'get_segment_neighbours.R' + 'get_tile_id.R' + 'get_upstream_catchment.R' + 'utils.R' + 'merge_tile.R' + 'read_geopackage.R' + 'reclass_raster.R' + 'report_no_data.R' + 'set_no_data.R' + 'snap_to_network.R' + 'snap_to_subc_segment.R' From b3a0d5cade85681fc32392010ce6b37f5a5b6364 Mon Sep 17 00:00:00 2001 From: afrogri37 Date: Mon, 4 Sep 2023 17:17:12 +0200 Subject: [PATCH 2/6] added snapping --- vignettes/example_other_stream_networks.Rmd | 440 ++++++++++++++++++++ 1 file changed, 440 insertions(+) create mode 100644 vignettes/example_other_stream_networks.Rmd diff --git a/vignettes/example_other_stream_networks.Rmd b/vignettes/example_other_stream_networks.Rmd new file mode 100644 index 00000000..95a26256 --- /dev/null +++ b/vignettes/example_other_stream_networks.Rmd @@ -0,0 +1,440 @@ +--- +title: "Usage of hydrographr with other stream networks" +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} + %\VignetteEncoding{UTF-8} + %\VignetteEngine{knitr::rmarkdown} +bibliography: bib/references_nhdp.bib +link-citations: yes +linkcolor: blue +# csl: copernicus.csl +nocite: "@*" +editor_options: + markdown: + wrap: 72 +--- + +```{r, include = FALSE, eval = F} +# writes out the references for all packages +knitr::write_bib(file = "bib/references_nhdp.bib") +``` + +```{r setup, include=F} +knitr::opts_chunk$set(echo = TRUE, warning=FALSE, message=FALSE, + out.width='50%', fig.align='center') +knitr::opts_knit$set(root.dir = "./data_cuba") +library(hydrographr) +library(nhdplusTools) +library(rgbif) +library(terra) +library(archive) +library(tools) +library(dplyr) +library(knitr) +library(kableExtra) +library(leaflet) +library(leafem) +library(htmlwidgets) +library(mapview) +library(here) +``` + +Load required libraries: + +```{r, eval=T} +library(hydrographr) +library(nhdplusTools) +library(rgbif) +library(dplyr) +library(archive) +library(terra) +library(leaflet) +``` + +Define working directory: + +```{r, eval=FALSE} +wdir <- paste0(here(), "/vignettes/data_other_networks") +if(!dir.exists(paste0(wdir, "/data"))) dir.create(paste0(wdir, "/data")) +``` + +```{r, eval=FALSE, include=FALSE} +wdir <- "/home/afroditi/Documents/PhD/hydrographr_case_study/nhdplus" +setwd(wdir) +dir.create("data") +``` + +```{r, eval = FALSE, include=FALSE} +wdir <- paste0("/home/mueblacker/Documents/Glowabio/Code/hydrographr/vignettes/data_hawaii") +``` + +### Download NHDPlus data + +```{r, eval=F} + +# Define the working directory for the region Hawaii +wdir <- "my/working/directory/data_hawaii" +wdir <- "/home/afroditi/Documents/PhD/hydrographr_case_study/nhdplus" +data_dir <- paste0(wdir, "/data") +rast_dir <- paste0(data_dir, "/raster") +vect_dir <- paste0(data_dir, "/vector") +sp_dir <- paste0(data_dir, "/species") +out_dir <- paste0(wdir, "/output") + +# Create a new folder in the working directories to store all the data +dir.create(data_dir) +dir.create(rast_dir) +dir.create(vect_dir) +dir.create(sp_dir) +dir.create(out_dir) +``` + +Data of the NHDPlus dataset can be downloaded using the R package +`nhdplusTools` or from this USGS website: + + +```{r download, eval=F} +# Download NHDPlus raster data for Hawaii +download_nhdplusv2(outdir = rast_dir, + url = paste0("https://prd-tnm.s3.amazonaws.com/StagedProducts/", + "Hydrography/NHDPlusHR/Beta/GDB/", + "NHDPLUS_H_2006_HU4_RASTER.7z"), + progress = TRUE) + +# Unzip -- can also be done manually +# When called with default arguments extracts all files in the archive. +archive_extract(archive = paste0(rast_dir, "/NHDPLUS_H_2006_HU4_RASTER.7z"), + dir = paste0(rast_dir, "/NHDPLUS_H_2006_HU4_RASTER") +) +``` + +```{r, eval=F} + +# Download NHDPlus vector data for Hawaii +download_nhdplusv2(outdir = vect_dir, + url = paste0("https://prd-tnm.s3.amazonaws.com/StagedProducts/", + "Hydrography/NHDPlusHR/Beta/GDB/", + "NHDPLUS_H_2006_HU4_GDB.zip"), + progress = TRUE) + +# Unzip -- can also be done manually +archive_extract(archive = paste0(vect_dir, "/NHDPLUS_H_2006_HU4_GDB.zip"), + dir = paste0(vect_dir, "/NHDPLUS_H_2006_HU4_GDB") +) +``` + +Check directory structure + +```{r, eval = F} +list.files(paste0(rast_dir, "/NHDPLUS_H_2006_HU4_RASTER"), recursive = T, + pattern = ".tif$") +``` + +### Species data + +Download occurrence records of *Stenogobius hawaiiensis* from GBIF.org +using the R package `rgbif`. + +```{r, eval = FALSE} +# Download occurrence data with coordinates from GBIF +gbif_data <- occ_data(scientificName = "Stenogobius hawaiiensis", + hasCoordinate = TRUE) + +# To cite the data use: +# gbif_citation(gbif_data) + +# Clean the data +spdata <- gbif_data$data %>% + select(decimalLongitude, decimalLatitude, species, occurrenceStatus, + country, stateProvince, year) %>% + filter(!is.na(year), + stateProvince == "Hawaii") %>% + distinct() %>% + mutate(occurrence_id = 1:nrow(.)) %>% + rename("longitude" = "decimalLongitude", + "latitude" = "decimalLatitude") %>% + select(8, 1:7) + +# Save the data +write.csv(spdata, paste0(sp_dir, "/stenogobius_hawaiiensis.csv"), row.names = F, + quote = F) + +spdata +``` + +### Visualise species data + +```{r, include = FALSE, eval = FALSE} +m <- leaflet() %>% + addProviderTiles('Esri.WorldShadedRelief') %>% + addCircles(data = spdata, color = "purple") + +saveWidget(m, file=paste0(wdir, "../../man/figures/hawaii_map.html")) +``` + +```{r, eval = FALSE} +m <- leaflet() %>% + addProviderTiles('Esri.WorldShadedRelief') %>% + addCircles(data = spdata, color = "purple") + +m +``` + +![](../man/figures/hawaii_map.png) + +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")) +``` + +Therefore we will need to convert the species points' coordinates to +EPSG:26904, so that they match with the spatial layers. + +```{r, eval=F} + +# Convert the coordinate columns to a matrix +spdata_coord <- spdata %>% + select(longitude, latitude) %>% + as.matrix() + +# Project to the coordinate reference system of the NHDPlus dataset +spdata_coord <- project(spdata_coord, + from = "+proj=longlat +datum=WGS84", + to = "EPSG:26904") + +# Replace the WGS84 coordinates with the EPSG:26904 coordinates in the species +# points' data frame +spdata <- spdata %>% + mutate(longitude = spdata_coord[,1], + latitude = spdata_coord[,2]) + +``` + +### Crop and merge raster data + +```{r, eval = FALSE} +# The coordinates of the bounding box should also be in the NAD_1983_UTM_Zone_4N projection system +bbox <- c(582638, 2361564, 594373, 2370255) +crop_to_extent(raster_layer = paste0(rast_dir, "/NHDPLUS_H_2006_HU4_RASTER/HRNHDPlusRasters2006/cat.tif"), + bounding_box = bbox, + out_dir = out_dir, + file_name = "cat_crop.tif", + compression = "high" + ) + +# We can crop to another bounding box, and then try to merge the two cropped +# layers +bbox2 <- c(595603, 2371564, 604373, 2380255) +crop_to_extent(raster_layer = paste0(rast_dir, "/NHDPLUS_H_2006_HU4_RASTER/HRNHDPlusRasters2006/cat.tif"), + bounding_box = bbox2, + out_dir = out_dir, + file_name = "cat_crop2.tif", + compression = "high" + ) +# Merge the two cropped layers +merge_tiles(tile_dir = out_dir, + tile_names = c("cat_crop.tif", "cat_crop2.tif"), + out_dir = out_dir, + file_name = "cat_merged.tif", + compression = "high") +``` + +### Extract IDs, Report no data value, Set no data value + +```{r, eval=F} +spdata_ids <- extract_ids(data = spdata, + lon = "longitude", lat = "latitude", + id = "occurrence_id", + basin_layer = paste0(rast_dir, "/NHDPLUS_H_2006_HU4_RASTER/HRNHDPlusRasters2006/cat.tif")) + +spdata_ids + +# We observe that some points obtain basin_id=NA, because they fall out of the +# extent of the downloaded raster layer. We will filter these points out, +# keeping only those which overlap with the raster. + +spdata_ids <- spdata_ids %>% + filter(!is.na(basin_id)) + +spdata_ids + +# On the other hand, some other points obtain basin_id=-32768. This should be +# the no data value of the raster layer, meaning a value attributed to the sea +# cells. We can double-check this using the report_no_data function: +no_data_val <- report_no_data(data_dir = paste0(rast_dir, "/NHDPLUS_H_2006_HU4_RASTER/HRNHDPlusRasters2006"), var_layer = "cat.tif") + +no_data_val +# We can filter out these points, too, so that we only keep the points sampled +# in rivers. +spdata_ids <- spdata_ids %>% + filter(basin_id != no_data_val$NoData) +spdata_ids + +``` + +```{r, eval = FALSE} +# Don't run: +# 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! + + +``` + +### Extract zonal statistics + +To extract zonal statistics of a raster layer using the NHDPlus basins +as zones, we first need to download and crop to our extent an example +raster layer. We will use the CHELSA Bioclim 12 (BIO12) layer of +1981-2010 for this example. + +For this operation, we also need to convert the projection system of the +cat.tif raster to EPSG:4326 using the 'terra' package, so that it is +compatible with the coordinate system of CHELSA Bioclim files. + +```{r, eval=F} + +# Download bio12 in the rast_dir +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", +destfile = paste0(rast_dir, "/bio12_1981-2010.tif"), mode = "wb") + + +# Crop bio12 to bounding box +# We first need to convert the bounding box coordinates to WGS84, that is the coordinate system of the CHELSA Bioclim layers. +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"), + bounding_box = bbox_wgs84, + out_dir = out_dir, + file_name = "bio12_crop.tif", + compression = "high" + ) + +# Load the cat.tif layer +cat_rast <- rast(paste0(rast_dir, "/NHDPLUS_H_2006_HU4_RASTER/HRNHDPlusRasters2006/cat.tif")) + +# Reproject to WGS84 and write out the file with a new name +cat_rast <- project(cat_rast, + 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) + +# Same as before, we need to filter out the no data values to get a better +# grasp of the resulting table +zonal_stats <- zonal_stats %>% + filter(bio12_crop_data_cells !=0 ) +zonal_stats +``` + +### Reclassify raster + +```{r, eval=F} + +# Create example dataframe with reclassification rules to reclassify the value +# of all basins containing species data to 1, and all the rest to -99999. +reclass_rules <- data.frame(cat = as.integer(spdata_ids$basin_id), + new_val = as.integer(1)) + +# Apply reclassification +reclass_raster( + data = reclass_rules, + rast_val = "cat", + new_val = "new_val", + raster_layer = paste0(rast_dir, "/NHDPLUS_H_2006_HU4_RASTER/HRNHDPlusRasters2006/cat.tif"), + recl_layer = paste0(rast_dir, "/NHDPLUS_H_2006_HU4_RASTER/HRNHDPlusRasters2006/cat_reclass.tif"), + read = FALSE, + no_data = no_data_val$NoData) +``` + +### + + + +```{r, eval=F, include = F} +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")) + +extract_from_gpkg(data_dir = gdb_dir, subc_id = "all", var_layer = "nhdplus_2006-HU4.gpkg", quiet = F) + +# I think that stream IDs don’t correspond to catchment ids, so it’s not +# possible to filter the .gpkg using the ids of the catchment raster +``` + +### Read a .GPKG + +See [Download NHDPlus data] to download NHDPlus vector data for Hawaii. + +```{r, eval=F} +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")) + +# Import flow line attribute table as a data.table +flow_line <- read_geopackage(gpkg = paste0(gdb_dir, "/nhdplus_2006-HU4.gpkg"), + layer_name = "NHDFlowline") + +# 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"), + layer_name = "NHDPlusCatchment", + name = "REACHCODE", + import_as = "SpatVect") + +``` + +### Snap points to the network + +```{r} +# 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 +flow_rast <- paste0(rast_dir, "/NHDPLUS_H_2006_HU4_RASTER/HRNHDPlusRasters2006/fac.tif") + +# Define thresholds for the flow accumulation of the stream segment, where +# the point location should be snapped to +accu_threshold <- 700 +# Define the distance radius +dist_radius <- 20 + +# Snap point locations to the stream network +point_locations_snapped <- snap_to_network(data = spdata_ids, + lon = "longitude", + lat = "latitude", + id = "occurrence_id", + stream_layer = stream_rast, + accu_layer = flow_rast, + method = "both", + distance = dist_radius, + accumulation = accu_threshold, + quiet = FALSE) + +``` + +Note: The columns "subc_id_snap_dist" and "subc_id_snap_accu" of the +table "point_locations_snapped" are only meaningful in case that the +stream and flow accumulation raster layers were derived from the +Hydrography90m. + +### References From ab173b19f7b3824f2cfafc95dabc6a221debd159 Mon Sep 17 00:00:00 2001 From: tomitho <113701037+tomitho@users.noreply.github.com> Date: Fri, 6 Dec 2024 09:25:53 +0100 Subject: [PATCH 3/6] Update reclass_raster function Function now handles differences between raster values of the tif file and reclass input table by introducing NAs for values in the raster not found in the input table and not using values of the input table not found in the raster. Note: that all changes in the bat and sh file are due to an old dev_reclass raster branch and all updates are based on main branch. --- DESCRIPTION | 2 +- R/reclass_raster.R | 218 +++++++++++++++++++++++++++++------- inst/bat/reclass_raster.bat | 8 +- inst/sh/reclass_raster.sh | 18 +-- man/reclass_raster.Rd | 40 +++++-- 5 files changed, 219 insertions(+), 67 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index df6ff057..d8f384e2 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -41,7 +41,7 @@ URL: https://glowabio.github.io/hydrographr/ BugReports: https://github.com/glowabio/hydrographr/issues Encoding: UTF-8 LazyData: true -RoxygenNote: 7.2.2 +RoxygenNote: 7.3.2 Imports: processx (>= 3.7.0), data.table (>= 1.14.2), diff --git a/R/reclass_raster.R b/R/reclass_raster.R index 5ee2e2e8..286f998c 100644 --- a/R/reclass_raster.R +++ b/R/reclass_raster.R @@ -1,6 +1,6 @@ #' @title Reclassify an integer raster layer #' -#' @decription Reclassifies an integer raster .tif layer using the r.reclass +#' @description Reclassifies an integer raster .tif layer using the r.reclass #' function of GRASS GIS. To reclassify the raster layer the present raster #' values and the new raster values have to be defined. #' @@ -14,6 +14,9 @@ #' @param rast_val character. The name of the column with the present #' raster values. #' @param new_val character. The name of the column with the new raster values. +#' @param reclass_value integer. Value to reclassify the entire raster. +#' Default is FALSE. Note in case reclass_value is provided, new_val parameter +#' needs to be empty, otherwise the raster is reclassified based on the new_val column. #' @param raster_layer Full path to the input raster .tif layer. #' @param recl_layer character. Full path of the output .tif layer #' of the reclassified raster file. @@ -32,7 +35,7 @@ #' @importFrom terra rast #' @export #' -#' @author Maria M. Üblacker +#' @author Maria M. Üblacker, Thomas Tomiczek #' #' @references #' https://grass.osgeo.org/grass82/manuals/r.reclass.html @@ -45,10 +48,8 @@ #' 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_1264942.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 @@ -68,45 +69,62 @@ #' 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, read = TRUE, +reclass_raster <- function(data, rast_val, new_val = FALSE, raster_layer, + recl_layer, reclass_value = FALSE, no_data = -9999, type = "Int32", - compress = "DEFLATE", quiet = TRUE) { + compression = "low", bigtiff = TRUE, + read = FALSE, quiet = TRUE) { # Check operating system system <- get_os() # Check if data.frame is defined 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.")) - # 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.")) + + # 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)) @@ -136,9 +154,29 @@ reclass_raster <- function(data, rast_val, new_val, raster_layer, stop("type: Has to be 'Byte', 'Int16', 'UInt16', 'Int32', 'UInt32', 'CInt16', or 'CInt32' ") - # Check if compress is "DEFLATE" or "LZW" - if (!(compress == "DEFLATE" || compress == "LZW")) - stop("compress: Has to be 'DEFLATE or 'LZW'.") + # 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") { + compression_type <- "NONE" + compression_level <- 0 + } else if (compression == "low") { + compression_type <- "DEFLATE" + compression_level <- 2 + } else if (compression == "high") { + compression_type <- "DEFLATE" + compression_level <- 9 + } else { + stop("'compression' must be one of 'none', 'low', or 'high'.") + } + + # Define whether BIGTIFF is used or not. BIGTIFF is required for + # tiff output files > 4 GB. + if(bigtiff) { + bigtiff <- "YES" + } else { + bigtiff <- "NO" + } + # Check if quiet is logical if (!is.logical(quiet)) @@ -147,12 +185,98 @@ 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]") @@ -161,13 +285,17 @@ 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 (system == "linux" || system == "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, compress), + no_data, type, compression_type, compression_level, bigtiff), echo = !quiet) } else { @@ -179,27 +307,33 @@ 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, compress, 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 file.remove(rules_path) + if (file.exists(recl_layer)) { + # Print message + cat("Reclassified raster saved under: ", recl_layer, "\n") + } else { + stop("Output file was not written. File size may have been larger than 4GB", + "\nSet bigtiff = TRUE, for writing large output files.") + } + if (read == TRUE) { # Read reclassified .tif layer recl_rast <- rast(recl_layer) - - } else { - # Print message - print(paste0("Reclassified raster saved under: ", recl_layer)) + return(recl_rast) } } diff --git a/inst/bat/reclass_raster.bat b/inst/bat/reclass_raster.bat index 96f334c7..dceaf728 100644 --- a/inst/bat/reclass_raster.bat +++ b/inst/bat/reclass_raster.bat @@ -1,9 +1,9 @@ @echo on -set WSLENV=%1 %2 %3 %4 %5 %6 +set WSLENV=%1 %2 %3 %4 %5 %6 %7 %8 %9 set RULES=%2 -set SHDIR=%7 +set SHDIR=%9 wsl dos2unix %RULES% -wsl dos2unix %SHDIR% -wsl bash %SHDIR% $WSLENV \ No newline at end of file +wsl dos2unix %SHDIR% +wsl bash %SHDIR% $WSLENV diff --git a/inst/sh/reclass_raster.sh b/inst/sh/reclass_raster.sh index 2079c5f6..cfcf88bf 100755 --- a/inst/sh/reclass_raster.sh +++ b/inst/sh/reclass_raster.sh @@ -5,21 +5,23 @@ export RULES=$2 export OUTPUT=$3 export NODATA=$4 export TYPE=$5 -export COMP=$6 +export COMPRESSION=$6 +export LEVEL=$7 +export BTIFF=$8 + # Start GRASS GIS session -grass -f --text --tmp-location $RASTER <<'EOF' - - # Load raster input file +grass -f --gtext --tmp-location $RASTER <<'EOF' + + # Load raster input file r.in.gdal --o input=$RASTER output=raster --overwrite - - # Reclassify the raster according to the rules + + # Reclassify the raster according to the rules r.reclass input=raster output=recl_raster rules=$RULES --overwrite # Export reclassified raster map - r.out.gdal input=recl_raster output=$OUTPUT type=$TYPE format=GTiff nodata=$NODATA --o -f -m -c createopt="COMPRESS=$COMP" + r.out.gdal input=recl_raster output=$OUTPUT type=$TYPE format=GTiff nodata=$NODATA --o -f -m -c createopt="COMPRESS=$COMPRESSION,ZLEVEL=$LEVEL,BIGTIFF=$BTIFF" -# done EOF exit diff --git a/man/reclass_raster.Rd b/man/reclass_raster.Rd index 36f02d0c..8e3c2d0a 100644 --- a/man/reclass_raster.Rd +++ b/man/reclass_raster.Rd @@ -7,13 +7,15 @@ reclass_raster( data, rast_val, - new_val, + new_val = FALSE, raster_layer, recl_layer, - read = TRUE, + reclass_value = FALSE, no_data = -9999, type = "Int32", - compress = "DEFLATE", + compression = "low", + bigtiff = TRUE, + read = FALSE, quiet = TRUE ) } @@ -31,9 +33,9 @@ raster values.} \item{recl_layer}{character. Full path of the output .tif layer of the reclassified raster file.} -\item{read}{logical. If TRUE, then the reclassified raster .tif layer -gets read into R as a SpatRaster (terra object). -If FALSE, the layer is only stored on disk. Default is TRUE.} +\item{reclass_value}{integer. Value to reclassify the entire raster. +Default is FALSE. Note in case reclass_value is provided, new_val parameter +needs to be empty, otherwise the raster is reclassified based on the new_val column.} \item{no_data}{numeric. The no_data value of the new .tif layer. Default is -9999.} @@ -41,10 +43,21 @@ Default is -9999.} \item{type}{character. Data type; Options are Byte, Int16, UInt16, Int32, UInt32,CInt16, CInt32. Default is Int32.} +\item{read}{logical. If TRUE, then the reclassified raster .tif layer +gets read into R as a SpatRaster (terra object). +If FALSE, the layer is only stored on disk. Default is TRUE.} + \item{compress}{Compression type: DEFLATE or LZW. Default is DEFLATE.} } \description{ -Reclassify an integer raster layer +Reclassifies an integer raster .tif layer using the r.reclass +function of GRASS GIS. To reclassify the raster layer the present raster +values and the new raster values have to be defined. + +If the input raster layer has floating point values, you should multiply +the input data by some factor (e.g. 1000) to achieve integer values, +otherwise the GRASS GIS r.reclass will round the raster values down to +the next integer which is not always desired. } \examples{ # Download test data into the temporary R folder @@ -53,10 +66,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_1264942.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 @@ -76,10 +87,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{ -Maria M. Üblacker +Maria M. Üblacker, Thomas Tomiczek } From 3705220b138927219c59a360cf05effce1c876c3 Mon Sep 17 00:00:00 2001 From: afrogri37 Date: Thu, 17 Jul 2025 11:08:05 +0200 Subject: [PATCH 4/6] added arguments --- R/reclass_raster.R | 25 +++++++++++++++++++++---- 1 file changed, 21 insertions(+), 4 deletions(-) diff --git a/R/reclass_raster.R b/R/reclass_raster.R index 208739a9..f50dccbf 100644 --- a/R/reclass_raster.R +++ b/R/reclass_raster.R @@ -17,8 +17,13 @@ #' 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 +#' 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. +#' @param remaining_value numeric. Value to assign if `remaining = "value"`. Default is -9999. #' @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. @@ -46,7 +51,7 @@ #' @importFrom terra rast #' @export #' -#' @author Marlene Schürz, Thomas Tomiczek +#' @author Marlene Schürz, Thomas Tomiczek, Afroditi Grigoropoulou #' #' @references #' https://grass.osgeo.org/grass82/manuals/r.reclass.html @@ -89,6 +94,7 @@ reclass_raster <- function(data, rast_val, new_val = FALSE, raster_layer, recl_layer, reclass_value = FALSE, + remaining = NULL, remaining_value = -9999, no_data = -9999, type = "Int32", compression = "low", bigtiff = TRUE, read = FALSE, quiet = TRUE) { @@ -212,7 +218,7 @@ reclass_raster <- function(data, rast_val, new_val = FALSE, raster_layer, 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:", + 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),]) @@ -224,7 +230,7 @@ reclass_raster <- function(data, rast_val, new_val = FALSE, raster_layer, 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:", + 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),] @@ -256,6 +262,17 @@ reclass_raster <- function(data, rast_val, new_val = FALSE, raster_layer, rules <- data.table::data.table(old = data[[rast_val]], equal = "=", new = data[[new_val]]) + + if (!is.null(remaining)) { + if (remaining == "same") { + rules <- rbind(rules, data.table(old = "*", equal = "=", new = "*")) + } else if (remaining == "value") { + rules <- rbind(rules, data.table(old = "*", equal = "=", new = remaining_value)) + } else { + stop("remaining: must be one of 'same', 'value', or NULL.") + } + } + } } From 3bd73e8ba7b6bf0b423cf3df9c6c3d68de585067 Mon Sep 17 00:00:00 2001 From: afrogri37 Date: Thu, 17 Jul 2025 11:09:01 +0200 Subject: [PATCH 5/6] updates of branch --- DESCRIPTION | 38 +----------------------------- hydrographr.Rproj | 42 +++++++++++++++++----------------- man/download_future_climate.Rd | 7 +++++- man/get_pfafstetter_basins.Rd | 2 +- man/get_segment_neighbours.Rd | 2 +- man/get_upstream_variable.Rd | 2 +- man/reclass_raster.Rd | 18 +++++++++++---- 7 files changed, 44 insertions(+), 67 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 5fb38f8a..8c495ba1 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,13 +1,8 @@ Package: hydrographr Type: Package Title: Scalable Hydrographic Data Processing in R -<<<<<<< HEAD -Date: 2023 -Version: 1.0.15 -======= Date: 2024 Version: 1.2.0 ->>>>>>> 51c1714db6631f1198315012d2bd5ab6cf66e1aa Authors@R: c(person("Marlene", "Schürz", email = "marlene.schuerz@igb-berlin.de", @@ -32,12 +27,9 @@ Authors@R: role = c("aut")), person("Vanessa", "Bremerich", email = "vanessa.bremerich@igb-berlin.de", -<<<<<<< HEAD -======= role = c("aut")), person("Merret", "Buurman", email = "merret.buurman@igb-berlin.de", ->>>>>>> 51c1714db6631f1198315012d2bd5ab6cf66e1aa role = c("aut")), person("Giuseppe", "Amatulli", email = "giuseppe.amatulli@gmail.com", @@ -45,13 +37,9 @@ Authors@R: person("Sami", "Domisch", email = "sami.domisch@igb-berlin.de", role = c("aut"))) -<<<<<<< HEAD -Maintainer: Marlene Schürz and -======= Maintainer: Thomas Tomiczek , Merret Buurman , Marlene Schürz and ->>>>>>> 51c1714db6631f1198315012d2bd5ab6cf66e1aa Afroditi Grigoropoulou Description: Scalable hydrographic geospatial data processing tools using open-source command-line utilities. The package provides functions to @@ -68,23 +56,15 @@ Description: Scalable hydrographic geospatial data processing tools using workflows. Please see the installation guide of the additional software at https://glowabio.github.io/hydrographr/articles/hydrographr.html. Windows users need to to first activate the Windows Subsystem for Linux -<<<<<<< HEAD - (WSL) feature. -======= (WSL) feature. Instructions on how to use hydrographr with other, finer resolution stream networks, can be found at https://glowabio.github.io/hydrographr/articles/example_other_stream_networks.html ->>>>>>> 51c1714db6631f1198315012d2bd5ab6cf66e1aa License: GPL-3 URL: https://glowabio.github.io/hydrographr/ BugReports: https://github.com/glowabio/hydrographr/issues Encoding: UTF-8 LazyData: true -<<<<<<< HEAD -RoxygenNote: 7.2.3 -======= -RoxygenNote: 7.3.1 ->>>>>>> 51c1714db6631f1198315012d2bd5ab6cf66e1aa +RoxygenNote: 7.3.2 Imports: processx (>= 3.7.0), data.table (>= 1.14.2), @@ -112,21 +92,12 @@ Collate: 'crop_to_extent.R' 'download_test_data.R' 'download_tiles.R' -<<<<<<< HEAD -======= 'download_env.R' 'download_future_climate.R' ->>>>>>> 51c1714db6631f1198315012d2bd5ab6cf66e1aa 'download_tiles_base.R' 'extract_from_gpkg.R' 'extract_ids.R' 'extract_zonal_stat.R' -<<<<<<< HEAD - 'get_catchment_graph.R' - 'get_distance.R' - 'get_distance_graph.R' - 'get_pfafstetter_basins.R' -======= 'get_all_upstream_distances.R' 'get_catchment_graph.R' 'get_centrality.R' @@ -136,15 +107,11 @@ Collate: 'get_modelfit_table.R' 'get_pfafstetter_basins.R' 'get_predict_table.R' ->>>>>>> 51c1714db6631f1198315012d2bd5ab6cf66e1aa 'get_regional_unit_id.R' 'get_segment_neighbours.R' 'get_tile_id.R' 'get_upstream_catchment.R' -<<<<<<< HEAD -======= 'get_upstream_variable.R' ->>>>>>> 51c1714db6631f1198315012d2bd5ab6cf66e1aa 'utils.R' 'merge_tile.R' 'read_geopackage.R' @@ -153,7 +120,4 @@ Collate: 'set_no_data.R' 'snap_to_network.R' 'snap_to_subc_segment.R' -<<<<<<< HEAD -======= 'split_table.R' ->>>>>>> 51c1714db6631f1198315012d2bd5ab6cf66e1aa 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/download_future_climate.Rd b/man/download_future_climate.Rd index 0945ac5f..98b0662c 100644 --- a/man/download_future_climate.Rd +++ b/man/download_future_climate.Rd @@ -11,6 +11,7 @@ download_future_climate( scenario = NULL, model = NULL, tile_id = NULL, + version = "V.2.1", download_dir = ".", delete_zips = TRUE ) @@ -35,6 +36,8 @@ provided.} \item{tile_id}{character vector. The IDs of the requested tiles.} +\item{version}{string. e.g. "V.2.1". Defaults to "V.2.1".} + \item{download_dir}{character. The directory where the files will be downloaded. Default is the working directory.} @@ -98,6 +101,8 @@ TODO: The CHELSA units also have scale and offset! Explain them here?\tabular{ll | mpi-esm1-2-hr | TODO: Model description | | psl-cm6a-lr | TODO: Model description | | ukesm1-0-ll | TODO: MOdel description | + +TODO: Also add something on versions! } \note{ If there is an error during the download of a file @@ -122,7 +127,7 @@ variable = c("bio1", "bio2", "bio3", "bio4", "bio5", "bio6", "bio7", "bio8", "bio16", "bio17", "bio18", "bio19") scenario = c("ssp370", "ssp585") model = c("ipsl-cm6a-lr", "mpi-esm1-2-hr", "ukesm1-0-ll") - time_period = time_period, +time_period = c("2071-2100") download_future_climate(variable = variable, file_format = "zip", scenario = scenario, model = model, time_period = time_period, diff --git a/man/get_pfafstetter_basins.Rd b/man/get_pfafstetter_basins.Rd index 043d5d4e..f5ceb3b1 100644 --- a/man/get_pfafstetter_basins.Rd +++ b/man/get_pfafstetter_basins.Rd @@ -75,7 +75,7 @@ my_graph <- read_geopackage(gpkg = paste0(my_directory, # a random ID, i.e. it does not need to be the real outlet of the basin. g_subset <- get_catchment_graph(g = my_graph, subc_id = "513867227", - outlet = FALSE, + use_outlet = FALSE, mode = "in", as_graph = TRUE) diff --git a/man/get_segment_neighbours.Rd b/man/get_segment_neighbours.Rd index 63904184..b2d9eac6 100644 --- a/man/get_segment_neighbours.Rd +++ b/man/get_segment_neighbours.Rd @@ -97,7 +97,7 @@ my_graph <- read_geopackage(gpkg= paste0(my_directory, # Subset the graph and get a smaller catchment my_graph <- get_catchment_graph(g = my_graph, subc_id = 513866048, mode = "in", - outlet = FALSE, as_graph = TRUE, n_cores = 1) + use_outlet = FALSE, as_graph = TRUE, n_cores = 1) # Get a vector of all segment IDs diff --git a/man/get_upstream_variable.Rd b/man/get_upstream_variable.Rd index 47db5519..3d9f2605 100644 --- a/man/get_upstream_variable.Rd +++ b/man/get_upstream_variable.Rd @@ -109,7 +109,7 @@ my_graph <- read_geopackage(gpkg = paste0(my_directory, my_graph <- get_catchment_graph(g = my_graph, subc_id = 513867228, mode = "in", - outlet = FALSE, + use_outlet = FALSE, as_graph = TRUE, n_cores = 1) diff --git a/man/reclass_raster.Rd b/man/reclass_raster.Rd index 90787953..35f48014 100644 --- a/man/reclass_raster.Rd +++ b/man/reclass_raster.Rd @@ -11,6 +11,8 @@ reclass_raster( raster_layer, recl_layer, reclass_value = FALSE, + remaining = NULL, + remaining_value = -9999, no_data = -9999, type = "Int32", compression = "low", @@ -32,13 +34,19 @@ multiplying the values e.g. by 1000 to keep three decimals.} \item{raster_layer}{Full path to the input raster .tif layer.} -\item{reclass_value}{integer. Value to reclassify the entire raster. -Default is FALSE. Note in case reclass_value is provided, new_val parameter -needs to be empty, otherwise the raster is reclassified based on the new_val column.} - \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.} @@ -110,5 +118,5 @@ str_ord_rast <- reclass_raster(data = str_ord$stream, https://grass.osgeo.org/grass82/manuals/r.reclass.html } \author{ -Marlene Schürz, Thomas Tomiczek +Marlene Schürz, Thomas Tomiczek, Afroditi Grigoropoulou } From 2a37ad75e77b7eab36e4e23e1717a646d071cc4f Mon Sep 17 00:00:00 2001 From: afrogri37 Date: Thu, 17 Jul 2025 12:24:16 +0200 Subject: [PATCH 6/6] =?UTF-8?q?fixed=20os=20checking=C2=B4?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- R/reclass_raster.R | 99 +++++++++++++++++++--------------------------- R/utils.R | 4 +- 2 files changed, 43 insertions(+), 60 deletions(-) diff --git a/R/reclass_raster.R b/R/reclass_raster.R index f50dccbf..0dfb0df6 100644 --- a/R/reclass_raster.R +++ b/R/reclass_raster.R @@ -19,11 +19,6 @@ #' @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 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. -#' @param remaining_value numeric. Value to assign if `remaining = "value"`. Default is -9999. #' @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. @@ -51,7 +46,7 @@ #' @importFrom terra rast #' @export #' -#' @author Marlene Schürz, Thomas Tomiczek, Afroditi Grigoropoulou +#' @author Marlene Schürz, Thomas Tomiczek #' #' @references #' https://grass.osgeo.org/grass82/manuals/r.reclass.html @@ -64,8 +59,8 @@ #' download_test_data(my_directory) #' #' # Read the stream order for each sub-catchment as a data.table -#' # my_dt <- read_geopackage(paste0(my_directory, "/order_vect_59.gpkg"), -#' type = "net", as_dt = T) +#' my_dt <- read_geopackage(paste0(my_directory, "/hydrography90m_test_data/order_vect_59.gpkg"), +#' import_as = "data.table") #' #' #' # Select the stream segment ID and and the Strahler stream order @@ -94,7 +89,6 @@ reclass_raster <- function(data, rast_val, new_val = FALSE, raster_layer, recl_layer, reclass_value = FALSE, - remaining = NULL, remaining_value = -9999, no_data = -9999, type = "Int32", compression = "low", bigtiff = TRUE, read = FALSE, quiet = TRUE) { @@ -115,13 +109,13 @@ reclass_raster <- function(data, rast_val, new_val = FALSE, raster_layer, # Check if rast_val column names exist if (is.null(data[[rast_val]])) stop(paste0("rast_val: Column name '", rast_val, - "' does not exist.")) + "' does not exist.")) # 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.")) + " have to be integers.")) # Check if new_val column names exist when no reclass_values is given if (isFALSE(new_val) && isFALSE(reclass_value)) @@ -137,11 +131,11 @@ reclass_raster <- function(data, rast_val, new_val = FALSE, raster_layer, } # 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.")) - } + 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)) @@ -154,8 +148,8 @@ reclass_raster <- function(data, rast_val, new_val = FALSE, 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)) @@ -192,7 +186,7 @@ reclass_raster <- function(data, rast_val, new_val = FALSE, 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)) @@ -208,12 +202,12 @@ reclass_raster <- function(data, rast_val, new_val = FALSE, raster_layer, 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]])) { + # 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 @@ -253,32 +247,21 @@ reclass_raster <- function(data, rast_val, new_val = FALSE, raster_layer, # 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))) + 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 (!is.null(remaining)) { - if (remaining == "same") { - rules <- rbind(rules, data.table(old = "*", equal = "=", new = "*")) - } else if (remaining == "value") { - rules <- rbind(rules, data.table(old = "*", equal = "=", new = remaining_value)) - } else { - stop("remaining: must be one of 'same', 'value', or NULL.") - } - } - + # 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 + # use reclass_value for reclassification data$reclass <- reclass_value data$reclass <- as.integer(data$reclass) # @@ -296,7 +279,7 @@ reclass_raster <- function(data, rast_val, new_val = FALSE, raster_layer, 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 @@ -318,15 +301,15 @@ reclass_raster <- function(data, rast_val, new_val = FALSE, raster_layer, rm(indx_miss_raster, miss_raster, same_val1, indx_miss_rast_val, miss_rast_val, missing_rast_values, same_val2, same_val) - if (system == "linux" || system == "osx") { + 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 @@ -342,11 +325,11 @@ reclass_raster <- function(data, rast_val, new_val = FALSE, raster_layer, # 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") }