Experimenting with reading straight from the GRASS database

[florisvdh]

With read_RAST() and read_VECT() GRASS is currently instructed to write a temp file using a {v,r}.out.* module, which may cause extra overhead and take disk space. It may be shortcut by using the GDAL-GRASS GIS driver to read from the GRASS database. This is a standalone driver.

Current results from experimenting at the GRASS Community Meeting in Prague:

# If the drivers were compiled & installed from source 
# (https://github.com/OSGeo/gdal-grass), run:
# Sys.setenv(GDAL_DRIVER_PATH="/usr/lib/gdalplugins")
#
# But if installed from the ubuntugis-unstable PPA (`apt install libgdal-grass`),
# then it will just work out-of-the-box on Ubuntu
# See https://repology.org/project/gdal-grass/versions for other distros.

library(terra)
#> terra 1.7.29
library(rgrass)
#> Loading required package: XML
#> GRASS GIS interface loaded with GRASS version: (GRASS not running)

gdal(drivers = TRUE) |> 
  tibble::as_tibble() |> 
  dplyr::filter(stringr::str_detect(name, "GRASS|grass"))
#> # A tibble: 3 × 5
#>   name           type   can   vsi   long.name           
#>   <chr>          <chr>  <chr> <lgl> <chr>               
#> 1 GRASS          raster read  FALSE GRASS Rasters (7+)  
#> 2 GRASSASCIIGrid raster read  TRUE  GRASS ASCII Grid    
#> 3 OGR_GRASS      vector read  FALSE GRASS Vectors (5.7+)

initGRASS(
  home=tempdir(), 
  gisDbase="/home/floris/grassdata", 
  location="nc_basic_spm_grass7", 
  mapset="PERMANENT", 
  override=TRUE
  )
#> No gisBase set. Trying to detect from the GRASS_INSTALLATION environment variable.
#> No GRASS_INSTALLATION environment variable was found.
#> Trying to set gisBase by running command `grass --config path` (requires grass in the system PATH).
#> Taking gisBase value from `grass --config path` output: /usr/lib/grass82
#> gisdbase    /home/floris/grassdata 
#> location    nc_basic_spm_grass7 
#> mapset      PERMANENT 
#> rows        1350 
#> columns     1500 
#> north       228500 
#> south       215000 
#> west        630000 
#> east        645000 
#> nsres       10 
#> ewres       10 
#> projection:
#>  PROJCRS["NAD83(HARN) / North Carolina",
#>     BASEGEOGCRS["NAD83(HARN)",
#>         DATUM["NAD83 (High Accuracy Reference Network)",
#>             ELLIPSOID["GRS 1980",6378137,298.257222101,
#>                 LENGTHUNIT["metre",1]]],
#>         PRIMEM["Greenwich",0,
#>             ANGLEUNIT["degree",0.0174532925199433]],
#>         ID["EPSG",4152]],
#>     CONVERSION["SPCS83 North Carolina zone (meters)",
#>         METHOD["Lambert Conic Conformal (2SP)",
#>             ID["EPSG",9802]],
#>         PARAMETER["Latitude of false origin",33.75,
#>             ANGLEUNIT["degree",0.0174532925199433],
#>             ID["EPSG",8821]],
#>         PARAMETER["Longitude of false origin",-79,
#>             ANGLEUNIT["degree",0.0174532925199433],
#>             ID["EPSG",8822]],
#>         PARAMETER["Latitude of 1st standard parallel",36.1666666666667,
#>             ANGLEUNIT["degree",0.0174532925199433],
#>             ID["EPSG",8823]],
#>         PARAMETER["Latitude of 2nd standard parallel",34.3333333333333,
#>             ANGLEUNIT["degree",0.0174532925199433],
#>             ID["EPSG",8824]],
#>         PARAMETER["Easting at false origin",609601.22,
#>             LENGTHUNIT["metre",1],
#>             ID["EPSG",8826]],
#>         PARAMETER["Northing at false origin",0,
#>             LENGTHUNIT["metre",1],
#>             ID["EPSG",8827]]],
#>     CS[Cartesian,2],
#>         AXIS["easting (X)",east,
#>             ORDER[1],
#>             LENGTHUNIT["metre",1]],
#>         AXIS["northing (Y)",north,
#>             ORDER[2],
#>             LENGTHUNIT["metre",1]],
#>     USAGE[
#>         SCOPE["Engineering survey, topographic mapping."],
#>         AREA["United States (USA) - North Carolina - counties of Alamance; Alexander; Alleghany; Anson; Ashe; Avery; Beaufort; Bertie; Bladen; Brunswick; Buncombe; Burke; Cabarrus; Caldwell; Camden; Carteret; Caswell; Catawba; Chatham; Cherokee; Chowan; Clay; Cleveland; Columbus; Craven; Cumberland; Currituck; Dare; Davidson; Davie; Duplin; Durham; Edgecombe; Forsyth; Franklin; Gaston; Gates; Graham; Granville; Greene; Guilford; Halifax; Harnett; Haywood; Henderson; Hertford; Hoke; Hyde; Iredell; Jackson; Johnston; Jones; Lee; Lenoir; Lincoln; Macon; Madison; Martin; McDowell; Mecklenburg; Mitchell; Montgomery; Moore; Nash; New Hanover; Northampton; Onslow; Orange; Pamlico; Pasquotank; Pender; Perquimans; Person; Pitt; Polk; Randolph; Richmond; Robeson; Rockingham; Rowan; Rutherford; Sampson; Scotland; Stanly; Stokes; Surry; Swain; Transylvania; Tyrrell; Union; Vance; Wake; Warren; Washington; Watauga; Wayne; Wilkes; Wilson; Yadkin; Yancey."],
#>         BBOX[33.83,-84.33,36.59,-75.38]],
#>     ID["EPSG",3358]]

## RASTERS

# g.findfile hint comes from Vaclav Petras
filename_ras <- "elevation"
elevation_path_for_gdalgrass <- execGRASS(
  "g.findfile", 
  element = "cellhd", 
  file = filename_ras, 
  intern = TRUE,
  mapset = "PERMANENT"
  )[4]
elevation_path_for_gdalgrass <- regmatches(
  elevation_path_for_gdalgrass, 
  regexpr("(?<==['\"]).+(?=['\"]$)", elevation_path_for_gdalgrass, perl = TRUE)
  )
elevation_path_for_gdalgrass
#> [1] "/home/floris/grassdata/nc_basic_spm_grass7/PERMANENT/cellhd/elevation"

res <- rast(elevation_path_for_gdalgrass)
res
#> class       : SpatRaster 
#> dimensions  : 1350, 1500, 1  (nrow, ncol, nlyr)
#> resolution  : 10, 10  (x, y)
#> extent      : 630000, 645000, 215000, 228500  (xmin, xmax, ymin, ymax)
#> coord. ref. : +proj=lcc +lat_0=33.75 +lon_0=-79 +lat_1=36.1666666666667 +lat_2=34.3333333333333 +x_0=609601.22 +y_0=0 +ellps=GRS80 +units=m +no_defs 
#> source      : elevation 
#> color table : 1 
#> name        : elevation 
#> min value   :  55.57879 
#> max value   : 156.32986
inMemory(res)
#> [1] FALSE
# describe(elevation_path_for_gdalgrass)

stars::read_stars(elevation_path_for_gdalgrass)
#> stars object with 2 dimensions and 1 attribute
#> attribute(s), summary of first 1e+05 cells:
#>    elevation    
#>  121    : 2113  
#>  126    : 2086  
#>  133    : 1757  
#>  134    : 1645  
#>  135    : 1462  
#>  (Other): 7823  
#>  NA's   :83114  
#> dimension(s):
#>   from   to offset delta                       refsys x/y
#> x    1 1500 630000    10 +proj=lcc +lat_0=33.75 +l... [x]
#> y    1 1350 228500   -10 +proj=lcc +lat_0=33.75 +l... [y]
stars::read_stars(elevation_path_for_gdalgrass, proxy = TRUE)
#> stars_proxy object with 1 attribute in 1 file(s):
#> $elevation
#> [1] "[...]/elevation"
#> 
#> dimension(s):
#>   from   to offset delta                       refsys x/y
#> x    1 1500 630000    10 +proj=lcc +lat_0=33.75 +l... [x]
#> y    1 1350 228500   -10 +proj=lcc +lat_0=33.75 +l... [y]

## VECTORS

filename_vec <- "zipcodes"
zipcodes_path_for_gdalgrass <- execGRASS(
  "g.findfile", 
  element = "vector", 
  file = filename_vec, 
  intern = TRUE,
  mapset = "PERMANENT"
)[4]
zipcodes_path_for_gdalgrass <- regmatches(
  zipcodes_path_for_gdalgrass, 
  regexpr("(?<==['\"]).+(?=['\"]$)", zipcodes_path_for_gdalgrass, perl = TRUE)
) |> 
  paste0("/head")
zipcodes_path_for_gdalgrass
#> [1] "/home/floris/grassdata/nc_basic_spm_grass7/PERMANENT/vector/zipcodes/head"

vect(zipcodes_path_for_gdalgrass, layer = "zipcodes")
#> Warning in p@ptr$read(x, layer, query, extent, filter, proxy, what): GDAL Error
#> 1: Cannot reset cursor.
#> Warning in p@ptr$read(x, layer, query, extent, filter, proxy, what): GDAL Error
#> 1: Attributes not found.
#>  class       : SpatVector 
#>  geometry    : polygons 
#>  dimensions  : 48, 12  (geometries, attributes)
#>  extent      : 610047.9, 677060.7, 196327.5, 258102.6  (xmin, xmax, ymin, ymax)
#>  source      : head (zipcodes)
#>  coord. ref. : +proj=lcc +lat_0=33.75 +lon_0=-79 +lat_1=36.1666666666667 +lat_2=34.3333333333333 +x_0=609601.22 +y_0=0 +ellps=GRS80 +units=m +no_defs 
#>  names       :   cat OBJECTID WAKE_ZIPCO PERIMETER ZIPCODE_ ZIPCODE_ID
#>  type        : <int>    <int>      <num>     <num>    <num>      <num>
#>  values      :     1        1  2.625e+08 8.075e+04        2         27
#>                    2        2  3.163e+07 3.106e+04        7         25
#>                    3        3   5.44e+08 1.153e+05        8         20
#>      ZIPNAME    ZIPNUM           ZIPCODE        NAME SHAPE_Leng SHAPE_Area
#>        <chr>     <num>             <chr>       <chr>      <num>      <num>
#>    CREEDMOOR 2.752e+04   CREEDMOOR 27522   CREEDMOOR  8.078e+04  2.624e+08
#>  YOUNGSVILLE  2.76e+04 YOUNGSVILLE 27596 YOUNGSVILLE  3.106e+04  3.163e+07
#>      RALEIGH 2.762e+04     RALEIGH 27615     RALEIGH  1.153e+05   5.44e+08

sf::st_read(zipcodes_path_for_gdalgrass, layer = "zipcodes")
#> Reading layer `zipcodes' from data source 
#>   `/home/floris/grassdata/nc_basic_spm_grass7/PERMANENT/vector/zipcodes/head' 
#>   using driver `OGR_GRASS'
#> Simple feature collection with 48 features and 12 fields
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 610047.9 ymin: 196327.5 xmax: 677060.7 ymax: 258102.6
#> CRS:           NA

sf::sf_extSoftVersion()
#>           GEOS           GDAL         proj.4 GDAL_with_GEOS     USE_PROJ_H 
#>       "3.11.1"        "3.6.2"        "9.1.1"         "true"         "true" 
#>           PROJ 
#>        "9.1.1"
terra::gdal()
#> [1] "3.6.2"

^{Created on 2023-06-05 with reprex v2.0.2}

Session info

sessioninfo::session_info()
#> ─ Session info ───────────────────────────────────────────────────────────────
#>  setting  value
#>  version  R version 4.3.0 (2023-04-21)
#>  os       Linux Mint 21.1
#>  system   x86_64, linux-gnu
#>  ui       X11
#>  language nl_BE:nl
#>  collate  nl_BE.UTF-8
#>  ctype    nl_BE.UTF-8
#>  tz       Europe/Brussels
#>  date     2023-06-05
#>  pandoc   2.19.2 @ /usr/lib/rstudio/resources/app/bin/quarto/bin/tools/ (via rmarkdown)
#> 
#> ─ Packages ───────────────────────────────────────────────────────────────────
#>  package     * version   date (UTC) lib source
#>  abind         1.4-5     2016-07-21 [3] CRAN (R 4.2.0)
#>  class         7.3-22    2023-05-03 [3] RSPM (R 4.2.0)
#>  classInt      0.4-9     2023-02-28 [3] RSPM (R 4.2.0)
#>  cli           3.6.1     2023-03-23 [3] RSPM (R 4.2.0)
#>  codetools     0.2-19    2023-02-01 [3] RSPM (R 4.2.0)
#>  DBI           1.1.3     2022-06-18 [3] RSPM (R 4.2.0)
#>  digest        0.6.31    2022-12-11 [3] RSPM (R 4.2.0)
#>  dplyr         1.1.2     2023-04-20 [3] RSPM (R 4.2.0)
#>  e1071         1.7-13    2023-02-01 [3] RSPM (R 4.2.0)
#>  evaluate      0.21      2023-05-05 [3] RSPM (R 4.2.0)
#>  fansi         1.0.4     2023-01-22 [3] RSPM (R 4.2.0)
#>  fastmap       1.1.1     2023-02-24 [3] RSPM (R 4.2.0)
#>  fs            1.6.2     2023-04-25 [3] RSPM (R 4.2.0)
#>  generics      0.1.3     2022-07-05 [3] RSPM (R 4.2.0)
#>  glue          1.6.2     2022-02-24 [3] RSPM (R 4.2.0)
#>  htmltools     0.5.5     2023-03-23 [3] RSPM (R 4.2.0)
#>  KernSmooth    2.23-21   2023-05-03 [3] RSPM (R 4.2.0)
#>  knitr         1.43      2023-05-25 [3] RSPM (R 4.2.0)
#>  lifecycle     1.0.3     2022-10-07 [3] RSPM (R 4.2.0)
#>  lwgeom        0.2-13    2023-05-22 [3] RSPM (R 4.2.0)
#>  magrittr      2.0.3     2022-03-30 [3] RSPM (R 4.2.0)
#>  pillar        1.9.0     2023-03-22 [3] RSPM (R 4.2.0)
#>  pkgconfig     2.0.3     2019-09-22 [3] CRAN (R 4.0.1)
#>  proxy         0.4-27    2022-06-09 [3] RSPM (R 4.2.0)
#>  purrr         1.0.1     2023-01-10 [3] RSPM (R 4.2.0)
#>  R.cache       0.16.0    2022-07-21 [3] RSPM (R 4.2.0)
#>  R.methodsS3   1.8.2     2022-06-13 [3] RSPM (R 4.2.0)
#>  R.oo          1.25.0    2022-06-12 [3] RSPM (R 4.2.0)
#>  R.utils       2.12.2    2022-11-11 [3] RSPM (R 4.2.0)
#>  R6            2.5.1     2021-08-19 [3] RSPM (R 4.2.0)
#>  Rcpp          1.0.10    2023-01-22 [3] RSPM (R 4.2.0)
#>  reprex        2.0.2     2022-08-17 [3] RSPM (R 4.2.0)
#>  rgrass      * 0.3-9     2023-06-04 [1] local
#>  rlang         1.1.1     2023-04-28 [3] RSPM (R 4.2.0)
#>  rmarkdown     2.22      2023-06-01 [3] RSPM (R 4.2.0)
#>  rstudioapi    0.14      2022-08-22 [3] RSPM (R 4.2.0)
#>  sessioninfo   1.2.2     2021-12-06 [3] RSPM (R 4.2.0)
#>  sf            1.0-12    2023-03-19 [1] CRAN (R 4.3.0)
#>  stars         0.6-1     2023-04-06 [3] RSPM (R 4.2.0)
#>  stringi       1.7.12    2023-01-11 [3] RSPM (R 4.2.0)
#>  stringr       1.5.0     2022-12-02 [3] RSPM (R 4.2.0)
#>  styler        1.10.0    2023-05-24 [3] RSPM (R 4.2.0)
#>  terra       * 1.7-29    2023-04-22 [1] CRAN (R 4.2.3)
#>  tibble        3.2.1     2023-03-20 [3] RSPM (R 4.2.0)
#>  tidyselect    1.2.0     2022-10-10 [3] RSPM (R 4.2.0)
#>  units         0.8-2     2023-04-27 [3] RSPM (R 4.2.0)
#>  utf8          1.2.3     2023-01-31 [3] RSPM (R 4.2.0)
#>  vctrs         0.6.2     2023-04-19 [3] RSPM (R 4.2.0)
#>  withr         2.5.0     2022-03-03 [3] RSPM (R 4.2.0)
#>  xfun          0.39      2023-04-20 [3] RSPM (R 4.2.0)
#>  XML         * 3.99-0.14 2023-03-19 [3] RSPM (R 4.2.0)
#>  yaml          2.3.7     2023-01-23 [3] RSPM (R 4.2.0)
#> 
#>  [1] /home/floris/lib/R/library
#>  [2] /usr/local/lib/R/site-library
#>  [3] /usr/lib/R/site-library
#>  [4] /usr/lib/R/library
#> 
#> ──────────────────────────────────────────────────────────────────────────────

sf doesn't get a CRS though, and the others appear to get a PROJ string, maybe derived from the PROJ_INFO file (by the driver?) even though the EPSG code is set in the GRASS side.

@neteler @wenzeslaus

To be discussed.