6

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 this image shows what it's plotting even after the fixes right now

souma
  • 63
  • 4

1 Answers1

3

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)

enter image description here

Robert Hijmans
  • 40,301
  • 4
  • 55
  • 63
  • So I have copied your code and it works, however, for some reason when I do a crop to the western US it breaks, my bounds go a bit north and south of the US so I wonder if that is causing it. When I crop to my state, it works. When i plot the full data it works. I also notice that reprojecting is causing similar issues. Catalyze is doing the same problem of seemingly only taking the first row of what would be an attribute table. The plotting definitely works better but with catalyze not working well that makes things trickier for what I need to do. – souma Oct 01 '21 at 21:59
  • That is useful feedback. "it breaks" is not very descriptive, it would be very useful if you showed what you do; and what messages you get, perhaps in a new question. It is also much better to ask one question at a time. E.g. catalyze is a different matter. – Robert Hijmans Oct 01 '21 at 22:18
  • So when I do `r <- rast("LC16_HDst_200.tif")` `r <- crop(r, ext(-3213000,583010, 764000,3947040))` `activeCat(r) <- 4` `plot(r)` I get the image that was added in an edit to the main post. (i Couldn't quickly find how to put images in comments). I get the same type of image, but for the whole us if I reproject. I'll do catalyze in a separate post later – souma Oct 01 '21 at 22:51
  • 1
    OK, I see that now. This is because it is a large area and the values have to be written to a (temporary) file. I need to add the ability to write ESRI VATs. – Robert Hijmans Oct 01 '21 at 23:37
  • I figured it was it was a large area problem. I am glad that it's now clarified and I can move forward with the catalyze question, thanks so much for your help by the way! – souma Oct 02 '21 at 00:34
  • Feel free to open an issue on https://github.com/rspatial/terra/issues (for better support for ESRI VATs) – Robert Hijmans Oct 02 '21 at 16:52
  • I will do that just so that there is a record of it, let me know if I can improve it, I'm kind of new to collaborative coding and the rules of boards. – souma Oct 03 '21 at 14:52