R terra retain labels for categorical raster after writing

Question

My goal is to read in a large landcover dataset (this one: https://www.mrlc.gov/data), crop it, and write the smaller resulting raster to file. This dataset, when downloaded, appears as 3 files with types .ige, .img and .xml.

The problem I am trying to solve is that I can crop the raster just fine, it lists the categories, and returns is.factor(r) = TRUE. But if I write it to file as a .tif, and read it back in, the categorical labels are gone, and there are just integers instead.

Before writing the cropped raster to file:

> xxCrop
class       : SpatRaster 
dimensions  : 6658, 9840, 1  (nrow, ncol, nlyr)
resolution  : 30, 30  (x, y)
extent      : 1755285, 2050485, 2115315, 2315055  (xmin, xmax, ymin, ymax)
coord. ref. : Albers_Conical_Equal_Area 
source(s)   : memory
color table : 1 
varname     : nlcd_2021_land_cover_l48_20230630 
categories  : NLCD Land Cover Class, Histogram, Red, Green, Blue, Opacity 
name        :        NLCD Land Cover Class 
min value   :                 Unclassified 
max value   : Emergent Herbaceous Wetlands 

... and after:

is.factor(zz)
[1] FALSE
zz
class       : SpatRaster 
dimensions  : 6658, 9840, 1  (nrow, ncol, nlyr)
resolution  : 30, 30  (x, y)
extent      : 1755285, 2050485, 2115315, 2315055  (xmin, xmax, ymin, ymax)
coord. ref. : Albers_Conical_Equal_Area 
source      : nlcd_2021_landCover_LongIsland.tif 
name        : NLCD Land Cover Class 
min value   :                     0 
max value   :                    95 

I tested this out with a dummy dataset:

cls <- data.frame(id=1:3, cover=c("forest", "water", "urban"))
levels(r) <- cls
is.factor(r)
[1] TRUE
r
class       : SpatRaster 
dimensions  : 10, 10, 1  (nrow, ncol, nlyr)
resolution  : 36, 18  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 
source(s)   : memory
categories  : cover 
name        :  cover 
min value   : forest 
max value   :  urban 

writeRaster(r, 'catRaster.tif')

and when I read it back in:

rr <- rast('catRaster.tif')
rr
class       : SpatRaster 
dimensions  : 10, 10, 1  (nrow, ncol, nlyr)
resolution  : 36, 18  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
source      : catRaster.tif 
categories  : cover 
name        :  cover 
min value   : forest 
max value   :  urban 

So in this example, the categorical labels are retained.

What do I need to do to not lose the labels in the large dataset I am working with?

Answer 1

You are not showing the code you use, but the below works as expected for me:

f <- "nlcd_2021_land_cover_l48_20230630.img"
if (!file.exists(f)) {
    url <- "https://s3-us-west-2.amazonaws.com/mrlc/nlcd_2021_land_cover_l48_20230630.zip"
    download.file(url, basename(url), mode="wb")
    unzip(basename(url))
}

library(terra)
#terra 1.7.55

r <- rast(f)
e <- ext(1755285, 2050485, 2115315, 2315055)

x <- crop(r, e, filename="test.tif", overwrite=TRUE))
x
#class       : SpatRaster 
#dimensions  : 6658, 9840, 1  (nrow, ncol, nlyr)
#resolution  : 30, 30  (x, y)
#extent      : 1755285, 2050485, 2115315, 2315055  (xmin, xmax, ymin, ymax)
#coord. ref. : Albers_Conical_Equal_Area 
#source      : test.tif 
#color table : 1 
#varname     : nlcd_2021_land_cover_l48_20230630 
#categories  : NLCD Land Cover Class, Histogram, Red, Green, Blue, Opacity 
#name        :        NLCD Land Cover Class 
#min value   :                 Unclassified 
#max value   : Emergent Herbaceous Wetlands

Or in two steps as you seem to be doing

y <- crop(r, e)
y <- writeRaster(y, "test.tif", overwrite=TRUE)
y
#class       : SpatRaster 
#dimensions  : 6658, 9840, 1  (nrow, ncol, nlyr)
#resolution  : 30, 30  (x, y)
#extent      : 1755285, 2050485, 2115315, 2315055  (xmin, xmax, ymin, ymax)
#coord. ref. : Albers_Conical_Equal_Area 
#source      : test.tif 
#color table : 1 
#categories  : NLCD Land Cover Class, Histogram, Red, Green, Blue, Opacity 
#name        :        NLCD Land Cover Class 
#min value   :                 Unclassified 
#max value   : Emergent Herbaceous Wetlands 

Perhaps you are using an old version of "terra"?