TransWikia.com

Attribute table in R terra package?

Geographic Information Systems Asked on April 8, 2021

I am interested in the terra R package as it appears to be a successor to the raster package that is being more actively developed than the raster package. But it seems like some features are missing.

For example I have an aggregated National Landcover Dataset raster in .GRD format.

If I call x <- raster::raster('nlcd_agg.grd') the resulting object has a slot x@data@attributes with the raster attribute table that includes the names of the land cover types corresponding to the integer values in the raster.

However if I call x <- terra::rast('nlcd_agg.grd') I cannot find that attribute table in the object nor do I see how to associate the object x with the attribute table of the raster.

I would like to be able to load the raster and have the associated attribute table with land cover type names and default plotting colors load as well. I understand that .GRD is a native format for the raster package but is there some way to achieve this with the terra or stars packages?

Here are two small files that can replicate the behavior (nlcd_agg.grd and the associated nlcd_agg.gri):

One Answer

See the cats and rats methods. Or levels which is the same as cats. cats is for categories (ID and label), rats for Raster Attribute Tables (ID and n labels). cats are more limited than rats, but generally easier to deal with. They are used to have rasters that behave like "factors".

With the specific file, the problem is that terra uses GDAL to read it, and the GDAL driver for GRD/GRI is rather incomplete. RATS are read fine, I think, when using e.g. GeoTIFF. I may write a translator to make it easier to extract this from GRD/GRI. But for now, what you can do is something along these lines:

library(raster)
f <- "nlcd_agg.grd"
r <- raster(f)

library(terra)
x <- rast(r)

# get the attributes
lev <- levels(r)[[1]]

Now set the attributes

# this fails because of a bug 
rats(x) <- lev
# work around 
rat <- lev
rats(x) <- lev
str(rats(x))

They are there, but they are ignored. I would use cats/levels.

lev <- lev[, c("ID", "Land.Cover.Class")]
lev[,2] <- as.character(lev[,2])

x <- rast(r)
levels(x) <- lev
is.factor(x)
#[1] TRUE

x
#class       : SpatRaster 
#dimensions  : 101, 121, 1  (nrow, ncol, nlyr)
#resolution  : 3750, 3750  (x, y)
#extent      : 1394535, 1848285, 1722765, 2101515  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=aea +lat_0=23 +lon_0=-96 +lat_1=29.5 +lat_2=45.5 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
#source      : memory 
#name        : nlcd_2011_landcover_2011_edition_2014_03_31 
#min value   :                                Unclassified 
#max value   :                                             

# legend shows class names    
plot(x)

Labels are returned as cell values

x[c(7357, 5047, 7360, 9307)]
#  nlcd_2011_landcover_2011_edition_2014_03_31
#1                          Perennial Snow/Ice
#2                                  Open Water
#3                          Perennial Snow/Ice
#4                                Unclassified

And the categories are stored if you save the raster to a GeoTIFF

z <- writeRaster(x, "test.tif", overwrite=TRUE)
z

You can also set the color-table

lev <- levels(r)[[1]]
rgb <- lev[,c("Red", "Green", "Blue", "Opacity")]
coltab(x) <- round(rgb*255)
plot(x)

But I did not consider the case where there is a color table and attributes. To be done.

Answered by Robert Hijmans on April 8, 2021

Add your own answers!

Ask a Question

Get help from others!

© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP