support SpatRaster with color table or RGB

[rhijmans]

Pull Request

This PR extends #807 (support for SpatRaster objects from the terra package). With this one there is support for "RGB(A)" SpatRasters (in which layers represent color channels), and for rasters with a color table (that match cell values with colors).

library(terra)
library(leaflet)

# RGB
s <- rast(system.file("ex/logo.tif", package="terra"))  
crs(s) <- "+proj=longlat" 
leaflet() |> addTiles() |> addRasterImage(s)

# color table
r <- rast("/vsicurl/https://geodata.ucdavis.edu/test/pr_nlcd.tif")
r[550:750, ] <- NA
leaflet() |> addTiles() |> addRasterImage(r, opacity = 0.75)

# difference with a continuous raster
x <- r * 1
leaflet() |> addTiles() |> addRasterImage(x, opacity = 0.75)

PR task list:

• Update NEWS (already covered by prior PR)
• Add tests (where appropriate)
  • R code tests: tests/testthat/
• Update documentation with devtools::document()

[jcheng5]

Hi, sorry I haven't gotten to this yet, I've been on vacation and then busy with rstudio::conf preparations. It's on my todo list though, thanks for your patience.

[amfriesz]

Hey all, I'm curious what the status of this PR is? I have a tutorial that has transitioned from the raster package to the terra package and now I'm having problems visualizing with leaflet. I've tested the forked update from @rhijmans and it seems to work well.

[jcheng5]

@amfriesz It's always taunting me from the top of my todo list 😞 I'll take a look at it now

[jcheng5]

Looks good overall!

How would you add a legend for these? (For both discrete and continuous) Currently, the "happy path" for doing so, is to pass a palette function to addLegend(). I think one easy fix for this would be to have a specific helper function that creates a palette function from a color-carrying SpatRaster.

[rhijmans]

If a SpatRaster has colors, either via a "color table" or because it has RGB layers, but it does not have levels (it is not categorical), then it is a "photo", e.g. a satellite image or an air-photo. There should not be a legend for these. For example object s above.

If a SpatRaster has levels it needs a categorical legend. The colors should be set by the color table, if there is one (if it is not NULL), and otherwise with a categorical palette. The former applies to SpatRaster r.

library(terra)
r <- rast("/vsicurl/https://geodata.ucdavis.edu/test/pr_nlcd.tif")

coltab(r)[[1]] |> head()
#  value red green blue alpha
#1     0 255   255  255   255
#2     1   0     0    0   255
#3     2   0     0    0   255
#4     3   0     0    0   255
#5     4   0     0    0   255
#6     5   0     0    0   255

levels(r)[[1]] |> head()
#  value       NLCD Land Cover Class
#1     0                Unclassified
#2    11                  Open Water
#3    12          Perennial Snow/Ice
#4    21       Developed, Open Space
#5    22    Developed, Low Intensity
#6    23 Developed, Medium Intensity

The levels and colors can be matched by the first column (in this case they are both called "value" but that is not guaranteed).

For continous rasters I use leaflet::addLegend and that seems to work fine (see terra::plot_let which is a base::plot-like wrapper around leaflet); but I have not tried to make that work for categorical rasters yet.

[jcheng5]

@rhijmans Thank you, that's extremely helpful. Give me a day or two to come up with a proposal for integrating the color table either with the color palette functions, or addLegend directly.

[amfriesz]

Hey @jcheng5, what's the status of this issue?

[jcheng5]

@amfriesz Thanks for the nudge.

@rhijmans How does this look? It's only intended for SpatRaster objects with a color table and categorical values; levels are optional. In particular, am I right in using numeric indexing into the coltab(x)[[layer]] and levels(x)[[layer]] objects to find the values and labels?

addRasterLegend <- function(map, x, layer = 1, ...) {
  stopifnot(inherits(x, "SpatRaster"))
  stopifnot(length(layer) != 1)

  ct <- coltab(x)[[layer]]
  if (is.null(ct)) {
    stop("addRasterLegend() can only be used on layers with color tables (see ?terra:coltab). Otherwise, use addLegend().")
  }
  colors <- grDevices::rgb(ct$red/255, ct$green/255, ct$blue/255, ct$alpha/255)
  lvls <- levels(x)[[layer]]
  if (is.data.frame(lvls)) {
    labels <- lvls[[2]]
    colors <- colors[match(lvls$value, ct[[1]])]
  } else {
    labels <- ct[[1]]
  }
  addLegend(map, colors = colors, labels = labels, ...)
}

Example:

library(terra)
library(leaflet)

r <- rast("/vsicurl/https://geodata.ucdavis.edu/test/pr_nlcd.tif")
r[550:750, ] <- NA
leaflet() |>
  addTiles() |>
  addRasterImage(r, opacity = 0.75) |>
  addRasterLegend(r, opacity = 0.75)

[rhijmans]

Two suggestions

addRasterLegend <- function(map, x, layer = 1, ...) {
  stopifnot(inherits(x, "SpatRaster"))
##RH1 should be ==, not !=
  stopifnot(length(layer) == 1)
## could be expanded to
  # stopifnot((layer>0) && (layer <= nlyr(x)) && (length(layer) == 1))


  ct <- coltab(x)[[layer]]
  if (is.null(ct)) {
    stop("addRasterLegend() can only be used on layers with color tables (see ?terra:coltab). Otherwise, use addLegend().")
  }
  colors <- grDevices::rgb(ct$red/255, ct$green/255, ct$blue/255, ct$alpha/255)
  lvls <- levels(x)[[layer]]
  if (is.data.frame(lvls)) {
    labels <- lvls[[2]]
    colors <- colors[match(lvls$value, ct[[1]])]
  } else {
##RH2 the below creates a better legend for the example raster as the values that are
## not used are excluded (and there are many). It does not affect the example much.
    labels <- unique(x[[layer]])[,1]
    colors <- colors[match(labels, ct[[1]])]	

## the downside is that computing the unique values may take some time for a large raster.
## an approximate but perhaps good enough and always fast solution would be 
#    d <- duplicated(ct[,-1])
#    ct <- ct[!d, ]
#    labels <- ct[,1]
#    colors <- colors[match(labels, ct[[1]])]
  }
  addLegend(map, colors = colors, labels = labels, ...)
}

Examples

library(terra)
library(leaflet)

r <- rast("/vsicurl/https://geodata.ucdavis.edu/test/pr_nlcd.tif")
r[550:750, ] <- NA
leaflet() |>
  addTiles() |>
  addRasterImage(r, opacity = 0.75) |> 
  addRasterLegend(r, opacity = 0.75)

rr <- r
levels(rr)  <- NULL 
leaflet() |>
  addTiles() |>
  addRasterImage(rr, opacity = 0.75) |> 
  addRasterLegend(rr, opacity = 0.75)

[rhijmans]

"?terra:coltab" should be "?terra::coltab"

[jcheng5]

Thank you for that feedback. Is it also worthwhile to remove unused levels in the case where levels(r) is non-NULL?

[rhijmans]

Yes, that situation is not uncommon. The below is one way to handle that (with changed test data)

addRasterLegend <- function(map, x, layer = 1, ...) {
  stopifnot(inherits(x, "SpatRaster"))
##RH1 should be ==, not !=
  stopifnot(length(layer) == 1)
## could be expanded to
  # stopifnot((layer>0) && (layer <= nlyr(x)) && (length(layer) == 1))


  ct <- coltab(x)[[layer]]
  if (is.null(ct)) {
    stop("addRasterLegend() can only be used on layers with color tables (see ?terra::coltab). Otherwise, use addLegend().")
  }
  colors <- grDevices::rgb(ct$red/255, ct$green/255, ct$blue/255, ct$alpha/255)
  lvls <- levels(x)[[layer]]
  if (is.data.frame(lvls)) {
##3 remove empty levels
    d <- duplicated(lvls[[2]])
	lvls <- lvls[!d, ]
    labels <- lvls[[2]]
    colors <- colors[match(lvls$value, ct[[1]])]
  } else {
##RH2 the below creates a better legend for the example raster as the values that are
## not used are excluded (and there are many). It does not affect the example much.
    labels <- unique(x[[layer]])[,1]
    colors <- colors[match(labels, ct[[1]])]	

## the downside is that computing the unique values may take some time for a large raster.
## an approximate but perhaps good enough and always fast solution would be 
#    d <- duplicated(ct[,-1])
#    ct <- ct[!d, ]
#    labels <- ct[,1]
#    colors <- colors[match(labels, ct[[1]])]
  }
  addLegend(map, colors = colors, labels = labels, ...)
}



library(terra)
library(leaflet)

r <- rast("/vsicurl/https://geodata.ucdavis.edu/test/pr_nlcd.tif")
r[550:750, ] <- NA
levels(r) <- merge(data.frame(value=0:255), levels(r)[[1]], all.x=TRUE)

leaflet() |>
  addTiles() |>
  addRasterImage(r, opacity = 0.75) |> 
  addRasterLegend(r, opacity = 0.75)

rr <- r
levels(rr)  <- NULL 
leaflet() |>
  addTiles() |>
  addRasterImage(rr, opacity = 0.75) |> 
  addRasterLegend(rr, opacity = 0.75)

[jcheng5]

@rhijmans I've addressed that situation (but using a different method--sorry, I didn't see your response until I already had it working). If you could give it one more look, I'd appreciate it!

[rhijmans]

@jcheng5 your approach is more elegant.

In my two commits I suggest using "unique(x)" instead of "values(x)" to avoid memory issues.
And using RGB(x) instead of x@ptr$rgb as the latter now fails in terra-devel.

It looks like I have also added assure_crs_terra --- that must have been a while ago, and I am not sure if it is useful.