3

I tried to read a DEM raster using getData (from the raster package), then convert the RasterLayer to a SpatRaster (terra package). First step worked, but the second one failed.

library(raster)
library(terra)
 
(alt <- getData('alt', country='DEU', mask=T))
#class      : RasterLayer 
#dimensions : 960, 1116, 1071360  (nrow, ncol, ncell)
#resolution : 0.008333333, 0.008333333  (x, y)
#extent     : 5.8, 15.1, 47.2, 55.2  (xmin, xmax, ymin, ymax)
#crs        : +proj=longlat +datum=WGS84 +no_defs 
#source     : D:/dummy/DEU_msk_alt.grd 
#names      : DEU_msk_alt 
#values     : -179, 2701  (min, max)

altT <- rast(alt)            
# rast is supposed to be able to read RasterLayer, but it fails. Why?
# Error : [rast] cannot read from D:/dummy/DEU_msk_alt.grd

Some hints? :

rast("DEU_msk_alt.grd")
# Error : [rast] cannot read from D:/dummy/DEU_msk_alt.grd

rast("DEU_msk_alt.vrt")
#class       : SpatRaster 
#dimensions  : 960, 1116, 1  (nrow, ncol, nlyr)
#resolution  : 0.008333333, 0.008333333  (x, y)
#extent      : 5.8, 15.1, 47.2, 55.2  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +ellps=WGS84 +no_defs 
#data source : DEU_msk_alt.vrt 
#names       : DEU_msk_alt

It looks like the rast function is looking for a .vrt file, while getData associated the raster to the grd file. Whatever, rast should work when applied to a RasterLayer, according to the documentation.

Any idea? How to convert such a RasterLayer object to a terra object? What do I miss? Thanks in advance,

JL

Robert Hijmans
  • 40,301
  • 4
  • 55
  • 63
Jean-Luc Dupouey
  • 341
  • 1
  • 6

1 Answers1

2

This occurs because the GDAL driver does not detect the correct datatype because there is a trailing space in the datatype description in the .grd file: "INT2S " instead of "INT2S"

The raster package uses its own code to read these files. terra has a stronger dependence on GDAL and uses it for reading all file types. As these are relatively small files, you can also transfer them like this:

library(terra)
alt <- raster::getData('alt', country='DEU', mask=TRUE)
x <- rast(alt * 1)

Or your work-around

y <- rast( gsub("grd$", "vrt", filename(alt)) )

The current development version of terra (version 1.1-18; April 2021) can now read these files, even though it still emits a warning

#Unhandled datatype=INT2S  (GDAL error 1)
Robert Hijmans
  • 40,301
  • 4
  • 55
  • 63
  • 1
    Thanks, it works and solved my problem. It could be nice to add something in the rast documentation, explaining that it does not work for all of Raster* objects. Currently, the documentation tells: "x: filename (character), SpatRaster ..., or Raster* object (from the "raster" package).". It is not exact and should perhaps be mentionned. – Jean-Luc Dupouey Dec 14 '20 at 23:05