Speed up and enhance reading from (and writing to) a GRASS GIS database

DESCRIPTION

@@ -12,11 +12,12 @@ Authors@R: c(
 Description: An interface between the 'GRASS' geographical information system ('GIS') and 'R', based on starting 'R' from within the 'GRASS' 'GIS' environment, or running a free-standing 'R' session in a temporary 'GRASS' location; the package provides facilities for using all 'GRASS' commands from the 'R' command line. The original interface package for 'GRASS 5' (2000-2010) is described in Bivand (2000) <doi:10.1016/S0098-3004(00)00057-1> and Bivand (2001) <https://www.r-project.org/conferences/DSC-2001/Proceedings/Bivand.pdf>. This was succeeded by 'spgrass6' for 'GRASS 6' (2006-2016) and 'rgrass7' for 'GRASS 7' (2015-present). The 'rgrass' package modernizes the interface for 'GRASS 8' while still permitting the use of 'GRASS 7'.
 Depends: R (>= 3.5.0)
 Imports: stats, utils, methods, xml2
-Suggests: terra (>= 1.6-16), sp (>= 0.9), knitr, rmarkdown, sf, stars, raster (>= 3.6-3)
+Suggests: terra (>= 1.6-16), sp (>= 0.9), knitr, rmarkdown, sf, stars, raster (>= 3.6-3), codetools
 VignetteBuilder: knitr
 SystemRequirements: GRASS (>= 7)
 License: GPL (>= 2)
 URL: https://rsbivand.github.io/rgrass/, https://grass.osgeo.org/, https://github.com/rsbivand/rgrass, https://lists.osgeo.org/mailman/listinfo/grass-stats
 BugReports: https://github.com/rsbivand/rgrass/issues/
-Collate: AAA.R options.R rgrass.R rast_link.R vect_link.R vect_link_ng.R initGRASS.R xml1.R
+Collate: AAA.R options.R rgrass.R rast_link.R vect_link.R vect_link_ng.R 
+  initGRASS.R xml1.R gdal_grass.R read_helpers.R
 

NEWS.md

@@ -1,6 +1,15 @@
 # **rgrass** version 0.4-3 (development)
 
-- see #87 - Windows QGIS standalone installations of GRASS can be used only if R is started in the OSGeo4W shell bundled with the installation
+- see #87 - Windows QGIS standalone installations of GRASS GIS can be used only if R is started in the OSGeo4W shell bundled with the installation
+
+- `write_VECT()`: when the `SpatVector` object already refers to a source file, an intermediate temporary file is no longer written to get the data into the GRASS GIS database (#93).
+A similar shortcut was already in place for `write_RAST()`.
+
+- `read_VECT()`: provide access to the standalone GDAL-GRASS driver to read vector data, which skips the step of writing a intermediate file (#93).
+Note that this standalone driver needs to be set up separately.
+More information is in the [driver's README](https://github.com/OSGeo/gdal-grass/blob/main/README.md).
+
+- `read_VECT()`: support reading as `SpatVectorProxy` class of `{terra}`, by providing a `proxy` argument (#93).
 
 # **rgrass** version 0.4-2 (2024-03-17)
 

R/gdal_grass.R

@@ -0,0 +1,43 @@
+gdal_has_grassraster_driver <- function() {
+  if (!(requireNamespace("terra", quietly = TRUE))) {
+    stop("terra is required to get the GDAL drivers list")
+  }
+  drv <- terra::gdal(drivers = TRUE)
+  "GRASS" %in% drv[drv$raster, ]$name
+}
+
+gdal_has_grassvector_driver <- function() {
+  if (!(requireNamespace("terra", quietly = TRUE))) {
+    stop("terra is required to get the GDAL drivers list")
+  }
+  drv <- terra::gdal(drivers = TRUE)
+  "OGR_GRASS" %in% drv[drv$vector, ]$name
+}
+
+generate_header_path <- function(name, type, ...) {
+  stopifnot(
+    is.character(type),
+    length(type) == 1L,
+    type %in% c("raster", "vector")
+  )
+  stopifnot(
+    is.character(name),
+    length(name) == 1L
+  )
+  element <- ifelse(type == "vector", "vector", "cellhd")
+  path <- execGRASS(
+    "g.findfile",
+    element = element,
+    file = name,
+    intern = TRUE,
+    ...
+  )[4]
+  path <- regmatches(
+    path,
+    regexpr("(?<==['\"]).+(?=['\"]$)", path, perl = TRUE)
+  )
+  if (type == "vector") {
+    path <- file.path(path, "head")
+  }
+  path
+}

R/rast_link.R

@@ -3,9 +3,10 @@
 #
 
 read_RAST <- function(
-    vname, cat = NULL, NODATA = NULL,
-    ignore.stderr = get.ignore.stderrOption(), return_format = "terra",
-    close_OK = return_format == "SGDF", flags = NULL) {
+    vname, cat = NULL, NODATA = NULL, return_format = "terra",
+    close_OK = return_format == "SGDF",
+    flags = NULL, Sys_ignore.stdout = FALSE,
+    ignore.stderr = get.ignore.stderrOption()) {
   if (!is.null(cat)) {
     if (length(vname) != length(cat)) {
       stop("vname and cat not same length")
@@ -17,7 +18,7 @@ read_RAST <- function(
   if (close_OK) {
     openedConns <- as.integer(row.names(showConnections()))
   }
-  stopifnot(is.logical(ignore.stderr))
+  stopifnot(is.logical(ignore.stderr), !is.na(ignore.stderr))
 
   if (!is.null(NODATA)) {
     if (any(!is.finite(NODATA)) || any(!is.numeric(NODATA))) {
@@ -32,10 +33,7 @@ read_RAST <- function(
     }
   }
 
-  msp <- unlist(strsplit(execGRASS("g.mapsets",
-    flags = "p",
-    intern = TRUE
-  ), " "))
+  msp <- get_mapsets()
 
   if (return_format == "SGDF") {
     if (!(requireNamespace("sp", quietly = TRUE))) {
@@ -77,37 +75,12 @@ read_RAST <- function(
       # 130422 at rgdal 0.8-8 GDAL.close(DS)
       # 061107 Dylan Beaudette NODATA
       # 071009 Markus Neteler's idea to use range
-      vca <- unlist(strsplit(vname[i], "@"))
-      if (length(vca) == 1L) {
-        exsts <- execGRASS("g.list",
-          type = "raster", pattern = vca[1],
-          intern = TRUE, ignore.stderr = ignore.stderr
-        )
-        if (length(exsts) > 1L) {
-          stop(
-            "multiple rasters named ", vca[1],
-            " found in in mapsets in search path: ",
-            paste(msp, collapse = ", "),
-            " ; use full path with @ to choose the required raster"
-          )
-        }
-        if (length(exsts) == 0L || exsts != vca[1]) {
-          stop(
-            vname[i], " not found in mapsets in search path: ",
-            paste(msp, collapse = ", ")
-          )
-        }
-      } else if (length(vca) == 2L) {
-        exsts <- execGRASS("g.list",
-          type = "raster", pattern = vca[1],
-          mapset = vca[2], intern = TRUE, ignore.stderr = ignore.stderr
-        )
-        if (length(exsts) == 0L || exsts != vca[1]) {
-          stop(vname[i], " not found in mapset: ", vca[2])
-        }
-      } else {
-        stop(vname[i], " incorrectly formatted")
-      }
+      vca <- sanitize_layername(
+        name = vname[i],
+        type = "raster",
+        mapsets = msp,
+        ignore.stderr = ignore.stderr
+      )
       typei <- NULL
       if (is.null(NODATA)) {
         tx <- execGRASS("r.info",
@@ -202,15 +175,17 @@ read_RAST <- function(
       if (!is.null(cat) && cat[i]) flags <- c(flags, "t")
       if (is.null(typei)) {
         execGRASS("r.out.gdal",
-          input = vname[i], output = tmplist[[i]],
-          format = drv, nodata = NODATAi, flags = flags,
-          ignore.stderr = ignore.stderr
+                  input = vname[i], output = tmplist[[i]],
+                  format = drv, nodata = NODATAi, flags = flags,
+                  ignore.stderr = ignore.stderr,
+                  Sys_ignore.stdout = Sys_ignore.stdout
         )
       } else {
         execGRASS("r.out.gdal",
-          input = vname[i], output = tmplist[[i]],
-          format = drv, nodata = NODATAi, type = typei, flags = flags,
-          ignore.stderr = ignore.stderr
+                  input = vname[i], output = tmplist[[i]],
+                  format = drv, nodata = NODATAi, type = typei, flags = flags,
+                  ignore.stderr = ignore.stderr,
+                  Sys_ignore.stdout = Sys_ignore.stdout
         )
       }
       reslist[[i]] <- getMethod("rast", "character")(tmplist[[i]])

R/read_helpers.R

@@ -0,0 +1,50 @@
+sanitize_layername <- function(name, type, mapsets, ignore.stderr) {
+  stopifnot(
+    is.character(type),
+    length(type) == 1L,
+    type %in% c("raster", "vector")
+  )
+  stopifnot(
+    is.character(name),
+    length(name) == 1L
+  )
+  vca <- unlist(strsplit(name, "@"))
+  if (length(vca) == 1L) {
+    exsts <- execGRASS("g.list",
+                       type = type, pattern = vca[1],
+                       intern = TRUE, ignore.stderr = ignore.stderr
+    )
+    if (length(exsts) > 1L) {
+      stop(
+        "multiple layers named ", vca[1],
+        " found in in mapsets in search path: ",
+        paste(mapsets, collapse = ", "),
+        " ; use full path with @ to choose the required raster"
+      )
+    }
+    if (length(exsts) == 0L || exsts != vca[1]) {
+      stop(
+        name, " not found in mapsets in search path: ",
+        paste(mapsets, collapse = ", ")
+      )
+    }
+  } else if (length(vca) == 2L) {
+    exsts <- execGRASS("g.list",
+                       type = type, pattern = vca[1],
+                       mapset = vca[2], intern = TRUE, ignore.stderr = ignore.stderr
+    )
+    if (length(exsts) == 0L || exsts != vca[1]) {
+      stop(name, " not found in mapset: ", vca[2])
+    }
+  } else {
+    stop(name, " incorrectly formatted")
+  }
+  vca
+}
+
+get_mapsets <- function() {
+  unlist(strsplit(
+    execGRASS("g.mapsets", flags = "p", intern = TRUE),
+    " "
+  ))
+}

R/vect_link_ng.R

@@ -2,20 +2,19 @@
 # Copyright (c) 2022 Roger S. Bivand
 #
 read_VECT <- function(
-    vname, layer, type = NULL, flags = "overwrite",
-    ignore.stderr = NULL) {
+    vname, layer = "", proxy = FALSE, use_gdal_grass_driver = TRUE, type = NULL,
+    flags = "overwrite", Sys_ignore.stdout = FALSE,
+    ignore.stderr = get.ignore.stderrOption()) {
   if (!(requireNamespace("terra", quietly = TRUE))) {
     stop("terra required for SpatVector output")
   }
-  if (is.null(ignore.stderr)) {
-    ignore.stderr <- get.ignore.stderrOption()
-  }
-  stopifnot(is.logical(ignore.stderr))
-  if (missing(layer)) layer <- "1"
-  layer <- as.character(layer)
+  stopifnot(is.logical(ignore.stderr), !is.na(ignore.stderr))
+  stopifnot(is.logical(use_gdal_grass_driver), !is.na(use_gdal_grass_driver))
   if (get.suppressEchoCmdInFuncOption()) {
     inEchoCmd <- set.echoCmdOption(FALSE)
   }
+  stopifnot(length(layer) == 1L)
+  if (!missing(layer)) layer <- as.character(layer)
   vinfo <- vInfo(vname)
   types <- names(vinfo)[which(vinfo > 0)]
   if (is.null(type)) {
@@ -24,13 +23,50 @@ read_VECT <- function(
     if (length(grep("areas", types)) > 0) type <- "area"
     if (is.null(type)) stop("Vector type not found")
   }
-  tf <- tempfile(fileext = ".gpkg")
-  execGRASS("v.out.ogr",
-    flags = flags, input = vname, type = type,
-    layer = as.character(layer), output = tf, output_layer = vname,
-    format = "GPKG", ignore.stderr = ignore.stderr
+  msp <- get_mapsets()
+  # in the v.out.ogr case we won't use vca, but this is done to run the checks
+  # on vname anyway:
+  vca <- sanitize_layername(
+    name = vname,
+    type = "vector",
+    mapsets = msp,
+    ignore.stderr = ignore.stderr
   )
-  res <- getMethod("vect", "character")(tf)
+  has_grassraster_drv <- gdal_has_grassraster_driver()
+
+  if (has_grassraster_drv && use_gdal_grass_driver) {
+    args <- list(name = vca[1], type = "vector")
+    if (length(vca) == 2L) args <- c(args, mapset = vca[2])
+    tf <- do.call(generate_header_path, args)
+    if (layer == "" && type == "area") {
+      layers <- terra::vector_layers(tf)
+      # Remove this condition once GDAL-GRASS driver issue
+      # has been solved (https://github.com/OSGeo/gdal-grass/issues/46).
+      # Then also move the type assignment code (from vInfo) to the
+      # v.out.ogr case, where it is used as an argument
+      layer <- layers[2]
+    }
+    # message(
+    #   "Will get data source ",
+    #   tf,
+    #   " (layername ",
+    #   ifelse(layer == "", "unknown, will get first layer", layer),
+    #   ")"
+    # )
+    res <- getMethod("vect", "character")(tf, layer, proxy = proxy)
+
+  } else {
+    if (layer == "") layer <- "1"
+    tf <- tempfile(fileext = ".gpkg")
+    execGRASS("v.out.ogr",
+              flags = flags, input = vname, type = type,
+              layer = layer, output = tf, output_layer = vname,
+              format = "GPKG", Sys_ignore.stdout = Sys_ignore.stdout,
+              ignore.stderr = ignore.stderr
+    )
+    # message("Reading ", tf)
+    res <- getMethod("vect", "character")(tf, proxy = proxy)
+  }
   if (!all(getMethod("is.valid", "SpatVector")(res))) {
     res <- getMethod("makeValid", "SpatVector")(res)
   }
@@ -40,27 +76,52 @@ read_VECT <- function(
   res
 }
 
-write_VECT <- function(x, vname, flags = "overwrite", ignore.stderr = NULL) {
+write_VECT <- function(x, vname, flags = "overwrite",
+                       ignore.stderr = get.ignore.stderrOption()) {
   if (!(requireNamespace("terra", quietly = TRUE))) {
     stop("terra required for SpatVector input")
   }
-  if (is.null(ignore.stderr)) {
-    ignore.stderr <- get.ignore.stderrOption()
-  }
-  stopifnot(is.logical(ignore.stderr))
+  stopifnot(is.logical(ignore.stderr), !is.na(ignore.stderr))
   if (get.suppressEchoCmdInFuncOption()) {
     inEchoCmd <- set.echoCmdOption(FALSE)
   }
+  srcs <- getMethod("sources", "SpatVector")(x)
+  if (length(srcs) == 1L) {
+    tf <- srcs
+  } else {
+    tf <- ""
+  }
+  # exit when the source is a GRASS database layer already:
+  if (grepl("[/\\\\]head::[^/\\\\]+$", tf)) {
+    grass_layername <- regmatches(
+      tf,
+      regexpr("(?<=[/\\\\]head::)[^/\\\\]+$", tf, perl = TRUE)
+    )
+    grass_dsn <- regmatches(
+      tf,
+      regexpr("(?<=[/\\\\])[^/\\\\]+(?=[/\\\\]head::)", tf, perl = TRUE)
+    )
+    stop(
+      "This SpatVector already links to layer '",
+      grass_layername,
+      "' of the data source '",
+      grass_dsn,
+      "' in the GRASS GIS database."
+    )
+  }
+
+  if (!file.exists(tf)) {
+    tf <- tempfile(fileext = ".gpkg")
+    getMethod("writeVector", c("SpatVector", "character"))(x, filename = tf,
+    filetype = "GPKG", overwrite = TRUE)
+  }
+
   type <- NULL
   if (getMethod("geomtype", "SpatVector")(x) == "points") type <- "point"
   if (getMethod("geomtype", "SpatVector")(x) == "lines") type <- "line"
   if (getMethod("geomtype", "SpatVector")(x) == "polygons") type <- "boundary"
   if (is.null(type)) stop("Unknown data class")
 
-  tf <- tempfile(fileext = ".gpkg")
-  getMethod("writeVector", c("SpatVector", "character"))(x, filename = tf,
-    filetype = "GPKG", overwrite = TRUE)
-
   execGRASS("v.in.ogr",
     flags = flags, input = tf, output = vname, type = type,
     ignore.stderr = ignore.stderr

R/xml1.R

@@ -536,8 +536,8 @@ execGRASS <- function(
       return(resOut)
     }
 
-    if (length(resOut) > 0) cat(resOut, sep = "\n")
-    if (length(resErr) > 0) cat(resErr, sep = "\n")
+    if (length(resOut) > 0 && !Sys_ignore.stdout) cat(resOut, sep = "\n")
+    if (length(resErr) > 0 && !ignore.stderr) cat(resErr, sep = "\n")
 
     attr(res, "resOut") <- resOut
     attr(res, "resErr") <- resErr