new addGeoRaster method

[tim-salabim]

This is a continuation of parts of r-spatial/mapview#305

The intent is to have a new method to use https://github.com/GeoTIFF/georaster-layer-for-leaflet to enable rendering of largish raster data.

This is currently limited to files, i.e. no in-memory objects are supported so far.

This is what we currently have

library(leaflet)
library(leafem)
library(stars)

# tst = read_stars("/media/timpanse/d8346522-ef28-4d63-9bf3-19fec6e13aab/bu_lenovo/software/data/global_elevation_0.tif")
tst = read_stars("/home/timpanse/software/data/srtm_39_03.tif")
st_is_longlat(tst)

if (!st_is_longlat(tst)) {
  tst = st_warp(tst, crs = 4326)
}

tst_ds = stars:::st_downsample(tst, n = 1)
dim(tst_ds)

fl = tempfile(fileext = ".tif")
write_stars(tst_ds, dsn = fl)

options(viewer = NULL)

leaflet() %>%
  addProviderTiles("Esri.WorldImagery") %>%
  leafem:::addGeotiff(
    file = fl
    , group = "test"
    , layerId = "testid"
    , resolution = 96
    , opacity = 1
    , colorOptions = leafem:::colorOptions(
      palette = viridisLite::inferno
      , breaks = seq(
        min(tst[[1]], na.rm = TRUE)
        , max(tst[[1]], na.rm = TRUE)
        , 500
      )
    )
  )

which yields:

and zoomed in:

[tim-salabim]

We can now also pass na.color to colorOptions

[tim-salabim]

We can now also pass custom JS rendering functions :-)

library(leaflet)
library(leafem)
library(stars)

# tst = read_stars("/media/timpanse/d8346522-ef28-4d63-9bf3-19fec6e13aab/bu_lenovo/software/data/global_elevation_0.tif")
# tst = read_stars("/home/timpanse/software/data/srtm_39_03.tif")
chrpsfl = "https://data.chc.ucsb.edu/products/CHIRPS-2.0/africa_monthly/tifs/chirps-v2.0.2020.05.tif.gz"
# chrpsfl = "https://data.chc.ucsb.edu/products/CHIRPS-2.0/global_annual/tifs/chirps-v2.0.2019.tif.gz"
dsn = file.path(tempdir(), basename(chrpsfl))

download.file(chrpsfl, dsn)

tiffl = gsub(".gz", "", dsn)
R.utils::gunzip(dsn, tiffl)

tst = read_stars(tiffl)
tst[[1]][tst[[1]] < 0] = NA

st_is_longlat(tst)

if (!st_is_longlat(tst)) {
  tst = st_warp(tst, crs = 4326)
}

tst_ds = stars:::st_downsample(tst, n = 1)
dim(tst_ds)

fl = tempfile(fileext = ".tif")
write_stars(tst_ds, dsn = fl)

pal = hcl.colors(256, "inferno")
brks = NULL

myCustomJSFunc = htmlwidgets::JS(
  "
    pixelValuesToColorFn = (raster, colorOptions) => {
      const cols = colorOptions.palette;
      var scale = chroma.scale(cols);

      if (colorOptions.breaks !== null) {
        scale = scale.classes(colorOptions.breaks);
      }
      var pixelFunc = values => {
        let clr = scale.domain([raster.mins, raster.maxs]);
        if (isNaN(values)) return colorOptions.naColor;
        if (values < 50) return chroma('pink').alpha(0.7).hex();
        return clr(values).hex();
      };
      return pixelFunc;
    };
  "
)

leaflet(options = leafletOptions(noWrap = TRUE)) %>%
  addProviderTiles("Esri.WorldImagery") %>%
  leafem:::addGeotiff(
    file = fl
    , group = "test"
    , layerId = "testid"
    , resolution = 96
    , opacity = 1
    , colorOptions = leafem:::colorOptions(
      palette = pal #viridisLite::cividis
      , breaks = brks
      , na.color = "transparent"
    )
    , pixelValuesToColorFn = myCustomJSFunc
  ) %>%
  addLayersControl(overlayGroups = "test")

[tim-salabim]

Proposal:

• rename the current method for files/urls to addGeotif (try to avoid file copy)
• use addGeoRaster for stars/raster/terra objects

Also:

• have two custome JS functions -> one for accessing the georaster object, one for pixel values

[tim-salabim]

We now have addGeotiff for files & addGeoRaster for stars, stars_proxy & raster objects. terra still to be figured out... Though, I guess that the fact that stars is much more flexible than raster/terra, i.e. the data model represented by raster/terra is a subset of the data model that stars can handle, we should focus on stars support and use st_as_stars for raster/terra objects. Is this a fair assumption @edzer?

[tim-salabim]

library(leaflet)
library(leafem)
library(stars)

chrpsfl = "https://data.chc.ucsb.edu/products/CHIRPS-2.0/africa_monthly/tifs/chirps-v2.0.2020.05.tif.gz"
dsn = file.path(tempdir(), basename(chrpsfl))

download.file(chrpsfl, dsn)

tiffl = gsub(".gz", "", dsn)
R.utils::gunzip(dsn, tiffl)

tst = read_stars(tiffl, proxy = TRUE)

pal = hcl.colors(256, "viridis")
brks = seq(0, 500, 50)

myCustomJSFunc = htmlwidgets::JS(
  "
    pixelValuesToColorFn = (raster, colorOptions) => {
      const cols = colorOptions.palette;
      var scale = chroma.scale(cols);

      if (colorOptions.breaks !== null) {
        scale = scale.classes(colorOptions.breaks);
      }
      var pixelFunc = values => {
        let clr = scale.domain([raster.mins, raster.maxs]);
        if (isNaN(values)) return colorOptions.naColor;
        if (values < 120) return colorOptions.naColor;
        return clr(values).hex();
      };
      return pixelFunc;
    };
  "
)

leaflet(options = leafletOptions(noWrap = TRUE)) %>%
  addProviderTiles("CartoDB.DarkMatter") %>%
  leafem:::addGeoRaster(
    x = tst
    , group = "test"
    , layerId = "testid"
    , resolution = 96
    , opacity = 1
    , colorOptions = leafem:::colorOptions(
      palette = pal #viridisLite::cividis
      , breaks = brks
      , na.color = "transparent"
    )
    , pixelValuesToColorFn = myCustomJSFunc
  ) %>%
  addLayersControl(overlayGroups = "test")

[edzer]

Though, I guess that the fact that stars is much more flexible than raster/terra, i.e. the data model represented by raster/terra is a subset of the data model that stars can handle, we should focus on stars support and use st_as_stars for raster/terra objects. Is this a fair assumption @edzer?

That would be great, as it would also concentrate "external" knowledge/understanding of the unwritten raster/terra API in one place.

[tim-salabim]

That would be great, as it would also concentrate "external" knowledge/understanding of the unwritten raster/terra API in one place.

@edzer I don't quite understand what you mean here. Do you mean that we avoid code duplication?

[edzer]

Yes.

[tim-salabim]

Here's one issue I am not sure how to solve.

In order to make addGeotiff/addGeoRaster work with leaflet.extras2::addSidebyside, we need to place 2 rasters on different panes. Creating and placing 1 raster on a custom pane works fine, but for 2 calls, both rasters get placed on the pane of the last call. Reprex below.

library(leaflet)
library(leafem)
library(stars)

chrpsfl_04 = "https://data.chc.ucsb.edu/products/CHIRPS-2.0/africa_monthly/tifs/chirps-v2.0.2020.04.tif.gz"
chrpsfl_05 = "https://data.chc.ucsb.edu/products/CHIRPS-2.0/africa_monthly/tifs/chirps-v2.0.2020.05.tif.gz"

dsn_04 = file.path(tempdir(), basename(chrpsfl_04))
dsn_05 = file.path(tempdir(), basename(chrpsfl_05))

download.file(chrpsfl_04, dsn_04)
download.file(chrpsfl_05, dsn_05)

tiffl_04 = gsub(".gz", "", dsn_04)
R.utils::gunzip(dsn_04, tiffl_04)

tiffl_05 = gsub(".gz", "", dsn_05)
R.utils::gunzip(dsn_05, tiffl_05)

## add 1 layer to custom pane - works, layer is correctly rendered on 'leaflet-left-pane'
leaflet() %>%
  addTiles() %>%
  addMapPane("left", 200) %>%
  leafem:::addGeotiff(
    file = tiffl_04
    , group = "april"
    , layerId = "april_id"
    , resolution = 96
    , opacity = 1
    , options = tileOptions(
      pane = "left"
    )
    , colorOptions = leafem:::colorOptions(
      palette = hcl.colors(256, "inferno")
      , breaks = seq(0, 700, 10)
      , na.color = "transparent"
    )
  ) %>%
  addLayersControl(overlayGroups = "april")

## add 2 layers to 2 custom panes - doesn't work, both rendered on 'leaflet-right-pane'
leaflet() %>%
  addTiles() %>%
  addMapPane("left", 200) %>%
  addMapPane("right", 201) %>%
  leafem:::addGeotiff(
    file = tiffl_04
    , group = "april"
    , layerId = "april_id"
    , resolution = 96
    , opacity = 1
    , options = tileOptions(
      pane = "left"
    )
    , colorOptions = leafem:::colorOptions(
      palette = hcl.colors(256, "inferno")
      , breaks = seq(0, 700, 10)
      , na.color = "transparent"
    )
  ) %>%
  leafem:::addGeotiff(
    file = tiffl_05
    , group = "may"
    , layerId = "may_id"
    , resolution = 96
    , opacity = 1
    , options = tileOptions(
      pane = "right"
    )
    , colorOptions = leafem:::colorOptions(
      palette = hcl.colors(256, "viridis")
      , breaks = seq(0, 700, 10)
      , na.color = "transparent"
    )
  ) %>%
  addLayersControl(overlayGroups = c("april", "may"))

My suspicion is that this is a result of the JavaScript fetch call that we use to read the data. This creates a promise that I think is not evaluated until the whole map is rendered at which point the current pane argument is the one of the last addGeotiff call.

One thing that baffles me though is that for addFgb we also use fetch and when adding more than one layer, they do get placed on the correct pane... https://github.com/r-spatial/leafem/blob/master/inst/htmlwidgets/lib/FlatGeoBuf/fgb.js#L39-L104

Ping @timelyportfolio

[timelyportfolio]

@tim-salabim the joys of scope... I think you can fix by changing lines to

  var pane;  // could also use let
  if (options.pane === undefined) {
    pane = 'tilePane';
  } else {
    pane = options.pane;
  }

If I understand correctly, this is the correct result. Let me know if I missed something.

[tim-salabim]

Thanks so much @timelyportfolio !!

Turns out I don't understand scope as much as I though, eh... I must admit I feel stupid now, but at the same time I think this issue just hammered home the point of initialising variables in JS for me (at least I hope so).

This now works nicely:

library(leaflet)
library(leafem)
library(stars)

chrpsfl_04 = "https://data.chc.ucsb.edu/products/CHIRPS-2.0/africa_monthly/tifs/chirps-v2.0.2020.04.tif.gz"
chrpsfl_05 = "https://data.chc.ucsb.edu/products/CHIRPS-2.0/africa_monthly/tifs/chirps-v2.0.2020.05.tif.gz"

dsn_04 = file.path(tempdir(), basename(chrpsfl_04))
dsn_05 = file.path(tempdir(), basename(chrpsfl_05))

download.file(chrpsfl_04, dsn_04)
download.file(chrpsfl_05, dsn_05)

tiffl_04 = gsub(".gz", "", dsn_04)
R.utils::gunzip(dsn_04, tiffl_04)

tiffl_05 = gsub(".gz", "", dsn_05)
R.utils::gunzip(dsn_05, tiffl_05)

pal = hcl.colors(256, "inferno")
brks = seq(0, 1000, 10)

myCustomJSFunc = htmlwidgets::JS(
  "
    pixelValuesToColorFn = (raster, colorOptions) => {
      const cols = colorOptions.palette;
      var scale = chroma.scale(cols);

      if (colorOptions.breaks !== null) {
        scale = scale.classes(colorOptions.breaks);
      }
      var pixelFunc = values => {
        let clr = scale.domain([raster.mins, raster.maxs]);
        if (isNaN(values)) return colorOptions.naColor;
        if (values < 120) return colorOptions.naColor;
        return clr(values).hex();
      };
      return pixelFunc;
    };
  "
)

pal = hcl.colors(256, "inferno")
brks = seq(0, 1000, 10)

## add 2 layers to 2 custom panes - doesn't work, both rendered on pane from last call
leaflet() %>%
  addTiles() %>%
  addMapPane("left", 200) %>%
  addMapPane("right", 201) %>%
  addProviderTiles(
    "CartoDB.DarkMatter"
    , group = "carto_left"
    , options = tileOptions(pane = "left")
    , layerId = "leftid"
  ) %>%
  addProviderTiles(
    "CartoDB.DarkMatter"
    , group = "carto_right"
    , options = tileOptions(pane = "right")
    , layerId = "rightid"
  ) %>%
  leafem:::addGeotiff(
    file = tiffl_04
    , group = "april"
    , layerId = "april_id"
    , resolution = 96
    , opacity = 1
    , options = tileOptions(
      pane = "left"
    )
    , colorOptions = leafem:::colorOptions(
      palette = pal
      , breaks = brks
      , na.color = "transparent"
    )
    , pixelValuesToColorFn = myCustomJSFunc
  ) %>%
  leafem:::addGeotiff(
    file = tiffl_05
    , group = "may"
    , layerId = "may_id"
    , resolution = 96
    , opacity = 1
    , options = tileOptions(
      pane = "right"
    )
    , colorOptions = leafem:::colorOptions(
      palette = pal
      , breaks = brks
      , na.color = "transparent"
    )
    , pixelValuesToColorFn = myCustomJSFunc
  ) %>%
  leaflet.extras2::addSidebyside(
    layerId = "sidebyside"
    , leftId = "leftid"
    , rightId = "rightid"
  ) %>%
  addLayersControl(overlayGroups = c("april", "may")) %>%
  addControl(htmltools::HTML("April 2020"), position = "bottomleft") %>%
  addControl(htmltools::HTML("May 2020"), position = "bottomright")

[tim-salabim]

Using this approach for RGB images. This should be an easy way to use in leafem::add*RGB functions. Maybe we should simply have one addRGB as a generic and have methods for Raster* & stars objects?

NOTE: for this example, we need to read/write via stars as passing the tif directly to addGeotiff places it somewhere in the middle of the Pacific Ocean (likely an issue with the UTM converter in georaster-layer-for-leaflet)

library(leaflet)
library(leafem)
library(stars)

tif = system.file("tif/L7_ETMs.tif", package = "stars")
x1 = read_stars(tif)
tmpfl = tempfile(fileext = ".tif")
write_stars(st_warp(x1, crs = 4326), tmpfl)

myCustomJSFunc = htmlwidgets::JS(
  "
    pixelValuesToColorFn = (raster, colorOptions) => {
    // helpers from https://stackoverflow.com/questions/5623838/rgb-to-hex-and-hex-to-rgb
      function componentToHex(c) {
        var hex = c.toString(16);
        return hex.length == 1 ? '0' + hex : hex;
      }
      
      function rgbToHex(r, g, b) {
        return '#' + componentToHex(r) + componentToHex(g) + componentToHex(b);
      }

      var pixelFunc = values => {
        if (isNaN(values[0])) return colorOptions.naColor;
        return rgbToHex(values[3], values[2], values[1]);
      };
      return pixelFunc;
    };
  "
)

leaflet() %>%
  addProviderTiles("Esri.WorldImagery") %>%
  addGeotiff(
    file = tmpfl
    , opacity = 0.9
    , colorOptions = colorOptions(
      , na.color = "transparent"
    )
    , pixelValuesToColorFn = myCustomJSFunc
    , group = "testRGB"
  ) %>%
  addLayersControl(overlayGroups = "testRGB")

[tim-salabim]

So, for baisc arithmetic operations on bands we can use something like the following approach:

library(leaflet)
library(leafem)
library(stars)

tif = system.file("tif/L7_ETMs.tif", package = "stars")
x1 = read_stars(tif)
tmpfl = tempfile(fileext = ".tif")
write_stars(st_warp(x1, crs = 4326), tmpfl)

pal = hcl.colors(256, "viridis")

# define band arithmetic as we would in R
ndvi = "(4 - 3) / (4 + 3)"

# convert to JS by subtracting 1 and pasting in the anonymous function arg `values`
bandCalc = function(band_calc) {
  idx_r = gregexpr("[0-9]+", band_calc)
  js_bands = as.numeric(unlist(regmatches(band_calc, idx_r))) - 1
  
  js_band_calc = gsub("[0-9]+", paste0("values[", "%s", "]"), band_calc)
  js_band_calc = do.call("sprintf", c(list(js_band_calc), js_bands))
  return(js_band_calc)
}

bandCalc(ndvi)
# [1] "(values[3] - values[2]) / (values[3] + values[2])"

myCustomJSFunc = htmlwidgets::JS(
  sprintf(
    "
      pixelValuesToColorFn = (raster, colorOptions) => {
        const cols = colorOptions.palette;
        var scale = chroma.scale(cols);
  
        if (colorOptions.breaks !== null) {
          scale = scale.classes(colorOptions.breaks);
        }
        var pixelFunc = values => {
          var vals = %s;
          let clr = scale.domain(colorOptions.domain);
          if (isNaN(vals)) return colorOptions.naColor;
          return clr(vals).hex();
        };
        return pixelFunc;
      };
    "
    , bandCalc(ndvi)
  )
)

leaflet() %>%
  addProviderTiles("Esri.WorldImagery") %>%
  addGeotiff(
    file = tmpfl
    , opacity = 0.9
    , colorOptions = colorOptions(
      palette = pal
      , na.color = "transparent"
      , domain = c(-0.8, 1)
    )
    , pixelValuesToColorFn = myCustomJSFunc
    , group = "testNDVI"
  ) %>%
  addLayersControl(overlayGroups = "testNDVI")

which simply converts the character representation of the band calculation ndvi = "(4 - 3) / (4 + 3)" into a valid JS representation by converting it to a 0-based index and inserting the values[] object to be indexed at the correct position. This should work for all calculations that use +, -, *, /, ( and ) which is nice, but rather limited... It means that the R user only needs to provide a band arithmetic as it would be valid in R, bandCalc() would take care of the rest. If you want to play around with the solution above, it should be enough to change the ndvi = "(4 - 3) / (4 + 3)" definition and the domain = c(-0.8, 1) in the colorOptions call in addGeotiff to account for the min/max values that the calculation may produce.

@timelyportfolio @edzer @lbusett how do we feel about an API like this? And how could we extend this to more complicated calculations involving function based calculations like max, sqrt, sd, var, mean, median, abs, ... to only name a few that are likely used in often in the remote sensing realm?

@timelyportfolio do you think it is realistic to come up with a general translator between R arithmetic expressions and their JS counterpart? I.e. can we come up with functionality translate at least the most common functions?

[tim-salabim]

to be continued in #29