Coercion between object formats
Roger Bivand
Introduction
The original R-GRASS interface (R. S. Bivand 2000; Neteler and Mitasova 2008) was designed to move raster and later vector data between R and GRASS GIS. To do this, use was made of intermediate files, often using the external GDAL library on both sides. On the R side, the rgdal package was used, interfacing GDAL and PROJ as GRASS GIS also did. The GRASS commands r.in.gdal, r.out.gdal, v.in.ogr and v.out.ogr were matched by rgdal functions using the same underlying external libraries:
GDAL was supplemented for raster data by simply reading and writing uncompressed binary files using r.in.bin and r.out.bin, with custom functions on the R side. As then written, the R-GRASS interface used sp classes for both raster and vector data, supplemented more recently with sf classes for vector data only.
The current version of the R-GRASS interface has been simplified to use the terra package because it, like sf and rgdal before it, links to the important external libraries. The workhorse driver is known as RRASTER, and has been widely used in raster and terra (see also (https://rspatial.org)). It uses GDAL but writes a flat uncompressed binary file. Using terra::rast() also appears to preserve category names and colour tables, but needs further testing (see (https://github.com/rsbivand/rgrass/issues/42)).
From GDAL 3.5.0, the RRASTER driver also supports WKT2_2019 CRS representations; in earlier versions of GDAL, the driver only supported the proj-string representation (https://github.com/rsbivand/rgrass/issues/51).
These changes mean that users transferring data between R and GRASS will need to coerce between terra classes SpatVector and SpatRaster and the class system of choice. In addition, SpatRaster is only read into memory from file when this is required, so requiring some care.
On loading and attaching, terra displays its version:
terra_available <- require("terra", quietly=TRUE)## terra 1.6.7sf_available <- require("sf", quietly=TRUE)## Linking to GEOS 3.11.0, GDAL 3.5.1, PROJ 9.0.1; sf_use_s2() is TRUESys.setenv("_SP_EVOLUTION_STATUS_"="2")
sp_available <- require("sp", quietly=TRUE)
stars_available <- require("stars", quietly=TRUE)
raster_available <- require("raster", quietly=TRUE)This description is constructed conditioning on the availability of suggested packages. terra::gdal() tells us the versions of the external libraries being used by terra:
gdal(lib="all")## gdal proj geos
## "3.5.1" "9.0.1" "3.11.0"When using CRAN binary packages built static for Windows and macOS, the R packages will use the same versions of the external libraries, but not necessarily the same versions as those against which GRASS was installed.
"SpatVector" coercion
In the terra package (Hijmans 2022b), vector data are held in "SpatVector" objects. This means that when read_VECT() is used, a "SpatVector" object is returned, and the same class of object is needed for write_VECT() for writing to GRASS.
fv <- system.file("ex/lux.shp", package="terra")
(v <- vect(fv))## class : SpatVector
## geometry : polygons
## dimensions : 12, 6 (geometries, attributes)
## extent : 5.74414, 6.528252, 49.44781, 50.18162 (xmin, xmax, ymin, ymax)
## source : lux.shp
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## names : ID_1 NAME_1 ID_2 NAME_2 AREA POP
## type : <num> <chr> <num> <chr> <num> <int>
## values : 1 Diekirch 1 Clervaux 312 18081
## 1 Diekirch 2 Diekirch 218 32543
## 1 Diekirch 3 Redange 259 18664These objects are always held in memory, so there is no inMemory() method:
try(inMemory(v))## Error in (function (classes, fdef, mtable) :
## unable to find an inherited method for function 'inMemory' for signature '"SpatVector"'The coordinate reference system is expressed in WKT2-2019 form:
cat(crs(v), "\n")## GEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["geodetic latitude (Lat)",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["geodetic longitude (Lon)",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4326]]"sf"
Most new work should use vector classes defined in the sf package (Pebesma 2022, 2018). In this case, coercion uses st_as_sf():
v_sf <- st_as_sf(v)
v_sf## Simple feature collection with 12 features and 6 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 5.74414 ymin: 49.44781 xmax: 6.528252 ymax: 50.18162
## Geodetic CRS: WGS 84
## First 10 features:
## ID_1 NAME_1 ID_2 NAME_2 AREA POP
## 1 1 Diekirch 1 Clervaux 312 18081
## 2 1 Diekirch 2 Diekirch 218 32543
## 3 1 Diekirch 3 Redange 259 18664
## 4 1 Diekirch 4 Vianden 76 5163
## 5 1 Diekirch 5 Wiltz 263 16735
## 6 2 Grevenmacher 6 Echternach 188 18899
## 7 2 Grevenmacher 7 Remich 129 22366
## 8 2 Grevenmacher 12 Grevenmacher 210 29828
## 9 3 Luxembourg 8 Capellen 185 48187
## 10 3 Luxembourg 9 Esch-sur-Alzette 251 176820
## geometry
## 1 POLYGON ((6.026519 50.17767...
## 2 POLYGON ((6.178368 49.87682...
## 3 POLYGON ((5.881378 49.87015...
## 4 POLYGON ((6.131309 49.97256...
## 5 POLYGON ((5.977929 50.02602...
## 6 POLYGON ((6.385532 49.83703...
## 7 POLYGON ((6.316665 49.62337...
## 8 POLYGON ((6.425158 49.73164...
## 9 POLYGON ((5.998312 49.69992...
## 10 POLYGON ((6.039474 49.44826...and the vect() method to get from sf to terra:
v_sf_rt <- vect(v_sf)
v_sf_rt## class : SpatVector
## geometry : polygons
## dimensions : 12, 6 (geometries, attributes)
## extent : 5.74414, 6.528252, 49.44781, 50.18162 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## names : ID_1 NAME_1 ID_2 NAME_2 AREA POP
## type : <num> <chr> <num> <chr> <num> <int>
## values : 1 Diekirch 1 Clervaux 312 18081
## 1 Diekirch 2 Diekirch 218 32543
## 1 Diekirch 3 Redange 259 18664all.equal(v_sf_rt, v, check.attributes=FALSE)## [1] TRUE"Spatial"
To coerce to and from vector classes defined in the sp package (Roger S. Bivand, Pebesma, and Gomez-Rubio 2013), use the classes defined in sf as an intermediate step:
v_sp <- as(v_sf, "Spatial")
print(summary(v_sp))## Object of class SpatialPolygonsDataFrame
## Coordinates:
## min max
## x 5.74414 6.528252
## y 49.44781 50.181622
## Is projected: FALSE
## proj4string : [+proj=longlat +datum=WGS84 +no_defs]
## Data attributes:
## ID_1 NAME_1 ID_2 NAME_2
## Min. :1.000 Length:12 Min. : 1.00 Length:12
## 1st Qu.:1.000 Class :character 1st Qu.: 3.75 Class :character
## Median :2.000 Mode :character Median : 6.50 Mode :character
## Mean :1.917 Mean : 6.50
## 3rd Qu.:3.000 3rd Qu.: 9.25
## Max. :3.000 Max. :12.00
## AREA POP
## Min. : 76.0 Min. : 5163
## 1st Qu.:187.2 1st Qu.: 18518
## Median :225.5 Median : 26097
## Mean :213.4 Mean : 50167
## 3rd Qu.:253.0 3rd Qu.: 36454
## Max. :312.0 Max. :182607v_sp_rt <- vect(st_as_sf(v_sp))
all.equal(v_sp_rt, v, check.attributes=FALSE)## [1] TRUE"SpatRaster" coercion
In the terra package, raster data are held in "SpatRaster" objects. This means that when read_RAST() is used, a "SpatRaster" object is returned, and the same class of object is needed for write_RAST() for writing to GRASS.
fr <- system.file("ex/elev.tif", package="terra")
(r <- rast(fr))## class : SpatRaster
## dimensions : 90, 95, 1 (nrow, ncol, nlyr)
## resolution : 0.008333333, 0.008333333 (x, y)
## extent : 5.741667, 6.533333, 49.44167, 50.19167 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source : elev.tif
## name : elevation
## min value : 141
## max value : 547In general, "SpatRaster" objects are files, rather than data held in memory:
try(inMemory(r))## [1] FALSE"stars"
The stars package (Pebesma 2021) uses GDAL through sf. A coercion method is provided from "SpatRaster" to "stars":
r_stars <- st_as_stars(r)
print(r_stars)## stars object with 2 dimensions and 1 attribute
## attribute(s):
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## elev.tif 141 291 333 348.3366 406 547 3942
## dimension(s):
## from to offset delta refsys point values x/y
## x 1 95 5.74167 0.00833333 WGS 84 FALSE NULL [x]
## y 1 90 50.1917 -0.00833333 WGS 84 FALSE NULL [y]which round-trips in memory.
(r_stars_rt <- rast(r_stars))## class : SpatRaster
## dimensions : 90, 95, 1 (nrow, ncol, nlyr)
## resolution : 0.008333333, 0.008333333 (x, y)
## extent : 5.741667, 6.533333, 49.44167, 50.19167 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source : memory
## name : lyr.1
## min value : 141
## max value : 547When coercing to "stars_proxy" the same applies:
(r_stars_p <- st_as_stars(r, proxy=TRUE))## stars_proxy object with 1 attribute in 1 file(s):
## $elev.tif
## [1] "[...]/elev.tif"
##
## dimension(s):
## from to offset delta refsys point values x/y
## x 1 95 5.74167 0.00833333 WGS 84 FALSE NULL [x]
## y 1 90 50.1917 -0.00833333 WGS 84 FALSE NULL [y]with coercion from "stars_proxy" also not reading data into memory:
(r_stars_p_rt <- rast(r_stars_p))## class : SpatRaster
## dimensions : 90, 95, 1 (nrow, ncol, nlyr)
## resolution : 0.008333333, 0.008333333 (x, y)
## extent : 5.741667, 6.533333, 49.44167, 50.19167 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source : elev.tif
## name : elevation
## min value : 141
## max value : 547"Spatial"
Coercion to the sp "SpatialGridDataFrame" representation is provided by stars:
r_sp <- as(r_stars, "Spatial")
summary(r_sp)## Object of class SpatialGridDataFrame
## Coordinates:
## min max
## x 5.741667 6.533333
## y 49.441667 50.191667
## Is projected: FALSE
## proj4string : [+proj=longlat +datum=WGS84 +no_defs]
## Grid attributes:
## cellcentre.offset cellsize cells.dim
## x 5.745833 0.008333333 95
## y 49.445833 0.008333333 90
## Data attributes:
## elev.tif
## Min. :141.0
## 1st Qu.:291.0
## Median :333.0
## Mean :348.3
## 3rd Qu.:406.0
## Max. :547.0
## NA's :3942and can be round-tripped:
(r_sp_rt <- rast(st_as_stars(r_sp)))## class : SpatRaster
## dimensions : 90, 95, 1 (nrow, ncol, nlyr)
## resolution : 0.008333333, 0.008333333 (x, y)
## extent : 5.741667, 6.533333, 49.44167, 50.19167 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source : memory
## name : lyr.1
## min value : 141
## max value : 547"RasterLayer"
As rgdal will retire soon, using rgdal in the raster package (Hijmans 2022a) is not a feasible way to coerce a "SpatRaster" object to a "RasterLayer" object.
tf <- tempfile(fileext=".grd")
terra::writeRaster(r, filename=tf, filetype="RRASTER")We can use a temporary file in the "RRASTER" format that terra::writeRaster() can write using GDAL, and raster::raster() can read using R code rather than GDAL:
(r_RL <- raster(tf))## class : RasterLayer
## dimensions : 90, 95, 8550 (nrow, ncol, ncell)
## resolution : 0.008333333, 0.008333333 (x, y)
## extent : 5.741667, 6.533333, 49.44167, 50.19167 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : file2079d36daa058.grd
## names : elevation
## values : 0, 547 (min, max)This object (held on file rather than in memory) can be round-tripped:
(r_RL_rt <- rast(r_RL))## class : SpatRaster
## dimensions : 90, 95, 1 (nrow, ncol, nlyr)
## resolution : 0.008333333, 0.008333333 (x, y)
## extent : 5.741667, 6.533333, 49.44167, 50.19167 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs
## source : file2079d36daa058.grd
## name : elevation
## min value : 0
## max value : 547