Match the raster values from a dbf table and filter based on raster

Question

I have a large raster covering Australia which I applied crop() and mask() to using a multi polygon sf (nested catchments). I read the dbf table using the foreign package, terra::rast() for the tif. I then ensured each file has the same CRS.

When I tried to match the dbf with cells in my raster, it doesn't give me what I expect because my raster doesn't have value. It just contain "forest categorize (native and non forest). I need to link data from the dbf table with my raster so I can filter cells.

For example, I want to select just the cells which are "Native forest" from the Forest_CATEGO column, and the "HT_CODE" and "COVER_CODE" columns %in% c(1,2,3). I then want to filter only those catchments which contain more than 60% of this condition in their area (here I should consider the catchments are nested). I don't know if this kind of filtering is possible?

This a view of nested catchments.

Here is my code:

library(terra)
#> terra 1.7.78

library(foreign)
library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:terra':
#> 
#>     intersect, union
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union

install.packages("foreign")
#> Warning: package 'foreign' is in use and will not be installed

# File paths
raster_file <- "C:/Users/for23.tif"
dbf_file <- "C:/Users/for23.tif.vat.dbf"
shapefile <- "C:/Users/boundry.shp"

# Load the raster, .dbf table, and shapefile
rast <- rast(raster_file)
rast
#> class       : SpatRaster 
#> dimensions  : 38370, 40100, 1  (nrow, ncol, nlyr)
#> resolution  : 100, 100  (x, y)
#> extent      : -1888000, 2122000, -4847000, -1010000  (xmin, xmax, ymin, ymax)
#> coord. ref. : GDA_1994_Albers 
#> source      : aus_for23.tif 
#> categories  : FOR_CATEGO, FOR_TYPE, FOR_CODE, COVER_CODE, HT_CODE, FOREST, STATE, FOR_SOURCE 
#> name        :    FOR_CATEGO 
#> min value   :    Non-forest 
#> max value   : Native forest

level <- levels(rast)
head(level)
#> [[1]]
#>     VALUE            FOR_CATEGO
#> 1       1            Non-forest
#> 2       2            Non-forest
#> 3       3         Native forest
#> 4       4         Native forest
#> 5       5         Native forest
#> 6       6         Native forest
#> 7       7         Native forest
#> 8       8         Native forest
#> 9       9         Native forest
#> 10     10         Native forest

dbf_table <- read.dbf(dbf_file)
head(dbf_table)
#>   VALUE     COUNT    FOR_CATEGO   FOR_TYPE FOR_CODE COVER_CODE HT_CODE FOREST
#> 1     1 121022073    Non-forest Non-forest        0          6       6      0
#> 2     2      5875    Non-forest Non-forest        0          6       6      0
#> 3     3     17703 Native forest   Mangrove       63          3       2      1
#> 4     4       335 Native forest   Mangrove       61          1       2      1
#> 5     5     16724 Native forest Rainforest       98          3       2      1
#> 6     6      4555 Native forest Rainforest       96          3       1      1
#>   STATE       FOR_SOURCE
#> 1   Qld             <NA>
#> 2  <NA>             <NA>
#> 3   Qld Global Mangroves
#> 4   Qld Global Mangroves
#> 5   Qld         NVIS 6.0
#> 6   Qld         NVIS 6.0

catchments <- vect(shapefile)
shapefile1<- st_read(shapefile)
#> Error in st_read(shapefile): could not find function "st_read"

if (st_crs(shapefile1) != crs(rast)) {
  shapefile_data <- st_transform(shapefile1, st_crs(rast))
}
#> Error in st_crs(shapefile1): could not find function "st_crs"

# Inspect the .dbf table
head(dbf_table)
#>   VALUE     COUNT    FOR_CATEGO   FOR_TYPE FOR_CODE COVER_CODE HT_CODE FOREST
#> 1     1 121022073    Non-forest Non-forest        0          6       6      0
#> 2     2      5875    Non-forest Non-forest        0          6       6      0
#> 3     3     17703 Native forest   Mangrove       63          3       2      1
#> 4     4       335 Native forest   Mangrove       61          1       2      1
#> 5     5     16724 Native forest Rainforest       98          3       2      1
#> 6     6      4555 Native forest Rainforest       96          3       1      1
#>   STATE       FOR_SOURCE
#> 1   Qld             <NA>
#> 2  <NA>             <NA>
#> 3   Qld Global Mangroves
#> 4   Qld Global Mangroves
#> 5   Qld         NVIS 6.0
#> 6   Qld         NVIS 6.0

croped<- crop(rast, vect(shapefile_data))
#> Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'vect': object 'shapefile_data' not found

masked_hrs<- mask(rast, vect(shapefile_data))
#> Error in h(simpleError(msg, call)): error in evaluating the argument 'mask' in selecting a method for function 'mask': error in evaluating the argument 'x' in selecting a method for function 'vect': object 'shapefile_data' not found

levels(masked_hrs)<- list(dbf_table)
#> Error: object 'masked_hrs' not found

# Define conditions
condition <- dbf_table$FOR_CATEGO == "Native forest" &
  dbf_table$COVER_CODE %in% c(1, 2, 3) &
  dbf_table$HT_CODE %in% c(1, 2, 3)

link to my data

Answer 1

Your SpatRaster has one layer, but you added multiple levels using list(dbf_table). Based on your description, it is probably better to simply assign levels using a two-column dataframe.

Here is a reprex that assumes:

• dbf_table is a summary of every combination of values in your SpatRaster e.g. 10 distinct categories in both your SpatRaster and dbf_table (or however many exist in your actual data)
• VALUE in your SpatRaster == VALUE in dbf_table

As head() only returned the first six rows of dbf_table, the solution below will only use the values as shown in your question (n = 6).

First, load required packages and create example data:

library(terra)
library(sf)
library(ggplot2)
library(tidyterra)

# Create example dbf file
dbf_table <- structure(list(VALUE = 1:6, COUNT = c(121022073L, 5875L, 17703L, 
335L, 16724L, 4555L), FOR_CATEGO = c("Non-forest", "Non-forest", 
"Native forest", "Native forest", "Native forest", "Native forest"
), FOR_TYPE = c("Non-forest", "Non-forest", "Mangrove", "Mangrove", 
"Rainforest", "Rainforest"), FOR_CODE = c(0L, 0L, 63L, 61L, 98L, 
96L), COVER_CODE = c(6L, 6L, 3L, 1L, 3L, 3L), HT_CODE = c(6L, 
6L, 2L, 2L, 2L, 1L), FOREST = c(0L, 0L, 1L, 1L, 1L, 1L), STATE = c("Qld", 
"<NA>", "Qld", "Qld", "Qld", "Qld"), FOR_SOURCE = c("<NA>", "<NA>", 
"Global Mangroves", "Global Mangroves", "NVIS 6.0", "NVIS 6.0"
)), row.names = c("1", "2", "3", "4", "5", "6"), class = "data.frame")

# Create SpatRaster based on your metadata (lower resolution example)
r <- rast(nrows = 3837, ncols = 4010,
          ext = ext(-1888000, 2122000, -4847000, -1010000),
          crs = "EPSG:3577")

# Assign values to SpatRaster in order to replicate your actual data
set.seed(123)
values(r) <- sample(dbf_table$VALUE,
                    size = ncell(r),
                    replace = TRUE,
                    prob = 1 / sum(dbf_table$COUNT) * dbf_table$COUNT)

# Set the levels for the SpatRaster
levels(r) <- dbf_table[, c("VALUE", "FOR_CATEGO")]

head(levels(r))
# [[1]]
#   VALUE    FOR_CATEGO
# 1     1    Non-forest
# 2     2    Non-forest
# 3     3 Native forest
# 4     4 Native forest
# 5     5 Native forest
# 6     6 Native forest

# Create example mask sf
sf_poly <- st_as_sfc(st_bbox(c(xmin = -1.5e6,
                               ymin = -3.5e6,
                               xmax = -1.44e6,
                               ymax = -3.44e6))) %>%
  st_as_sf(crs = 3577)

# Mask and crop r
r1 <- crop(r, vect(sf_poly))
r1 <- mask(r1, vect(sf_poly))

Given that you want to subset your data (based on your defined conditions), it makes sense to convert all other values to NA. If you need to retain all values/levels, skip this next step (note that if you do, the ggplot() examples below may also need to be modified).

# Keep only wanted values, change others to NA
r1[!(r1 %in% dbf_table[dbf_table$FOR_CATEGO == "Native forest" &
                         dbf_table$COVER_CODE %in% 1:3 &
                         dbf_table$HT_CODE %in% 1:3, 1])] <- NA

# Change levels to desired variable e.g. FOR_TYPE
levels(r1) <- dbf_table[, c("VALUE", "FOR_TYPE")]

r1
# class       : SpatRaster 
# dimensions  : 60, 60, 1  (nrow, ncol, nlyr)
# resolution  : 1000, 1000  (x, y)
# extent      : -1500000, -1440000, -3500000, -3440000  (xmin, xmax, ymin, ymax)
# coord. ref. : GDA94 / Australian Albers (EPSG:3577) 
# source(s)   : memory
# categories  : FOR_TYPE 
# name        :   FOR_TYPE 
# min value   :   Mangrove 
# max value   : Rainforest 

ggplot() +
  geom_spatraster(data = r1, maxcell = Inf) +
  scale_fill_discrete(na.translate = FALSE,
                    name = "Type") +
  coord_sf(expand = FALSE) +
  theme(panel.background = element_blank())

If you need more nuance in your output, you can do do something like:

# Change levels to desired variable e.g. FOR_CODE
levels(r1) <- dbf_table[, c("VALUE", "FOR_CODE")]

r1
# class       : SpatRaster 
# dimensions  : 60, 60, 1  (nrow, ncol, nlyr)
# resolution  : 1000, 1000  (x, y)
# extent      : -1500000, -1440000, -3500000, -3440000  (xmin, xmax, ymin, ymax)
# coord. ref. : GDA94 / Australian Albers (EPSG:3577) 
# source(s)   : memory
# categories  : FOR_CODE 
# name        : FOR_CODE 
# min value   :       63 
# max value   :       96 

# For legend (illustrative purposes only)
legend <- subset(dbf_table, dbf_table$VALUE %in% as.numeric(unique(values(r1))))

ggplot() +
  geom_spatraster(data = r1, maxcell = Inf) +
  scale_fill_manual(name = "Type/Code",
                    labels = paste(legend$FOR_TYPE, legend$FOR_CODE),
                    values = c("#61D04F" ,"#009292", "darkgreen"),
                    na.translate = FALSE) +
  coord_sf(expand = FALSE) +
  theme(panel.background = element_blank())

Note that even though your defined condition returns four rows from the dbf_table used in this example, only three distinct levels were within the extent of sf_polygon: