r Terra issue with multicategorical raster. How to properly extract the categories and their values into layers without losing data?

Question

I am working with rTerra and having an issue with the CONUS historical disturbance dataset from LANDFIRE found here:https://landfire.gov/version_download.php (HDist is the name). To summarize what I want to do, I want to take this dataset, crop and project to my extent, then take the values of the cells and separate them as layers. So I want a layer for severity, one for disturbance type, etc. The historical disturbance data has these things all in one attribute table. In terra, this attribute table is set up under categories and this is providing a lot of problems. I have not had issues with the crop nor reproject, it is getting into the values and separating the categories into layers. I have the following code

library(terra)


setwd("your pathway to historical disturbance tif here")

h1 <- terra::rast("LC16_HDst_200.tif")  #read in the Hdist tif
h2 <- terra::project(h1, "EPSG:5070", method = "near") #project it using nearest neighbor
h3 <- crop(h2, ext([xmin,xmax,ymin,ymax]) #crop to the extent
h3

This then gives the output in the extent and projection I want but the main focus is the categories

categories  : Count, HDIST_ID, DISTCODE_V, DIST_TYPE, TYPE_CONFI, SEVERITY, SEV_CONFID, HDIST_CAT, FDIST, R, G, B

So I learned that with these kinds of datasets, the values are stored under these categories. if I plot with plot(h3) I only get the first row of the count category. In order to switch that category I can use

activeCat(h3) <- 4
h3

and I would get

name       :   DIST_TYPE
min value  :   Clearcut
max value  : Wildland Fire Use

The default active category was count, but now its DIST_TYPE, the fourth category, nothing too crazy. I try plotting

plot(h3)

I only get NoData plotted. None of the others. There is a function called catalyze() That claims to take your categories and converts them all into numerical layers

h4 <- catalyze(h3)

which gave me a thirteen layer dataset, which makes sense because there are 13 categories and it takes them and converts them into numeric layers. I tried plotting

plot(h4, 4) #plot h4 layer 4, which would correspond to DIST_TYPE category

it only plots a value of 8, and it looks to only show what is likely noData values. The map is mostly green, which is inline with the NoData from HDist. Anytime I try directly accessing values, it crashes. When I look at the min and max values I get 8 and 8 for min and max for that 'name" names: DIST_TYPE min values: 8 max values: 8. Other categories show a similar pattern. So it appeared to just take the first row of values for each category and make that the entire layer.

In summary, it is clear that terra stores all of the values that would easily be seen in an attribute table if the dataset were brought into arcgis. However, whenever I try to plot it or work with it, even before any real manipulation, it only accesses the top row of that attribute table, and when I catalyze, it just seems to mess everything up even more. I know this is really easy to solve in arcgis pro, but I want to keep everything in r from a documentation coherency standpoint. Any terra whizzes know what to do about this? I figure it has to be something very easy, but I don't really know what else to try. Maybe it is some major issue too. I have the same issue with LANDFIRE evt data. I have not had this issue with simple rasters such as dem, canopy cover, etc. It is only with these rasters with multiple categories (or columns in an attribute table)

edit this is the break image

Answer 1

That failed because the (ESRI) VAT IDs are not in the expected (for GDAL categories) 0..255 range. This has now been fixed and I get:

library(terra)
#terra version 1.4.6
r <- rast("LC16_HDst_200.tif")
activeCat(r) <- 4
r <- crop(r, ext(-93345, -57075, 1693125, 1716735))
plot(r)