Search code examples
rgeospatialterra

Terra extracting incorrect value from categorical raster

I am attempting to update old code I have where I previously used raster to terra. I encounter an issue extracting dominant vegetation type from a LANDFIRE raster layer. In ArcGIS I created a new column titled "CONDENSED" for my condensed vegetation types (i.e., 1="Aspen", 2="Other", etc.)and this is the active category I would like to extract from. I have 8 categories. For some reason, the extracted value is one category "higher" than it should be. For example, a point should have a value of 1 for aspen, but the extracted value is 2/Other. This is consistent across all veg types and doesn't matter if my points are sf or SpatVector. I am using terra version 1.4.11.

I have tried creating a reproducible example from scratch and it works perfectly fine when compared to values from raster::extract(), so I'm not sure if the issue has to do with how I'm specifying the active layer? I have uploaded a small sample of the raster and points I am working with here: https://github.com/Cara-Thompson/Elk-resource-selection

Anyone know what might be going on? Below is the description of my SpatRaster and my steps.

BONUS: I noticed terra::extract() works with sf objects in full form when a Raster* object is used. So I can get it to extract correctly not using SpatVect/SpatRast formats, but is this any faster than using raster::extract()?

# Read in raster and define crs
Veg <- terra::rast("./Vegetation/veg_cond.tif")
crs(Veg) <- "EPSG:4326"

#> Veg
# class       : SpatRaster 
# dimensions  : 9401, 21368, 1  (nrow, ncol, nlyr)
# resolution  : 0.0003398576, 0.0003398576  (x, y)
# extent      : -113.7326, -106.4705, 32.44176, 35.63676  (xmin, xmax, ymin, ymax)
# coord. ref. : lon/lat WGS 84 (EPSG:4326) 
# source      : veg_cond.tif 
# categories  : COUNT, CONDENSED 
# name        :    COUNT 
# min value   :   135890 
# max value   : 77959803 

# Make the condensed category the active category so count isn't extracted
activeCat(Veg) <- is.factor("CONDENSED")

# Extract
Values <- terra::extract(Veg, st_coordinates(EP.sf))

# Verify using raster::extract()
Veg2 <- raster(Veg)
Values$raster <- raster::extract(Veg2, EP.sf)
Solution

  • This is a bug related to the difference between ESRI VAT and GDAL categories (see this issue) that now seems fixed. I now get

    library(terra)
#terra version 1.4.12
elk <- vect("elk subset/elk subset.shp")
veg <- rast("veg_sample/veg sample.tif")
# note the way to use activeCat!
activeCat(veg) <- 2

veg
#class       : SpatRaster 
#dimensions  : 933, 1272, 1  (nrow, ncol, nlyr)
#resolution  : 0.0003398576, 0.0003398576  (x, y)
#extent      : -109.554, -109.1217, 33.8277, 34.14478  (xmin, xmax, ymin, ymax)
#coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#source      : veg sample.tif 
#categories  : COUNT, CONDENSED 
#name        : CONDENSED 
#min value   :     ASPEN 
#max value   : PONDEROSA 

plot(veg, mar=c(2,2,2,8))
points(elk)

    enter image description here

    head(extract(veg, elk))
#  ID      CONDENSED
#1  1          ASPEN
#2  2  OAK/SHRUBLAND
#3  3  OAK/SHRUBLAND
#4  4 PINYON-JUNIPER
#5  5 PINYON-JUNIPER
#6  6      PONDEROSA
 
e <- ext(-109.18688, -109.18483,   34.11579,   34.11798)
cveg <- crop(veg, e)
activeCat(cveg) <- 2
celk <- crop(elk, e)

plot(cveg, mar=c(2,2,2,8))
points(celk)

    enter image description here

    extract(veg, celk)
#  ID     CONDENSED
#1  1 OAK/SHRUBLAND

    You can install the development version with

    install.packages('terra', repos='https://rspatial.r-universe.dev')

    BONUS:

    extract is a generic function. This means that there is no difference between calling terra::extract(x) or raster::extract(x) because the version of the function you get depends on the class of x. So if you do terra::extract(x) and x is a RasterLayer, you get the function that is defined in the raster package. In this case it is therefore better just to use extract(x) --- as that would not be misleading.