1

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?

1
  • One question per post please, then we have a nice set of single questions with best answers! I've answered your first, if you can cut your second Q out of here and paste it into a new post we can look at that.
    – Spacedman
    May 28, 2021 at 21:04

2 Answers 2

2

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]

4
  • Thanks @Spacedman. Can you explain why adding "0" works? I think ultimately I want the raster to read in correctly to be sure no errors when we go to actually work with the values in the raster... May 28, 2021 at 22:21
  • I can't definitively explain why it works, but I can imagine that at some point terra is converting the factor back to numeric to do the addition. The only docs on factors in terra I can find is in help(terra::levels) and there's not much in the code either. I strongly suggest you file an issue in the terra bug tracker, explain the problem and link to this Q. The terra dev does use SO so he may see this anyway.
    – Spacedman
    May 29, 2021 at 6:36
  • terra is still in development so you have to expect a bit of weirdness and minor hacks to fix things. Adding 0 seems minor enough to not break stuff, except it might result in a large on-disk raster being read into RAM, I guess. Finding the right way to do it, if that exists, would be better, but I don't see any documentation, and other guesses, like as.numeric(r) don't work.
    – Spacedman
    May 29, 2021 at 6:38
  • 2
    as.numeric was a good guess --- but I only added that a week ago! May 31, 2021 at 5:01
1

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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.