reading rasters in r using terra package

Question

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...

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?

Answer 1

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)

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

[ask your second Q in a separate post]

Answer 2

Both raster and terra show that this is a categorical raster because it has a Raster Attribute Table (RAT). But these are different for AI GRID than for the GTiff that I have tested and developed with. The AI GRID attributes are not correctly interpreted and they are not that useful. All the RAT adds is the counts of the unique values. If I recall well, that is the case for all integer AI GRID files (my memory is blurry, I rarely encounter one of them these days). In the current development version of terra the RAT is ignored in cases like this and the SpatRaster is not categorical.

@Spacedman trick to add zero works fine, because the levels are dropped when doing arithmetic.

You can install the development version with

install.packages('terra', repos='https://rspatial.r-universe.dev')

plot(x) still looks a bit odd, because taking a sample with GDAL seems to mess up the NODATA value for this file. It appears to work well for all other purposes.