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"
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.
With
read_RAST()andread_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:
Created on 2023-06-05 with reprex v2.0.2
Session info
``` r 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.