Nothing
#' Import data into 'GRASS.'
#'
#' This function loads a DEM (digital elevation model) and sites data (both required) into the 'GRASS' session.
#' Optionally, prediction sites and streams data can be loaded and the streams
#' may be corrected by snapping to prevent lose ends. Likewise, potential predictor
#' maps (raster or vector format) can be loaded.
#'
#' @param dem character; path to DEM (digital elevation model) raster file.
#' @param band integer (optional); defines which band of the dem to use
#' @param sites character string or object; path to sites vector file (ESRI shape)
#' or sp or sf data object.
#' @param streams character string or object (optional); path to network vector
#' file (ESRI shape) or sp or sf data object. If available this can be burnt into the DEM
#' in \code{\link{derive_streams}}
#' @param snap_streams boolean (optional); snap line ends.
#' If TRUE line ends of the streams are snapped to the next feature if they are
#' unconnected with threshold of 10 m using 'GRASS' function v.clean.
#' @param pred_sites character string vector or object(s) (optional); path(s) to prediction sites
#' vector files (ESRI shape) or sp or sf data object.
#' Different formats (i.e. path and objects) must not be mixed; more than one sf or sp
#' object must be provided as a list, not concatenated with \code{c}.
#' @param predictor_raster character vector (optional); paths to raster data to
#' import as predictors.
#' @param predictor_r_names character string vector (optional); names for potential predictor
#' variables in raster format; if not provided \code{perdictor_raster} is used.
#' @param predictor_vector character string vector of object(s) (optional); path(s)
#' to vector data (ESRI shape) or sp or sf object names to import as predictors.
#' Different formats (i.e. path and objects) must not be mixed; more than one sf or sp
#' object must be provided as a list, not concatenated with \code{c}.
#' @param predictor_v_names character vector (optional); names for potential predictor
#' variables in vector format ; if not provided \code{perdictor_vector} is used.
#'
#' @return Nothing, the data is loaded into the 'GRASS' session (mapset PERMANENT).
#' The DEM is stored as raster 'dem', sites as vector 'sites_o', prediction sites
#' as vector using the original file names with an appended '_o' (without extension),
#' streams as vector 'streams_o' in the 'GRASS' location. Additionally, predictor
#' raster map(s) can be read in and are stored in 'GRASS' using either the
#' original file names (without extension) or using the names provides in
#' \code{predictor_r_names}. The latter option may be useful if ArcGIS grid data
#' (typically stored as 'grid_name/w001001.adf') are used. Likewise, predictor
#' vector maps can be read in from Esri Shape file (given as the full file path)
#' or as sf or sp objects. Potential predictor data can also be read in later, e.g.
#' using GRASS commands \href{https://grass.osgeo.org/grass78/manuals/v.import.html}{v.import}
#' or \href{https://grass.osgeo.org/grass78/manuals/r.in.gdal.html}{r.in.gdal}
#' (see examples below).
#'
#' @details All vector data (sites, streams and potential predictors) is imported
#' into the current location using \href{https://grass.osgeo.org/grass78/manuals/v.import.html}{v.import}.
#' Hence, if the projections does not match to the one of the DEM (which was used
#' to specify the location in \code{\link{setup_grass_environment}}) the maps
#' are projected and imported on the fly.
#' All raster data are not transformed but it is assumed that they have the same
#' projection as the current location. Hence, it is important to make sure that
#' they all have indeed the same projection (and same cell size) and that the correct
#' one is set in \code{\link{setup_grass_environment}}. If this condition is not met,
#' the raster data should be preprocessed before importing.
#' Use \code{\link{check_projection}} to compare the projection of a raster data set and
#' the one of the current location (i.e the one of the dem).
#'
#' @note A GRASS session must be initiated and setup before, see \code{\link{setup_grass_environment}}.
#'
#' If sites, pred_sites and / or streams are sp objects it is important that they
#' have a datum defined otherwise the import will not work. Hence, it is e.g.
#' better to use proj4string = CRS("+proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=3500000 +y_0=0 +datum=potsdam +units=m +no_defs")
#' instead of proj4string = CRS("+proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=3500000 +y_0=0 +ellps=bessel +towgs84=598.1,73.7,418.2,0.202,0.045,-2.455,6.7 +units=m +no_defs"))
#' when defining sp objects. However, please note that proj4 is outdated
#' (\code{link{https://www.r-spatial.org/r/2020/03/17/wkt.html}}) and will at least rise warnings.
#'
#' @author Eduard Szoecs, \email{eduardszoecs@@gmial.com}, Mira Kattwinkel
#' \email{mira.kattwinkel@@gmx.net}
#' @export
#'
#' @examples
#' \donttest{
#' dem_path <- system.file("extdata", "nc", "elev_ned_30m.tif", package = "openSTARS")
#' if(.Platform$OS.type == "windows"){
#' grass_program_path = "c:/Program Files/GRASS GIS 7.6"
#' } else {
#' grass_program_path = "/usr/lib/grass78/"
#' }
#'
#' setup_grass_environment(dem = dem_path,
#' gisBase = grass_program_path,
#' remove_GISRC = TRUE,
#' override = TRUE
#' )
#' gmeta()
#'
#' # Load files into GRASS
#' dem_path <- system.file("extdata", "nc", "elev_ned_30m.tif", package = "openSTARS")
#' sites_path <- system.file("extdata", "nc", "sites_nc.shp", package = "openSTARS")
#' streams_path <- system.file("extdata", "nc", "streams.shp", package = "openSTARS")
#' preds_v_path <- system.file("extdata", "nc", "pointsources.shp", package = "openSTARS")
#' preds_r_path <- system.file("extdata", "nc", "landuse_r.tif", package = "openSTARS")
#'
#' import_data(dem = dem_path, sites = sites_path, streams = streams_path,
#' predictor_vector = preds_v_path, predictor_raster = preds_r_path)
#'
#' # Plot data
#' library(sp)
#' dem <- readRAST("dem", ignore.stderr = TRUE, plugin = FALSE)
#' sites_orig <- readVECT("sites_o", ignore.stderr = TRUE)
#' lu <- readRAST("landuse_r", ignore.stderr = TRUE, plugin = FALSE)
#' # import additional vector data
#' fp <- system.file("extdata", "nc", "pointsources.shp", package = "openSTARS")
#' execGRASS("v.import", flags = c("overwrite", "quiet"),
#' parameters = list(
#' input = fp,
#' output = "psources",
#' extent = "region"), # to import into current regien
#' intern = TRUE, ignore.stderr = TRUE)
#'
#' #plot(dem, col = terrain.colors(20))
#' plot(dem, col = grey.colors(20))
#' points(sites_orig, pch = 4)
#' ps <- readVECT("psources")
#' points(ps, bg = "red", pch = 21, col = "grey", cex = 1.5)
#'
#' # plot landuse data
#' library(raster)
#' par(mfcol = c(1,1), mar = c(5,4,4,2))
#' op <- par()
#' par(xpd = FALSE)
#' plot(raster(lu), legend = FALSE, xaxt = "n", yaxt = "n", bty = "n",
#' col = c("red", "goldenrod", "green", "forestgreen","darkgreen", "blue", "lightblue"))
#' par(xpd = TRUE)
#' legend("bottom", cex = 0.75,
#' legend = c("developed", "agriculture", "herbaceous",
#' "shrubland", "forest", "water", "sediment"),
#' fill = c("red", "goldenrod", "green",
#' "forestgreen","darkgreen", "blue", "lightblue"),
#' horiz = TRUE, inset = -0.175)
#' par <- op
#' }
import_data <- function(dem, band = 1, sites, streams = NULL, snap_streams = FALSE,
pred_sites = NULL, predictor_raster = NULL, predictor_r_names = NULL,
predictor_vector = NULL, predictor_v_names = NULL) {
if (nchar(get.GIS_LOCK()) == 0)
stop("GRASS not initialised. Please run setup_grass_environment().")
if (is.null(dem) | is.null(sites))
stop("DEM and sites are needed.")
# Import data -------------------
message("Loading DEM into GRASS as 'dem' ...")
# import raster
# MiKatt: it is necassary to set the region with g.region in setup_grass_environment;
# using flag "e" when importing the dem does not work
# (r.hydrodem produces very huge raster)
execGRASS("r.in.gdal",
flags = c("overwrite", "quiet", "o"), # overwrite projection check
parameters = list(
input = dem,
band = band,
output = "dem"),ignore.stderr=T)
message("Loading sites into GRASS as 'sites_o' ...")
# sites data
# flag "-r": only current region
# 20180216: not flag "o", because if projection is wrong I want to see an error
import_vector_data(data = sites, name = "sites_o", proj_ref_obj = dem)
# prediction sites data
if (!is.null(pred_sites)) {
if(is.character(pred_sites)){
pred_sites_names <- do.call(rbind,base::strsplit(sapply(pred_sites,basename,USE.NAMES=F), split="[.]"))[,1]
} else {
if(length(class(pred_sites)) == 1 && class(pred_sites) == "list"){
pred_sites_names <- paste("pred_sites", 1:length(pred_sites))
} else{
pred_sites_names <- "pred_sites"
}
}
if(any(pred_sites_names == "sites")){
message("Prediction sites cannot be named 'sites'. Please rename.")
} else{
pred_sites_names <- paste0(pred_sites_names, "_o")
message(paste0("Loading preditions sites into GRASS as ",
paste("'", pred_sites_names, "'", collapse=", ", sep=""), " ..."))
for(i in 1:length(pred_sites_names)){
import_vector_data(data = pred_sites[i], name = pred_sites_names[i], proj_ref_obj = dem)
}
}
}
# predictor raster maps
if (!is.null(predictor_raster)) {
if(is.null(predictor_r_names)){
predictor_r_names <- do.call(rbind,base::strsplit(sapply(predictor_raster,basename,USE.NAMES=F), split="[.]"))[,1]
}
#message(writeLines(strwrap(paste0("Loading raster predictor variables into GRASS as ",paste("'",predictor_r_names, "'", collapse = ", ", sep=""), " ..."),
# width = 80)))
message(paste0("Loading raster predictor variables into GRASS as ",paste("'",predictor_r_names, "'", collapse = ", ", sep=""), " ..."))
for(i in 1:length(predictor_r_names)){
execGRASS("r.in.gdal",
flags = c("overwrite","quiet","o"),
parameters = list(
input = predictor_raster[i],
output = predictor_r_names[i]),ignore.stderr=T)
}
}
# predictor vector maps
if (!is.null(predictor_vector)) {
if(is.null(predictor_v_names)){
if(is.character(predictor_vector)){
predictor_v_names <- do.call(rbind,base::strsplit(sapply(predictor_vector,basename,USE.NAMES=F), split="[.]"))[,1]
} else {
if(length(class(predictor_vector)) == 1 && class(predictor_vector) == "list"){
predictor_v_names <- paste("predictor_v", 1:length(predictor_vector))
} else{
predictor_v_names <- "predictor_v"
}
}
}
#message(writeLines(strwrap(paste0("Loading vector predictor variables into GRASS as ",paste("'", predictor_v_names, "'", collapse = ", ", sep=""), " ..."),
# width = 80), con = stderr()), appendLF = TRUE)
message(paste0("Loading vector predictor variables into GRASS as ",paste("'", predictor_v_names, "'", collapse = ", ", sep=""), " ..."))
for(i in 1:length(predictor_v_names)){
import_vector_data(data = predictor_vector[i], name = predictor_v_names[i], proj_ref_obj = dem)
}
}
# streams data
if (!is.null(streams)) {
message("Loading streams into GRASS as 'streams_o' ...")
# flag "-r": only current region
import_vector_data(data = streams, name = "streams_o", proj_ref_obj = dem)
# MiKatt: snapp line ends to next vertext to prevent loose ends/ unconnected streams and to build topography
if(snap_streams){
execGRASS("v.clean",
flags=c("c","overwrite","quiet"),
parameters = list(
input = "streams_o",
type = "line",
output = "streams_oc",
tool = "snap",
threshold = 10),ignore.stderr=T)
#sink("temp.txt")
execGRASS("g.copy",
flags = c("overwrite", "quiet"),
parameters = list(
vector = "streams_oc,streams_o"), intern = TRUE, ignore.stderr = TRUE)
#sink()
execGRASS("g.remove",
flags = c("quiet", "f"),
parameters = list(
type = "vector",
name = "streams_oc"
))
}
} else {
message("No streams available, skipping.\n")
}
try(unlink("temp.txt"), silent = TRUE)
}
#' Import vector data into GRASS.
#'
#' Generic function to import vector data of various formats into the GRASS environment.
#'
#' @param data character string or object; path to data vector file (shape), postgis
#' data source name (dsn; see details), or sp or sf data object.
#' @param name string giving the base name of the vector data within the GRASS environment (i.e. output)
#' @param layer character string; default 1, particularly needed if data is a dsn for
#' importing postgis data (see details)
#' @param proj_ref_obj character; path to a georeferenced data file to be used as reference; only
#' used if \code{data} is an sf of sp object, then, the two projections are compared to check if
#' on-the-fly reprojection is needed for importing; typically the dem raster file used in this project.
#' @param snap float; snapping threshold in map units. If != -1 (default) vertices are snapped to other vertices
#' in this snapping distance during import. If used, the features are automatically cleaned afterwards
#' (see GRASS tools \href{https://grass.osgeo.org/grass78/manuals/v.import.html}{v.import}
#' and \href{https://grass.osgeo.org/grass78/manuals/v.clean.html}{v.clean} )
#'
#' @return Nothing.
#'
#' @details For importing data from Postgis, all data base credentials must be supplied
#' in \code{data} and the correct \code{layer} and, if the table containing the polygons
#' are in a specific schema also that one (see example)
#'
#'@examples
#'# import data from Postgis
#'\dontrun{
#' import_vector_data(data = "PG: 'pgname=postgit_DB', 'host=123.45.67.890',
#' 'port='1234', 'user=username', 'password=password'",
#' name = "forest", layer = "landuse_schema.forest")
#'}
#'
#' @author Mira Kattwinkel, \email{mira.kattwinkel@@gmx.net}
import_vector_data <- function(data, name, layer = NULL, proj_ref_obj = NULL, snap = -1){
# flag "-r": only current region
# import_flag <- TRUE
if(inherits(data, 'sf')){
data <- as(data, 'Spatial')
}
if(inherits(data, 'Spatial')) {
# if(!is.null(proj_ref_obj)){
# MIKatt 20200720
# writeVECT also " moves a Spatial*DataFrame object through a temporary shapefile to a GRASS vector object file".
# Hence, just write temp shapefile and use v.import without projection check here!
# import_flag not needed
# MiKatt 20200717 wkt option but not much better than proj4
# temp_dir <- tempdir()
# wkt_file <- file.path(temp_dir, "wkt_data.prj")
# write(paste(wkt(data)), file = wkt_file)
# proj_data <- execGRASS("g.proj", flags = c("w", "f"),
# parameters = list(
# wkt = wkt_file
# ), intern = TRUE)
# proj_ref <- execGRASS("g.proj", flags = c("w", "f"),
# parameters = list(
# georef = proj_ref_obj
# ), intern = TRUE, ignore.stderr = TRUE)
# invisible(file.remove(file.path(temp_dir,"wkt_data.prj")))
# proj_data <- execGRASS("g.proj", flags = c("j", "f"),
# parameters = list(
# proj4 = proj4string(data)
# ), intern = TRUE)[1] # sites_in@proj4string@projargs does not work!
# proj_ref <- execGRASS("g.proj", flags = c("j", "f"),
# parameters = list(
# georef = proj_ref_obj
# ), intern = TRUE, ignore.stderr = TRUE)
# if(identical(proj_ref, proj_data)) {
# writeVECT(data, name, v.in.ogr_flags = c("overwrite", "quiet", "r", "o"), # "o": overwrite projection check
# ignore.stderr=TRUE)
# import_flag <- FALSE
# }
# }
#if(import_flag) {
rgdal::writeOGR(obj = data, dsn = tempdir(), layer = name, driver="ESRI Shapefile", overwrite_layer = TRUE)
data <- file.path(tempdir(), paste0(name, ".shp"))
#}
}
#if(import_flag) {
# gives error on Linux if projection is not identical but it works!
# 'intern' and 'ignore.stderr' to suppress error messages
if(is.null(layer)){
execGRASS("v.import", flags = c("overwrite", "quiet"),
parameters = list(
input = data,
output = name,
snap = snap,
extent = "region"), # to import into current region (= flags("r") in v.in.ogr)
intern = TRUE, ignore.stderr = TRUE)
} else {
execGRASS("v.import", flags = c("overwrite", "quiet"),
parameters = list(
input = data,
layer = layer,
output = name,
snap = snap,
extent = "region"), # to import into current region (= flags("r") in v.in.ogr)
intern = TRUE, ignore.stderr = TRUE)
#}
if(snap != -1){
execGRASS("v.clean", flags = c("c", "quiet"),
parameters = list(
input = name,
output = paste0(name, "_c"),
tool = "break,rmdupl,rmsa"
))
execGRASS("v.copy", flags = c("overwrite", "quiet"),
parameters = list(
vector(paste0(name, "_c"), name)
))
execGRASS("g.remove", flags = c("f", "quiet"),
parameters = list(
type = "vector",
name = paste0(name, "_c")
))
}
if(file.exists(file.path(tempdir(), paste0(name, ".shp")))){
invisible(file.remove(file.path(tempdir(), list.files(path = tempdir(), pattern = paste0(name, ".")))))
}
}
}
#' Show the projection of raster data and compare to the current GRASS location.
#'
#' @description Check if the projection of raster files matches the region of the
#' current location
#'
#' @param path character string vector; path raster data file(s)
#' @param format character string; how to format the output (see details)
#'
#' @return Nothing.
#'
#' @details This is a wrapper for the GRASS function
#' \href{https://grass.osgeo.org/grass78/manuals/g.proj.html}{g.proj}.
#' It prints out the projection information of the current location
#' and of the raster file(s). Based on this information it can be decided if the data
#' can be read into GRASS (\code{\link{import_data}}) without prior processing, i.e.
#' if all raster data are of the same projection.
#' Different output options can be chosen:
#' * \code{wkt} WKT format (default)
#' * \code{grass} conventional GRASS format
#' * \code{shell} shell script style
#' * \code{proj4} PROJ.4 format (note that this format is deprecated)
#'
#' @export
check_projection <- function(path, format = "wkt"){
format_flag <- "w"
if(format == "grass"){
format_flag <- "p"
} else if(format == "shell"){
format_flag <- "g"
} else if(format == "proj4"){
format_flag <- "j"
}
loc_proj <- execGRASS("g.proj", flags = c(format_flag),
intern = TRUE)
message("current location")
print(loc_proj)
n <- length(path)
for(i in 1:n){
rast_proj <- execGRASS("g.proj", flags = c(format_flag),
parameters = list(
georef = path[i]
), intern = TRUE)
message(paste0(basename(path[i])))
if(identical(loc_proj, rast_proj)){
message("Projection seems to match current location")
} else {
message("Projection seem to be different from current location; please consider reprojection before importing.")
}
print(rast_proj)
}
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.