0

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?

1 Answer 1

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"?

1
  • 1
    I realize now that the problem was that I was ignoring the .tif.aux.xml file that was also being written. I had moved the cropped raster tif file and ignored the xml file, and apparently that .xml file is important in this context! thank you for your help.
    – Pascal
    Oct 26, 2023 at 18:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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