How to style rasters added to a leaflet map with `addTileFolder`?

[jpd-defra]

Hello,

I've created a set of tiles with gdal2tiles.py I can get them to display on the leaflet map with addTileFolder, but from the docs can't work out how to style the raster(s). Which options can I use to provide a colour palette?

My code looks like this:

leaflet() %>%
  addTiles() %>%
  addTileFolder(
    "/tmp/30by30_tiles_WGS84",
    tms = TRUE,
    layerId = NULL,
    group = NULL,
    attribution = NULL,
    options = leaflet::tileOptions(),
    data = leaflet::getMapData(.)
  )

And my map looks like this - the raster extent is correct.

If there's a bit more detail to add to the docs, I would be happy to do this in a PR.

Any help much appreciated!

[tim-salabim]

Hi, its a long time since i have looked into addTileFolder. Any chance you can supply the data?

[jpd-defra]

Hi @tim-salabim , I can't share my actual data, so I've created some dummy data to illustrate the problem with the following code:

library(raster)
set.seed(42) 

xmin <- 400000  
xmax <- 600000
ymin <- 0
ymax <- 300000 
resolution <- 500

r <- raster(xmn = xmin, xmx = xmax, ymn = ymin, ymx = ymax, 
            resolution = resolution, crs = CRS("+init=epsg:27700"))
values(r) <- sample(c(0, 1), ncell(r), replace = TRUE)
r_reprojected <- projectRaster(r, crs = CRS("+init=epsg:4326"), method = "ngb")
writeRaster(r_reprojected, "/tmp/binary_raster_epsg4326.tif", format = "GTiff", overwrite = TRUE)

Then I fed this into gdal and used:

gdal_translate -of VRT -ot Byte -scale /tmp/binary_raster_epsg4326.tif temp.vrt
gdal2tiles.py --zoom=3-10 --processes=8 temp.vrt /tmp/dummy_tiles_WGS84

And the tiled data is here.

I then add this dummy data to the map with:

leaflet() %>%
  addTiles() %>%
  setView(-0.91431, 51.59701, zoom = 7) %>%
  addTileFolder(
    "/tmp/dummy_tiles_WGS84",
    tms = TRUE,
    layerId = NULL,
    group = NULL,
    attribution = NULL,
    options = leaflet::tileOptions(),
    data = leaflet::getMapData(.)
)

And the map looks like this:

So it is styling a bit differently to my previous image, but my question still remains around how to change this styling.

[tim-salabim]

Thanks, I'll have a look

[tim-salabim]

Hi @jpd-defra
sorry for the longer silence on this. Been busy with other leafem related issues...
Turns out, addTileFolder() is not able to pass along any styling as it simply uses leaflet::addTiles() to render stuff. This means you need to do the styling when/prior to creating the tiles. Here's one way that seems to do what you are after (at least if I understand you request correctly):

library(leaflet)
library(leafem)
library(raster)
library(colourvalues)

set.seed(42) 

xmin <- 400000  
xmax <- 600000
ymin <- 0
ymax <- 300000 
resolution <- 500

r <- raster(xmn = xmin, xmx = xmax, ymn = ymin, ymx = ymax, 
            resolution = resolution, crs = CRS("+init=epsg:27700"))
# values(r) <- sample(c(0, 1), ncell(r), replace = TRUE)
values(r) = rnorm(ncell(r), 50, 10)
r_reprojected <- projectRaster(r, crs = CRS("+init=epsg:4326"), method = "ngb")
writeRaster(r_reprojected, "/tmp/binary_raster_epsg4326.tif", format = "GTiff", overwrite = TRUE)

clrs = colourvalues::color_values(values(r))
clrs = t(col2rgb(clrs))

clr_mat = cbind(elevation_value = values(r), clrs)
write.table(clr_mat, file = "/tmp/color_lookup.txt", sep = " ", row.names = FALSE)

## now run these gdal commands 
# gdaldem color-relief /tmp/binary_raster_epsg4326.tif /tmp/color_lookup.txt /tmp/binary_raster_epsg4326_clrs.tif
# gdal_translate -of VRT -ot Byte -scale -a_nodata none /tmp/binary_raster_epsg4326_clrs.tif /tmp/temp.vrt
# gdal2tiles.py --zoom=3-10 --processes=8 /tmp/temp.vrt /tmp/dummy_tiles_WGS84

leaflet() %>%
  addTiles() %>%
  setView(-0.91431, 51.59701, zoom = 7) %>%
  addTileFolder(
    "/tmp/dummy_tiles_WGS84",
    tms = TRUE,
    layerId = NULL,
    group = NULL,
    attribution = NULL,
    options = leaflet::tileOptions(),
    data = leaflet::getMapData(.)
  )

If you want to be able to style stuff from the R call, probably best to have a look at addGeoTiff() or we could consider adding a raster method to the addPM* functions, as PMTiles, supports raster data, but I haven't looked into that yet.

addGeotiff() example:

leaflet() %>%
  addTiles() %>%
  setView(-0.91431, 51.59701, zoom = 7) %>%
  addGeotiff(
    file = "/tmp/binary_raster_epsg4326.tif"
    , colorOptions = colorOptions(
      palette = hcl.colors(256, palette = "viridis")
      , na.color = "#333333"
    )
  )

Let me know what you think and how I can help.

[jpd-defra]

Hey @tim-salabim, thanks very much for looking into this and getting back to me. I've got a few thoughts/followups in two areas.

addTileFolder

I was definitely drawn to addTileFolder for use when I have raster data that is a bit too big for a shiny app to handle with addGeotiff. Your code works well for me.

I've tried to adapt it for a binary raster containing only values of 1 (and no data) but not having any luck. It always renders as black and white, regardless of the colour matrix I'm providing to GDAL. Ideally I wanted the 1 values to be a colour, and for no data to be transparent (apologies if this is more of a GDAL issue). For example, I'm making a raster using raster, burning in a value of 1 using rasterize, resulting in a raster with these variables:

class      : RasterLayer 
dimensions : 624, 422, 263328  (nrow, ncol, ncell)
resolution : 0.00716, 0.00449  (x, y)
extent     : -2.0358, 0.98572, 49.8316, 52.63336  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
source     : memory
names      : layer 
values     : 1, 1  (min, max) 

The color lookup I'm providing to GDAL looks like this:

"elevation_value" "red" "green" "blue"
1 62 132 177

If helpful, once we have a working continuous and binary example, I could create a PR adding both examples the documentation. Let me know.

Possible pmtile development for rasters

Given the above is a bit of a long way to go about styling 'big' rasters, and especially as the GDAL commands are tailored to elevation data and not other types of data, looking at pmtiles for raster would be great if you have the appetite for it. I haven't used pmtiles much before, but would be happy to help where I can!