Merge pull request #808 from rhijmans/main

support SpatRaster with color table or RGB

Author: jcheng5

NAMESPACE

@@ -76,6 +76,7 @@ export(addPolylines)
 export(addPopups)
 export(addProviderTiles)
 export(addRasterImage)
+export(addRasterLegend)
 export(addRectangles)
 export(addScaleBar)
 export(addSimpleGraticule)

NEWS

@@ -206,7 +206,9 @@ epsg3857 <- "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y
 #' @param x a \code{\link[terra]{SpatRaster}} or a \code{RasterLayer} object--see \code{\link[raster]{raster}}
 #' @param colors the color palette (see \code{\link{colorNumeric}}) or function
 #'   to use to color the raster values (hint: if providing a function, set
-#'   \code{na.color} to \code{"#00000000"} to make \code{NA} areas transparent)
+#'   \code{na.color} to \code{"#00000000"} to make \code{NA} areas transparent).
+#'   The palette is ignored if \code{x} is a SpatRaster with a color table or if
+#'   it has RGB channels.
 #' @param opacity the base opacity of the raster, expressed from 0 to 1
 #' @param attribution the HTML string to show as the attribution for this layer
 #' @param layerId the layer id
@@ -225,15 +227,20 @@ epsg3857 <- "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y
 #'   (before base64 encoding); defaults to 4MB.
 #' @template data-getMapData
 #'
+#' @seealso \code{\link{addRasterLegend}} for an easy way to add a legend for a
+#'   SpatRaster with a color table.
+#'
 #' @examples
 #' \donttest{library(raster)
 #'
 #' r <- raster(xmn = -2.8, xmx = -2.79, ymn = 54.04, ymx = 54.05, nrows = 30, ncols = 30)
 #' values(r) <- matrix(1:900, nrow(r), ncol(r), byrow = TRUE)
 #' crs(r) <- CRS("+init=epsg:4326")
 #'
+#' pal <- colorNumeric("Spectral", domain = c(0, 1000))
 #' leaflet() %>% addTiles() %>%
-#'   addRasterImage(r, colors = "Spectral", opacity = 0.8)
+#'   addRasterImage(r, colors = pal, opacity = 0.8) %>%
+#'   addLegend(pal = pal, values = c(0, 1000))
 #' }
 #' @export
 addRasterImage <- function(
@@ -283,6 +290,93 @@ addRasterImage <- function(
 }
 
 
+#' Add a color legend for a SpatRaster to a map
+#'
+#' A function for adding a [legend][addLegend()] that is specifically designed
+#' for [terra::SpatRaster] objects, with categorical values, that carry their
+#' own [color table][terra::coltab()].
+#'
+#' @param map a map widget object
+#' @param x a [SpatRaster][terra::SpatRaster] object with a color table
+#' @param layer the layer of the raster to target
+#' @param ... additional arguments to pass through to [addLegend()]
+#' @seealso [addRasterImage()]
+#' @examples
+#'
+#' library(terra)
+#'
+#' r <- rast("/vsicurl/https://geodata.ucdavis.edu/test/pr_nlcd.tif")
+#' leaflet() %>%
+#'   addTiles() %>%
+#'   addRasterImage(r, opacity = 0.75) %>%
+#'   addRasterLegend(r, opacity = 0.75)
+#'
+#' plot.new() # pause in interactive mode
+#'
+#' rr <- r
+#' levels(rr)  <- NULL
+#' leaflet() %>%
+#'   addTiles() %>%
+#'   addRasterImage(rr, opacity = 0.75) %>%
+#'   addRasterLegend(rr, opacity = 0.75)
+#'
+#' @md
+#' @export
+addRasterLegend <- function(map, x, layer = 1, ...) {
+  stopifnot(inherits(x, "SpatRaster"))
+  stopifnot(length(layer) == 1 && layer > 0 && layer <= terra::nlyr(x))
+  
+## might as well do this here and only once. Subsetting would otherwise have been necessary in 
+##  color_info <- base::subset(color_info, value %in% terra::values(x))
+  x <- x[[layer]]
+  
+  # Retrieve the color table from the layer. If one doesn't exist, that means
+  # the raster was colored some other way, like using colorFactor or something,
+  # and the regular addLegend() is designed for those cases.
+  ct <- terra::coltab(x)[[1]]
+  if (is.null(ct)) {
+    stop("addRasterLegend() can only be used on layers with color tables (see ?terra::coltab). Otherwise, use addLegend().")
+  }
+
+  # Create a data frame that has value and color columns
+  # Extract the colors in #RRGGBBAA format
+  color_info <- data.frame(
+    value = ct[[1]],
+    color = grDevices::rgb(ct$red/255, ct$green/255, ct$blue/255, ct$alpha/255)
+  )
+
+  lvls <- terra::levels(x)[[1]]
+
+  # Drop values that aren't part of the layer
+## unlike "values",  "unique" is memory-safe; it does not load all values 
+## into memory if the raster is large. So instead of:
+
+#  color_info <- base::subset(color_info, value %in% terra::values(x))
+
+## remove the levels to get the raw cell values
+  levels(x) <- NULL  
+  color_info <- base::subset(color_info, value %in% terra::unique(x)[[1]])
+
+  res <- if (is.data.frame(lvls)) {
+    # Use the labels from levels(x), and look up the matching colors in the
+    # color table
+
+    # The levels data frame can have varying colnames, just normalize them
+    colnames(lvls) <- c("value", "label")
+    base::merge(color_info, lvls, by.x = "value", by.y = 1)
+  } else {
+    # No level labels provided; use the values as labels
+    cbind(color_info, label = color_info$value)
+  }
+
+  # At this point, res is a data frame with `value`, `color`, and `label` cols,
+  # and values/colors not present in the raster layer have been dropped
+
+  addLegend(map, colors = res[["color"]], labels = res[["label"]], ...)
+}
+
+
+
 addRasterImage_RasterLayer <- function(
   map,
   x,
@@ -370,50 +464,70 @@ addRasterImage_SpatRaster <- function(
   data = getMapData(map)
 ) {
 
-  if (terra::nlyr(x) > 1) {
+  # terra 1.5-50 has terra::has.RGB()
+  if (has.RGB(x)) {
+    # RGB(A) channels to color table
+    x <- terra::colorize(x, "col")
+  } else if (terra::nlyr(x) > 1) {
     x <- x[[1]]
     warning("using the first layer in 'x'", call. = FALSE)
   }
 
   raster_is_factor <- terra::is.factor(x)
+
+  # there 1.5-50 has terra::has.colors(x)
+  ctab <- terra::coltab(x)[[1]]
+  has_colors <- !is.null(ctab)
+
   method <- match.arg(method)
   if (method == "ngb") method = "near"
   if (method == "auto") {
-    if (raster_is_factor) {
+    if (raster_is_factor || has_colors) {
       method <- "near"
     } else {
       method <- "bilinear"
     }
   }
 
-  if (project) {
-    # if we should project the data
-    projected <- projectRasterForLeaflet(x, method)
-  } else {
-    # do not project data
-    projected <- x
-  }
-
   bounds <- terra::ext(
     terra::project(
       terra::project(
         terra::as.points(terra::ext(x), crs=terra::crs(x)),
         epsg3857),
       epsg4326)
   )
+## can't the above be simplified to this?
+#  bounds <- terra::ext(
+#    terra::project(
+#        terra::as.points(terra::ext(x), crs=terra::crs(x)),
+#        epsg4326)
+#  )
+
+  if (project) {
+    # if we should project the data
+    x <- projectRasterForLeaflet(x, method)
+    if (method=="bilinear") {
+      has_colors <- FALSE
+    }
+  }
 
   if (!is.function(colors)) {
-    if (method == "near") {
+    if (method == "near" || has_colors) {
       # 'factors'
-      colors <- colorFactor(colors, domain = NULL, na.color = "#00000000", alpha = TRUE)
+      domain <- NULL
+      if (has_colors) {
+        colors <- rgb(ctab[,2], ctab[,3], ctab[,4], ctab[,5], maxColorValue=255)
+        domain <- ctab[,1]
+      }
+      colors <- colorFactor(colors, domain = domain, na.color = "#00000000", alpha = TRUE)
     } else {
       # 'numeric'
       colors <- colorNumeric(colors, domain = NULL, na.color = "#00000000", alpha = TRUE)
     }
   }
 
-  tileData <- terra::values(projected) %>% as.vector() %>% colors() %>% col2rgb(alpha = TRUE) %>% as.raw()
-  dim(tileData) <- c(4, ncol(projected), nrow(projected))
+  tileData <- terra::values(x) %>% as.vector() %>% colors() %>% col2rgb(alpha = TRUE) %>% as.raw()
+  dim(tileData) <- c(4, ncol(x), nrow(x))
   pngData <- png::writePNG(tileData)
   if (length(pngData) > maxBytes) {
     stop(

R/layers.R

@@ -24,8 +24,6 @@ pointData.SpatVector <- function(obj) {
 polygonData.SpatVector <- function(obj) {
   check_crs_terra(obj)
 
-  # this is a bit convoluted. I will add a simpler
-  # and more efficient method to terra to replace the below
   xy = data.frame(terra::geom(obj))
   names(xy)[3:4] = c("lng", "lat")
   xy = split(xy[,2:5], xy[,1]) # polygons
@@ -40,6 +38,9 @@ polygonData.SpatVector <- function(obj) {
     })
   })
 
+  # with terra >= 1.5-50 you can do this instead
+  # xy = terra::geom(obj, list=TRUE, xnm="lng", ynm="lat")
+
   structure(
     xy,
     bbox = terra_bbox(obj)
@@ -49,6 +50,21 @@ polygonData.SpatVector <- function(obj) {
 
 
 # helpers -----------------------------------------------------------------
+assure_crs_terra <- function(x) {
+  prj <- crs(x, proj=TRUE)
+  if (is.lonlat(x, perhaps=TRUE, warn=FALSE)) {
+    if (!grepl("+datum=WGS84", prj, fixed = TRUE)) {
+	  x <- project(x, "+proj=longlat +datum=WGS84")
+    }
+	return(x)
+  }
+  # Don't have enough information to check
+  if (is.na(crs) || (crs=="")) {
+    warning("SpatVector layer is not long-lat data", call. = FALSE)	
+    return(x)
+  }
+  project(x, "+proj=longlat +datum=WGS84")
+}
 
 check_crs_terra <- function(x) {
   crs <- crs(x)

R/normalize-terra.R

@@ -37,4 +37,11 @@ test_that("rasters", {
 
   expect_maps_equal(r1, r2)
 
+ # test with color map
+  r <- rast(ncols=10, nrows=10, vals=rep_len(10:15, length.out=100), xmin=0, xmax=10^6, ymin=0, ymax=10^6, crs=pmerc)
+  r[5,] <- NA
+  coltab(r) <- c(rep("#FFFFFF", 10), rainbow(6, end=.9))
+  (r3 <- rtest(r))
+
 })
+