[GIS] reading rasters in r using terra package

projectrrasterterra

Using the R terra package for the first time after many years using the raster package. Stuck on the basics.

Trying to read DEM raster into r and reproject it. The raster is in ESRI grid format and can be downloaded here.

I previously asked this question on r-sig-geo.

And had working solution with raster

dem_raster<-raster("/Users/mypath/w001001.adf")

Attempting the same with terra::rast

dem_terra<-rast("/Users/mypath/w001001.adf")

And plotting the results side by side…

enter image description here

Printing the details for the two rasters

dem_raster

class      : RasterLayer 
dimensions : 8827, 10224, 90247248  (nrow, ncol, ncell)
resolution : 1000, 1000  (x, y)
extent     : -5762000, 4462000, -3920000, 4907000  (xmin, xmax, ymin, ymax)
crs        : +proj=laea +lat_0=-100 +lon_0=6370997 +x_0=45 +y_0=0 +datum=WGS84 +units=m +no_defs 
source     : w001001.adf 
names      : w001001 
values     : -76, 5930  (min, max)
attributes :
         ID COUNT
 from:  -76     2
  to : 5930     1

dem_terra

class       : SpatRaster 
dimensions  : 8827, 10224, 1  (nrow, ncol, nlyr)
resolution  : 1000, 1000  (x, y)
extent      : -5762000, 4462000, -3920000, 4907000  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=laea +lat_0=-100 +lon_0=6370997 +x_0=45 +y_0=0 +datum=WGS84 +units=m +no_defs 
source      : w001001.adf 
name        : VALUE 
min value   :   -76 
max value   :  5930

How do I get rast to read in the ESRI grid values correctly?

Best Answer

The raster package fails for me with that data:

> r = raster(f)
proj_create: Error -22: lat_0, lat_1 or lat_2 >= 90
proj_create: Error -22: lat_0, lat_1 or lat_2 >= 90
Error in .rasterObjectFromFile(x, band = band, objecttype = "RasterLayer",  : 
  Cannot create a RasterLayer object from this file.

Oh well. terra works:

> r = rast(f)
> r
class       : SpatRaster 
dimensions  : 8827, 10224, 1  (nrow, ncol, nlyr)
resolution  : 1000, 1000  (x, y)
extent      : -5762000, 4462000, -3920000, 4907000  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=laea +lat_0=-100 +lon_0=6370997 +x_0=45 +y_0=0 +datum=WGS84 +units=m +no_defs 
source      : w001001.adf 
name        : VALUE 
min value   :   -76 
max value   :  5930

The problem appears to be due to terra thinking this is a categorical raster, with levels. That's why the plot appears differently. Also:

> is.factor(r)
[1] TRUE

which shouldn't be true for a numeric raster:

> is.factor(rast(matrix(1:12,3,4)))
[1] FALSE

but is for one made of letters:

> is.factor(rast(matrix(letters[1:12],3,4)))
[1] TRUE

I can't find any documentation on how to make rast read in as numeric, or how to convert from one to the other (as.numeric doesn't work) but adding "0" seems to work:

 plot(r+0)

enter image description here

There may be documentation somewhere or this may be changed in a later release, I guess.

[ask your second Q in a separate post]

Related Question