`readRAST` and `writeRAST` support for `raster` and `terra` #42

bniebuhr · 2022-02-02T23:28:05Z

Hi,

Is there any plan to support reading and writing rasters to/from R using the raster and terra packages?

For instance, if I have to read a raster from GRASS GIS and I want to use terra within R, so far I need the following:

library(rgrass7)
library(raster)
library(terra)
library(dplyr)

my_map <- rgrass7::readRAST("my_map_name_within_GRASS") %>%
    raster::raster() %>%
    terra::rast()

Even though I could drop the use of dplyr here, that is quite some code for a simple reading command.
Any support in this direction would be a super nice feature!

The text was updated successfully, but these errors were encountered:

rsbivand · 2022-02-03T09:30:19Z

Indeed, dplyr etc. are completely irrelevant, and should be avoided, as should non-native pipes, because they may create extra copies of large objects. For now, simply coerce and ask terra to provide a method from SpatialGridDataFrame to SpatRaster directly.

Support for terra might be provided if someone conntributes a full PR, for example using terra to read the temporary file written as part of readRAST(). Any PR should be for the rgrass7 branch of the rgrass github repository.

It is more likely that support for stars will emerge, but both really need a proxy approach, leaving the data in GRASS without bringing it into the R workspace. However, the GDAL plugin driver for GRASS is not under active maintenance. For this reason there are no specific plans now.

bniebuhr · 2022-02-03T10:23:23Z

Thanks, @rsbivand . I'll raise an issue at the terra repo, then, at least for the time being.
In any case, I keep my comment here that, in my view, further integration directly within readRAST() would ease the integration R-GRASS (whether it is via raster, terra, or stars).

rsbivand · 2022-02-03T11:14:14Z

Yes, but all are welcome to contribute code. I am retired, and will cease supporting this package soon. Contribute a PR, please. Consider such contributions as a reasonable charge on your time.

rhijmans · 2022-02-03T21:43:02Z

Do you refer to the temporary file that is produced with rgrass7:::.read_rast_non_plugin ? There I see

rOutBinFlags <- "b"
execGRASS("r.out.bin", flags = rOutBinFlags,  input = vname[i], output = gtmpfl11, bytes = rOutBinBytes, 
                  null = as.integer(NODATA), ignore.stderr = ignore.stderr)`

And according to the GRASS manual that creates a flat binary with hdr and world files. Surely raster/terra/stars can read that file?

rsbivand · 2022-02-04T10:58:47Z

Currently read by

rgrass/R/bin_link.R

Lines 283 to 288 in c10d347

    
           readBinGridData <- function(fname, what, n, size, endian, nodata) { 
        
           	map <- readBin(fname, what=what, n=n, size=size, signed=TRUE, 
        
           		endian=endian) 
        
           	is.na(map) <- map == nodata 
        
           	map 
        
           }

, and post-processed to SpatialGridDataFrame for all rasters imported in the same call. Flat binary was found by Rainer Krug @rkrug to be much faster than export to GTiff and using GDAL GTiff (probably unneeded compression, and custom processing permitted multiple bands to be handled in one call (GRASS only has image objects as multiband).

So not simple to unpick, someone should contribute a back-end perhaps within readRAST() to create a "SpatRaster" or appropriate multiband/multilayer representation, postprocessing strings as factors, and integers and numerics appropriately.

Since "SpatRaster" can (?) handle integers other than signed 4-byte, and numeric other than 8-byte, I think this isn't simple, and that the PR should come with tinytest tests that run when a specific GRASS location is present.

rhijmans · 2022-02-04T17:21:30Z

Could readRast have an option to return the filename(s) instead? I would expect you could just to do rast(filename) with most files, and where not, the burden to deal with it would be somewhere elseBI believe that both can read various signed and unsigned integer (8,16,32,64 bytes) except GDAL cannot read signed Byte and perhaps raster does not read 64 Ints. Both read 4 and 8 byte numeric (float and double).

I wondered if readRAST(... return_SGDF=FALSE) already does that; based on the guess that "a list of bands" refers to filenames; I realize that this may not be the case. So I ventured into rgrass for the first time and tried:

library(raster)
f <- system.file("ex/elev.tif", package="terra")
r <- raster(f)
g <- as(r, "SpatialGridDataFrame")

library(rgrass7)
initGRASS("c:/OSGeo4W/apps/grass/grass78")
writeRAST(g, "test", zcol="elev", flags=c("quiet", "overwrite"))
execGRASS("r.info", map="test")
# +----------------------------------------------------------------------------+
# | Map:      test                           Date: Fri Feb 04 09:01:30 2022    |

x <- readRAST(c("test", "elev"),  useGDAL=FALSE, return_SGDF=FALSE)
#Error in sp::CRS(getLocationProj()) : NA

I first updated the package from CRAN, and then from github, and in both cases I get:

writeRAST(g, "test", zcol="elev", flags=c("quiet", "overwrite"))
#Error in writeRAST(g, "test", zcol = "elev", flags = c("quiet", "overwrite")) :
#  no stars support yet

class(g)
#[1] "SpatialGridDataFrame"
#attr(,"package")
#[1] "sp"

rsbivand · 2022-02-04T17:38:34Z

Either use_sp() or use_sf() is needed (the latter to use sf for vector file transfer).

rhijmans · 2022-02-04T17:42:34Z

Ay, sorry, missed that after the update! Now getting this error again:

x <- readRAST(c("test", "elev"),  useGDAL=FALSE, return_SGDF=FALSE)
#Error in sp::CRS(getLocationProj()) : NA

rsbivand · 2022-02-04T17:43:17Z

Yes, seeing the same here:

Browse[3]> getLocationProj()
[1] "XY location (unprojected)"

so sp::CRS() fails. GRASS locations have a defined CRS.

 execGRASS("g.proj", flags="p")
XY location (unprojected)
> execGRASS("g.proj", georef=f, flags=c("j", "f"))
+proj=longlat +datum=WGS84 +no_defs +type=crs

rsbivand · 2022-02-04T18:14:46Z

This works - the problem was that GRASS locations must have a defined CRS, but initGRASS() doesn't set one automatically. So a number of extra steps are needed when using a temporary location. When R is started and the package loaded inside an existing location, all this does happen transparently.

library(raster)
f <- system.file("ex/elev.tif", package="terra")
r <- raster(f)
g <- as(r, "SpatialGridDataFrame")
library(rgrass7)
use_sp()
loc <- initGRASS("/home/rsb/topics/grass/g800/grass80", SG=g, override=TRUE)
execGRASS("g.mapset", mapset="PERMANENT", flag="quiet")
execGRASS("g.proj", flag="c", epsg=4326)
execGRASS("g.mapset", mapset=loc$MAPSET, flag="quiet")
execGRASS("g.region", flag="d")
writeRAST(g, "test", zcol="elev", flags=c("quiet", "overwrite"))
execGRASS("g.list", type="raster")
x <- readRAST("test",  useGDAL=FALSE, return_SGDF=FALSE)
str(x)

rhijmans · 2022-02-04T19:49:35Z

Thanks very much and apologies for stumbling over these novice issues. I see that return_SGDF=FALSE does not what I was hoping.

I added an argument return_filenames to my fork ( here: rhijmans@cca115d )
(I add the extension ".bil" to the files but I do not think that is required).

This appears to work fine for this example, insofar that I get the same result as with return_SGDF=TRUE which in this example is not correct (only 1 row and 1 column); but I am guessing that is a different matter. I have no idea about other cases or other implications. For example, is it likely that these files will get overwritten?

library(raster)
#Loading required package: sp
f <- system.file("ex/elev.tif", package="terra")
r <- raster(f)
g <- as(r, "SpatialGridDataFrame")
library(rgrass7)
#Loading required package: XML
#GRASS GIS interface loaded with GRASS version: (GRASS not running)
use_sp()
loc <- initGRASS("c:/OSGeo4W/apps/grass/grass78")

execGRASS("g.mapset", mapset="PERMANENT", flag="quiet")
#WARNING: Concurrent mapset locking is not supported on Windows
execGRASS("g.proj", flag="c", epsg=4326)
#Default region was updated to the new projection, but if you have multiple
#mapsets `g.region -d` should be run in each to update the region from the
#default
#Projection information updated
execGRASS("g.mapset", mapset=loc$MAPSET, flag="quiet")
#WARNING: Concurrent mapset locking is not supported on Windows
execGRASS("g.region", flag="d")

writeRAST(g, "test", zcol="elev", flags=c("quiet", "overwrite"))
execGRASS("r.info", map="test")
# +----------------------------------------------------------------------------+
# | Map:      test                           Date: Fri Feb 04 11:32:48 2022    |
# | Mapset:   file262825a979cc               Login of Creator: rhijm           |
# | Location: file26286ec84517                                                 |
# | DataBase: c:/temp/RtmpMjFdbD                                               |
# | Title:                                                                     |
# | Timestamp: none                                                            |
# |----------------------------------------------------------------------------|
# |                                                                            |
# |   Type of Map:  raster               Number of Categories: 0               |
# |   Data Type:    CELL                                                       |
# |   Rows:         90                                                         |
# |   Columns:      95                                                         |
# |   Total Cells:  8550                                                       |
# |        Projection: Latitude-Longitude                                      |
# |            N: 50:11:30.12N    S: 49:26:30.12N   Res: 0:00:30               |
# |            E: 6:31:59.88E    W: 5:44:30.12E   Res: 0:00:29.997474          |
# |   Range of data:    min = 141  max = 547                                   |
# |                                                                            |
# |   Data Description:                                                        |
# |    generated by RINBIN~1                                                   |
# |                                                                            |
# |   Comments:                                                                |
# |    RINBIN~1 --overwrite --quiet input="c:/temp/RtmpMjFdbD/file26286ec84\   |
# |    517/file262825a979cc/.tmp/unknown/X581" output="test" bytes=4 header\   |
# |    =0 bands=1 order="native" north=50.1917 south=49.4417 east=6.5333 we\   |
# |    st=5.7417 rows=90 cols=95 anull=140                                     |
# |                                                                            |
# +----------------------------------------------------------------------------+


x <- readRAST("test",  useGDAL=FALSE, return_SGDF=TRUE)
#Creating BIL support files...
#Exporting raster as integer values (bytes=4)
# 100%
#Warning message:
#In system(syscmd, intern = intern, ignore.stderr = ignore.stderr,  :
#  running command 'g.proj.exe -w' had status 884

x
#Object of class SpatialGridDataFrame
#Object of class SpatialGrid
#Grid topology:
#  cellcentre.offset cellsize cells.dim
#1               0.5        1         1
#2               0.5        1         1
#SpatialPoints:
#      s1  s2
#[1,] 0.5 0.5
#Coordinate Reference System (CRS) arguments: +proj=longlat +datum=WGS84
#+no_defs
#
#Data summary:
#   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's
#     NA      NA      NA     NaN      NA      NA       1
#Warning message:
#In proj4string(x) :
#  CRS object has comment, which is lost in output; in tests, see
#https://cran.r-project.org/web/packages/sp/vignettes/CRS_warnings.html

fout <- readRAST("test",  useGDAL=FALSE, return_filenames=TRUE)
#Creating BIL support files...
#Exporting raster as integer values (bytes=4)
# 100%
#Warning message:
#In system(syscmd, intern = intern, ignore.stderr = ignore.stderr,  :
#  running command 'g.proj.exe -w' had status 884

raster(fout)
#class      : RasterLayer
#dimensions : 1, 1, 1  (nrow, ncol, ncell)
#resolution : 1, 1  (x, y)
#extent     : 0, 1, 0, 1  (xmin, xmax, ymin, ymax)
#crs        : NA
#source     : test.bil
#names      : test

rsbivand · 2022-02-04T20:25:52Z

I just pushed some edits to rgrass@rgrass7, that is the rgrass7 branch. With these, the location CRS is set from the WKT of the SG= "SpatialGrid" object.

library(raster)
f <- system.file("ex/elev.tif", package="terra")
r <- raster(f)
g <- as(r, "SpatialGridDataFrame")
library(rgrass7)
use_sp()
loc <- initGRASS("/home/rsb/topics/grass/g800/grass80", SG=g, override=TRUE)
writeRAST(g, "test", zcol="elev", flags=c("quiet", "overwrite"))
execGRASS("g.list", type="raster")
x <- readRAST("test",  useGDAL=FALSE, return_SGDF=FALSE)
str(x)

Yes, it should be possible to return the temporary file name as you suggest for numerical layers. A lot of the extra code is special casing systems with odd paths, like cygwin, added years ago when it was commonly used. However, I think that:

tf <- tempfile()
execGRASS("r.out.bin", flags=c("i", "b"), input="test", output=tf)
raster(tf)

is simpler and more idiomatic - the user would have to choose "i" or "f". terra::rast(tf) also works in the same setting. readRAST() etc. were written long before execGRASS(), @bniebuhr does this help?

rsbivand · 2022-02-04T20:29:14Z

I probably also should look at whether GRASS 8 supports the instantiation of a location from a file name, rather than needing to read it into R first.

rhijmans · 2022-02-04T23:18:04Z

Applying the same changes to the "rgrass7" branch (rhijmans@96f9d58), and using your code (except for the GRASS path), I now get

library(raster)
f <- system.file("ex/elev.tif", package="terra")
r <- raster(f)
g <- as(r, "SpatialGridDataFrame")
library(rgrass7)
use_sp()
loc <- initGRASS("c:/OSGeo4W/apps/grass/grass78", SG=g, override=TRUE)
writeRAST(g, "test", zcol="elev", flags=c("quiet", "overwrite"))
execGRASS("g.list", type="raster")

x <- readRAST("test",  useGDAL=FALSE, return_SGDF=TRUE)
str(x)
#Formal class 'SpatialGridDataFrame' [package "sp"] with 4 slots
#  ..@ data       :'data.frame': 1 obs. of  1 variable:
#  .. ..$ test: int NA
#  ..@ grid       :Formal class 'GridTopology' [package "sp"] with 3 slots
#  .. .. ..@ cellcentre.offset: num [1:2] 0.5 0.5
#  .. .. ..@ cellsize         : num [1:2] 1 1
#  .. .. ..@ cells.dim        : int [1:2] 1 1
#  ..@ bbox       : num [1:2, 1:2] 0 0 1 1
#  .. ..- attr(*, "dimnames")=List of 2
#  .. .. ..$ : NULL
#  .. .. ..$ : chr [1:2] "min" "max"
#  ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
#  .. .. ..@ projargs: chr "+proj=longlat +datum=WGS84 +no_defs"
#  .. .. ..$ comment: chr "GEOGCRS[\"unknown\",\n    DATUM[\"World Geodetic System 1984\",\n        ELLIPSOID[\"WGS 84\",6378137,298.25722"| __truncated__

library(terra)
rast(x)

tmpf <- readRAST("test",  useGDAL=FALSE, return_SGDF=FALSE, return_filenames=TRUE)
rast(tmpf)

rast(tmpf)
#class       : SpatRaster
#dimensions  : 1, 1, 1  (nrow, ncol, nlyr)
#resolution  : 1, 1  (x, y)
#extent      : 0, 1, 0, 1  (xmin, xmax, ymin, ymax)
#coord. ref. : lon/lat WGS 84
#source      : test
#name        : test

So, it seems that the filename method works, but that there is still something else that goes wrong because the BIL file only has a single cell.

rsbivand · 2022-02-05T12:46:09Z

I haven't checked, but readRAST() unlinks the temporary file(s) - one for each layer/variable in SGDF - on exit.

rhijmans · 2022-02-05T17:51:25Z

The files are there, I took care of that so in that sense it "works". But note that the SGDF is the same. It only has once grid cell and that is not correct.

veroandreo · 2022-02-08T12:18:10Z

The files are there, I took care of that so in that sense it "works". But note that the SGDF is the same. It only has once grid cell and that is not correct.

GRASS will export rasters with the resolution and extent of the computational region. When you create a new location, say a latlon one with EPSG 4326, the default region will be of one cell (one degree, which is the unit of the CRS). If it were a metric CRS, the one cell would be one meter.

To get from GRASS the same raster that you imported into it, you need to set the computational region to the extent and res of your raster first: execGRASS("g.region", raster="test", flags="p"); after that the readRAST() call should give you the same raster that you imported. This happens only because you created a new location, when R is called within a location or you open GRASS from R in an existing location, it automatically reads whatever region is set.

HTH :)

rhijmans · 2022-02-09T05:21:09Z

Thanks @veroandreo, now I got it to work:

library(raster)
f <- system.file("ex/elev.tif", package="terra")
r <- raster(f)
g <- as(r, "SpatialGridDataFrame")
library(rgrass7)
use_sp()
loc <- initGRASS("c:/OSGeo4W/apps/grass/grass78", SG=g, override=TRUE)
writeRAST(g, "test", zcol="elev", flags=c("quiet", "overwrite"))
execGRASS("g.region", raster="test", flags="p")
tmpf <- readRAST("test",  useGDAL=FALSE, return_SGDF=FALSE, return_filenames=TRUE)

library(terra)
x <- rast(tmpf)
x
#class       : SpatRaster
#dimensions  : 90, 95, 1  (nrow, ncol, nlyr)
#resolution  : 0.008333, 0.008333  (x, y)
#extent      : 5.7417, 6.533335, 49.44173, 50.1917  (xmin, xmax, ymin, ymax)
#coord. ref. : lon/lat WGS 84
#source      : test
#name        : test

By having readRAST (or another function) return a filename or filenames, it works well with raster/terra/stars/sp.
If there is important additional metadata that is not written the the hdr/wrld files, than perhaps returning a list would be an option.

rsbivand · 2022-02-09T09:14:11Z

Yes but not at all is this simple. Rasters with category labels are harder than numeric rasters, so any implementation has to handle them.

The nc_basic_spm_grass7 geology layer based on digitized https://files.nc.gov/ncdeq/Energy+Mineral+and+Land+Resources/Geological+Survey/1985_state_geologic_map_500000_scale.pdf has multiple instances of geology labels. The current code to recreate the raster is execGRASS("v.to.rast", input="geology", output="geology_10m", use="attr", type="point,line,area", attribute_column="GEOL250_ID", label_column="GEO_NAME") rather than the GRASS 6 (?) version in execGRASS("r.info", map="geology"). This makes the RAT from r.out.* odd - in the original readRAST() implementation, it was assumed that the category numbers (matching the raster values) were meaningful, wrongly - the labels are correct. So:

tf0 <- tempfile(fileext=".tif")
execGRASS("r.out.gdal", input="geology", output=tf0, format="GTiff", flags=c("overwrite", "c", "m", "t"))
geol <- rast(tf0)
str(levels(geol))
# List of 1
# $ : chr [1:1832] "Zml" "Zmf" "Zml" "Zml" ...
plot(geol)

is correct, and the original implementation is wrong. I can't read this *.grd correctly with read_stars(), and coercion from "SpatRaster" is unsuccessful in the same way (incorrect factor label assignment). The RAT contains all of the 1832 digitized geology polygons despite the extent of the raster being very much smaller than that of the vector.

This also works:

tf0 <- tempfile(fileext=".grd")
execGRASS("r.out.gdal", input="geology", output=tf0, format="RRASTER", flags=c("overwrite", "c", "t"))
geol <- rast(tf0)
str(levels(geol))
# List of 1
# $ : chr [1:1832] "Zml" "Zmf" "Zml" "Zml" ...
plot(geol)

with the raster-native format, but with slightly changed flags.

It would be interesting to see whether the "RRASTER" driver could be the workhorse using r.out.gdal, rather than r.out.bin (and from R to GRASS), because it seems to support RAT. It also supports multiple bands of the same kind. The readRAST() legacy implementation supported bands of different kinds (float/integer/factor) because the R representation was a data.frame, by looping over the GRASS rasters to be read in the current region resolution and extent.

rhijmans · 2022-02-14T22:37:03Z

RRASTER could be an option, I think. But If using r.out.gdal can't the output format be GTiff? (I think I saw something about r.out.bin being used for speed; but are there big differences between GDAL drivers?). GDAL also has the generic PAM XML approach for RATs. I do not think it would be terribly difficult to generate these with R, if r.out.bin is to be preferred. r.out.gdal might add them when writing to the equivalent (?) GenBin format.

rsbivand · 2022-02-14T23:01:10Z

I'll try to push a branch tomorrow with a skeleton implementation for GRASS to R. So far SpatRaster is well-behaved wrt. to cats, based on RRASTER and the RAT representation. I use the GRASS nc basic location for finding corner cases.

rsbivand · 2022-02-15T21:21:06Z

The rastNG branch has been pushed to github. The data used for trying things out may be installed in the working directory like this:

url <- "https://grass.osgeo.org/sampledata/north_carolina/nc_basic_spm_grass7.zip"
to0 <- getOption("timeout")
options(timeout = max(300, to0))
download.file(url, "nc_basic_spm_grass7.zip")
unzip("nc_basic_spm_grass7.zip")
options(timeout = to0)

Using this standard location but subsetting the region, we can start GRASS in R:

library(rgrass7)
my_GRASS <- "/home/rsb/topics/grass/g800/grass80"
loc <- initGRASS(my_GRASS, home=tempdir(), gisDbase=".", 
 location="nc_basic_spm_grass7", mapset="user1", override=TRUE)
execGRASS("g.region", raster="elevation")
execGRASS("g.region", flags="p")

Reading numeric rasters from GRASS to R seems OK:

execGRASS("r.slope.aspect", flags="overwrite", elevation="elevation",
 slope="slope", aspect="aspect")
use_sp()
u00 <- readRAST(c("elevation", "slope", "aspect")) # legacy
object.size(u00) # in memory
u0 <- rgrass7:::read_RAST(c("elevation", "slope", "aspect"),
 return_format="SGDF") # re-written SGDF
object.size(u0) # in memory
library(sp)
all.equal(u00[["elevation"]], u0[["elevation"]])
u1 <- rgrass7:::read_RAST(c("elevation", "slope", "aspect"),
 return_format="terra") # re-written terra
object.size(u1) # on file
library(terra)
u1
all.equal(u00[["elevation"]], c(values(u1[["elevation"]])))
library(stars)
u2 <- st_as_stars(u1) # coerce to stars
u2
object.size(u2) # in memory
all.equal(u00[["elevation"]], c(dplyr::pull(u2[,,,"elevation"])))
u2a <- read_stars(u1@ptr$filenames, proxy=TRUE)
u2a
object.size(u2a) # on file
all.equal(u00[["elevation"]], c(st_as_stars(u2a)[[1]])) # moved to memory

Here, the "SGDF" output re-uses legacy code (without mismatched category handling), and the temorary files are unlinked on return. The "terra" output is new. Rather than add a "stars" output format, the user may coerce - typically, "terra" holds data on file (not unlinked), "SGDF" as before in memory, as is "stars" by cercion. It is possible to make a "stars_proxy" object as shown, but computing on the values means moving to memory.

Categories are perhaps more difficult, only some have correctly matched codes, labels and colours.

# landuse categories and colors (1-to-1 match between codes, labels and colors
v0 <- rgrass7:::read_RAST("landuse", cat=TRUE, return_format="SGDF")
cl0 <- rgrass7:::read_cat_colors("landuse")
x11()
spplot(v0, "landuse", col.regions=cl0$pal, main="SGDF")
v1 <- rgrass7:::read_RAST("landuse", cat=TRUE, return_format="terra")
x11()
plot(v1[["landuse"]], col=cl0$pal, main="terra")
v2 <- st_as_stars(v1)
x11()
plot(v2, col=cl0$pal)
library(tmap)
x11()
tm_shape(v1) + tm_raster(pal=cl0$pal, title="Land use") +
 tm_layout(legend.outside=TRUE, title="tmap")

The "geology" raster is harder, as there is one code for each polygon digitized, and non-unique labels are recorded for each polygon. The polygon count covered the whole state, so the raster attribute table provided by default through GDAL includes many unused entries, and multiple codes for the same label. The legacy readRAST() function always got this wrong, because it assumed that the codes, not the labels, should be preserved, and catenated code_label strings were returned as factor levels. An excerpt from the state-wide geological map for SW Wade County is upper left in the screenshot below.

# raster from v.to.rast from scanned polygons from:
# https://files.nc.gov/ncdeq/Energy+Mineral+and+Land+Resources/Geological+Survey/1985_state_geologic_map_500000_scale.pdf
v0 <- rgrass7:::read_RAST("geology", cat=TRUE, return_format="SGDF")
# warning about categories - polygons with same category label have different
# codes, multipolygons split between several codes
ci <- rgrass7:::cats_internals("geology") # attempt to reclassify in GRASS
a <- split(ci$rcp_out[,1], ci$rcp_out[,2])
tf <- tempfile()
write.table(data.frame(input=sapply(a, paste, collapse=" "), eq="=", ci$levs),
 file=tf, row.names=FALSE, col.names=FALSE, quote=FALSE)
file.show(tf)
execGRASS("r.reclass", input="geology", output="geology_rc", rules=tf,
 flags="overwrite") # reclassifying in GRASS
v0a <- rgrass7:::read_RAST("geology_rc", TRUE, return_format="SGDF")
# palette taken from source legend colors, subset SW Wake County (Raleigh)
pal <- c("#73A4CF", "#6F87A9", "#ABA9B9", "#5D8DC3", "#BC3F28", "#6D2449",
 "#DA8B9B", "#C6C41B")
x11()
spplot(v0a, "geology_rc", col.regions=pal)

sp output upper right.

It appears that "terra" handles the RAT well (lower left):

v1 <- rgrass7:::read_RAST("geology", cat=TRUE, return_format="terra")
str(levels(v1))
x11()
plot(v1, col=pal)# terra handles multiple codes correctly

This output agrees with reading the reclassified GRASS raster, and with reclassification using methods in terra:

v1a <- rgrass7:::read_RAST("geology_rc", cat=TRUE, return_format="terra")
str(levels(v1a))
x11()
plot(v1a, col=pal)
v1r <- rgrass7:::repair_cats(v1, "geology", "geology")# reclassifying in terra
str(levels(v1r))
x11()
plot(v1r, col=pal)
x11()

Problems arise with "stars" coercion in tmap (lower right output in screenshot):

# tmap calls stars::st_as_stars() internally
# fails to handle code/label mismatch
tm_shape(v1, raster.downsample=FALSE) + tm_raster(pal=pal, title="Geology") +
 tm_layout(legend.outside=TRUE, title="tmap")
# offset code/label/key
x11()
tm_shape(v1a, raster.downsample=FALSE) + tm_raster(pal=pal, title="Geology") +
 tm_layout(legend.outside=TRUE, title="tmap")
# offset code/label/key
x11()
v1as <- st_as_stars(v1a)
plot(v1as, col=pal)
x11()
tm_shape(v1as, raster.downsample=FALSE) + tm_raster(pal=pal, title="Geology") +
 tm_layout(legend.outside=TRUE, title="tmap")
# offset code/label/key
x11()
tm_shape(v1r, raster.downsample=FALSE) + tm_raster(pal=pal, title="Geology") +
 tm_layout(legend.outside=TRUE, title="tmap")

Since the coercion method is not heavily exercised yet, and this is a corner case, coercion for numeric rasters and categorical rasters with 1-to-1 code/label relationships should be OK. After more testing, perhaps an issue can be raised on stars.

rsbivand · 2022-02-17T19:04:09Z

The rastNG branch now also has write_RAST() and initGRASS() adapted for "SpatRaster" objects:

library(terra)
f <- system.file("ex/elev.tif", package="terra")
r <- rast(f)
plot(r, col=grDevices::terrain.colors(50))
library(rgrass7)
my_GRASS <- "/home/rsb/topics/grass/g801/grass80"
loc <- initGRASS(my_GRASS, home=tempdir(), SG=r, override=TRUE)
rgrass7:::write_RAST(r, "elev", flags="overwrite")
execGRASS("r.info", map="elev")
execGRASS("r.slope.aspect", flags="overwrite", elevation="elev", slope="slope", aspect="aspect")
u1 <- rgrass7:::read_RAST(c("elev", "slope", "aspect"), return_format="terra")
plot(u1[["elev"]], col=grDevices::terrain.colors(50))

Comments please @bniebuhr ! For now use rgrass7::: to access the new functions. "SpatVector" to follow, giving a route to freeing the package from rgdal, which needed doing anyway. After that, temporal needs attention.

rsbivand · 2022-02-18T16:11:39Z

And "SpatVector":

v1 <- rgrass7:::read_VECT("census")
v1
rgrass7:::write_VECT(v1, "census_sV")
execGRASS("v.info", map="census_sV")

rsbivand · 2022-03-05T10:49:58Z

rgrass7_0.2-8 submitted to CRAN.

rsbivand · 2022-03-05T17:13:20Z

rgrass7_0.2-8 on CRAN.

bniebuhr · 2022-03-06T22:20:04Z

@rsbivand , that is fantastic!
Sorry for not participating actively after the initial suggestion, but at some point the level of discussion between you and @rhijmans was above my understanding of both R and GRASS. But now I just tested everything and works smoothly. I see you guys have throughly tested it as well. I am happy to see this unfolding of this initial suggestion. Thanks for all this hard work!

I followed the discussion in PR #45, but just to make sure if I understand. Now using read_RAST(), for instance, do I still need use_sp() before it?
Besides, is there a particular reason for keeping "SGDF" as the standard value for the return_format parameter? Wouldn't it make sense to use "terra" instead as the default?

rsbivand · 2022-03-06T22:33:57Z

I do not think use_sp() is needed. However, this is a process, so maybe wating for reactions to changes so far is fair to legacy users.

bniebuhr · 2022-03-07T03:17:21Z

Really fair, indeed!

bniebuhr mentioned this issue Feb 3, 2022

Support to read rspatial/terra#524

Closed

rsbivand mentioned this issue Feb 18, 2022

FOSS4G 2022 r-spatial workshops ? r-spatial/discuss#53

Open

rsbivand mentioned this issue Feb 18, 2022

Rast ng probably ready to merge #45

Merged

11 tasks

rsbivand closed this as completed Mar 5, 2022

rsbivand mentioned this issue Mar 18, 2022

Change terra intermediate file format? #51

Closed

4 tasks

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

`readRAST` and `writeRAST` support for `raster` and `terra` #42

`readRAST` and `writeRAST` support for `raster` and `terra` #42

bniebuhr commented Feb 2, 2022

rsbivand commented Feb 3, 2022

bniebuhr commented Feb 3, 2022

rsbivand commented Feb 3, 2022

rhijmans commented Feb 3, 2022 •

edited

rsbivand commented Feb 4, 2022

rhijmans commented Feb 4, 2022 •

edited

rsbivand commented Feb 4, 2022

rhijmans commented Feb 4, 2022

rsbivand commented Feb 4, 2022 •

edited

rsbivand commented Feb 4, 2022

rhijmans commented Feb 4, 2022

rsbivand commented Feb 4, 2022

rsbivand commented Feb 4, 2022

rhijmans commented Feb 4, 2022

rsbivand commented Feb 5, 2022

rhijmans commented Feb 5, 2022

veroandreo commented Feb 8, 2022

rhijmans commented Feb 9, 2022

rsbivand commented Feb 9, 2022 •

edited

rhijmans commented Feb 14, 2022

rsbivand commented Feb 14, 2022

rsbivand commented Feb 15, 2022 •

edited

rsbivand commented Feb 17, 2022

rsbivand commented Feb 18, 2022

rsbivand commented Mar 5, 2022

rsbivand commented Mar 5, 2022

bniebuhr commented Mar 6, 2022

rsbivand commented Mar 6, 2022

bniebuhr commented Mar 7, 2022

readRAST and writeRAST support for raster and terra #42

readRAST and writeRAST support for raster and terra #42

Comments

bniebuhr commented Feb 2, 2022

rsbivand commented Feb 3, 2022

bniebuhr commented Feb 3, 2022

rsbivand commented Feb 3, 2022

rhijmans commented Feb 3, 2022 • edited

rsbivand commented Feb 4, 2022

rhijmans commented Feb 4, 2022 • edited

rsbivand commented Feb 4, 2022

rhijmans commented Feb 4, 2022

rsbivand commented Feb 4, 2022 • edited

rsbivand commented Feb 4, 2022

rhijmans commented Feb 4, 2022

rsbivand commented Feb 4, 2022

rsbivand commented Feb 4, 2022

rhijmans commented Feb 4, 2022

rsbivand commented Feb 5, 2022

rhijmans commented Feb 5, 2022

veroandreo commented Feb 8, 2022

rhijmans commented Feb 9, 2022

rsbivand commented Feb 9, 2022 • edited

rhijmans commented Feb 14, 2022

rsbivand commented Feb 14, 2022

rsbivand commented Feb 15, 2022 • edited

rsbivand commented Feb 17, 2022

rsbivand commented Feb 18, 2022

rsbivand commented Mar 5, 2022

rsbivand commented Mar 5, 2022

bniebuhr commented Mar 6, 2022

rsbivand commented Mar 6, 2022

bniebuhr commented Mar 7, 2022

`readRAST` and `writeRAST` support for `raster` and `terra` #42

`readRAST` and `writeRAST` support for `raster` and `terra` #42

rhijmans commented Feb 3, 2022 •

edited

rhijmans commented Feb 4, 2022 •

edited

rsbivand commented Feb 4, 2022 •

edited

rsbivand commented Feb 9, 2022 •

edited

rsbivand commented Feb 15, 2022 •

edited