Mapview fails over "stars" proxy objects

[lbusett]

I noticed that using mapview over a "proxy" stars object currently fails:

f <- system.file("external/test.grd", package="raster")
r <- stars::read_stars(f, proxy = T)
mapview(r)

 Error in (function (classes, fdef, mtable)  : 
  unable to find an inherited method for function ‘mapView’ for signature ‘"stars_proxy"’ 

Therefore, also this fails:

r <- raster(f)
r <- stars::st_as_stars(r)
mapview(r)

 Error in (function (classes, fdef, mtable)  : 
  unable to find an inherited method for function ‘mapView’ for signature ‘"stars_proxy"’ 

Consider that I fear therefore that mapview will always fail over large stars objects unless "proxy = FALSE" is specified explicitly.

[edzer]

I can't reproduce your first code snippet:

> f <- system.file("external/test.grd", package="raster")
> r <- stars::read_stars(f, proxy = T)
trying to read file: /home/edzer/R/x86_64-pc-linux-gnu-library/4.0/raster/external/test.grd
Error in CPL_read_gdal(as.character(x), as.character(options), as.character(driver),  : 
  file not found
In addition: Warning message:
In CPL_read_gdal(as.character(x), as.character(options), as.character(driver),  :
  GDAL Error 1: Unhandled datatype=FLT4S 

but can with the second, when f is a file that can be read. It is a bug in stars, which always wanted to return a proxy object if a Raster object had a file or filename slot filled (unequal to ""). This should now be corrected. st_as_stars.Raster still returns a proxy object if the object is big, but this can be overridden by proxy=FALSE.

[lbusett]

@edzer thanks!

@tim-salabim: I was thinking that mapview could maybe check if a "proxy" stars object was passed and convert it automatically to a "standard" one to avoid problems (Although it would require loading the entire raster, so maybe is not a great idea). @edzer Is there a simple way to do that ?

PS: Strange - this works for me:

> f <- system.file("external/test.grd", package="raster")
> r <- stars::read_stars(f, proxy = T)
> r
stars_proxy object with 1 attribute in file:
$test.grd
[1] "[...]/test.grd"

dimension(s):
  from  to offset delta                       refsys point values    
x    1  80 178400    40 +proj=sterea +lat_0=52.15...    NA   NULL [x]
y    1 115 334000   -40 +proj=sterea +lat_0=52.15...    NA   NULL [y]

Maybe a different raster version? I am currently on the latest github one. although test.grd has been there for a long time.

[edzer]

st_as_stars will work for that, but didn't work (and currently won't work without proxy=FALSE) for huge Raster objects. You would need

r = st_as_stars(st_as_stars(raster(f))) # in case f is a huge file

or directly

st_as_stars(raster(f), proxy=FALSE)

I don't like this, but see no simple alternative ATM.

For mapviewing massive rasters, you'd like a pattern that re-reads and rerenders when you zoom & pan, but that's not in place yet. What you can do now is st_as_stars(stars_proxy_obj, downsample = x) to downsample, or call read_stars with RasterIO parameters that specify a subregion and sampling rate.

[lbusett]

What you can do now is st_as_stars(stars_proxy_obj, downsample = x) to downsample, or call read_stars with RasterIO parameters that specify a subregion and sampling rate.

AFAICT, this is what tmap is currently doing: (https://github.com/mtennekes/tmap/blob/626177dbee57ca62050f45cc798b5f79e1285266/R/misc_stars.R)

Seems reasonable to me. Downsampling could be governed by the already available maxpixels argument.

[tim-salabim]

https://github.com/GeoTIFF/georaster-layer-for-leaflet/tree/master

seems to have capabilities that could make the whole raster experience with leaflet in R a lot nicer, though I'm not too convinced about using KNN for on-the-fly interpolation. I'll try to mash-up a test repo for this when I find the time. Currently trying to get mapview et. al. in good shape for a summer school :-)

[tim-salabim]

Seems reasonable to me. Downsampling could be governed by the already available maxpixels argument.

For raster we use sampleRegular I think using maxpixels

[lbusett]

I can see if i can manage a PR on this in the coming days, if you wish.

[tim-salabim]

Of course. That would be much appreciated

[edzer]

@tim-salabim the georaster-layer-for-leaflet uses simple nearest neighbor interpolation, not KNN (well, K=1, yes); this is also what stars' downsample uses or default for RasterIO stuff.

[tim-salabim]

@edzer thanks for the clarification! And good to know that it's somewhat aligned. I guess it's a matter of getting used to that the raster constantly changes appearance when you zoom in...

[edzer]

If you resample to a cell size that corresponds to screen pixels it should look "natural".

[tim-salabim]

Not sure if that can be set. I will investigate when I find the time

[lbusett]

@tim-salabim Having a look at the issue. "Dispatching" mapview also over proxy objects is straightforward, and then as a start even something simple like this at the beginning of the "stars" method could work:

  # check if downsampling is needed
  n_pixels <- dim(x)[1] * dim(x)[2]
  if (n_pixels > maxpixels) {
    do.downscale = TRUE
  } else {
    do.downscale = FALSE
  }
    
  if (inherits(x, "stars_proxy")) {
    # convert to proper "stars" if proxy and downscale if needed
    if (!do.downscale) {
      x <- st_as_stars(x)
    } else {
      # compute the resampling factor
      resamp_factor <- 2 # (TODO : function to compute this  - approach can be "stolen" from tmap)
      x <- stars::st_as_stars(x, downsample = resamp_factor)
    }
  } else {
    # how to downscale a non-proxy object ? 
  }

Only inconvenient could be that the "downsample" argument in st_as_stars only affects proxy objects. Therefore, large "standard" stars images would not be downsampled.

@edzer Am I missing something here ? Is there a simple way to easily downsample also a non-proxy stars raster?

[edzer]

Yes, try

tif = system.file("tif/L7_ETMs.tif", package = "stars")
plot(stars:::st_downsample(read_stars(tif), c(10,10,1)))

st_downsample is not exported, but we could do that if mapview would need it - within stars, it's only used in the plotting routines; @mtennekes - tmap may also want it?

[tim-salabim]

@lbusett just so you know, I am currently trying to add georaster-for-leaflet support to leafem. If that works nicely, the downsampling approach may not warrant the effort. Not sure what to advise as we will need to see how the performance is, but I don't want you to put work into something that may not be needed.

[mtennekes]

Yes, try

tif = system.file("tif/L7_ETMs.tif", package = "stars")
plot(stars:::st_downsample(read_stars(tif), c(10,10,1)))

st_downsample is not exported, but we could do that if mapview would need it - within stars, it's only used in the plotting routines; @mtennekes - tmap may also want it?

Yes, please! Actually, I'm already using it as a hard copy. Also happy with exports of: select_band, is_regular_grid, has_raster, has_rotate_or_shear, is_curvilinear, is_rectilinear, regular_intervals and get_downsample ;-)

[tim-salabim]

So, I've had success to get it to work. See https://twitter.com/TimSalabim3/status/1280950560892518402 (gif too large for github).

The code for the gif:

library(leaflet)
library(leafem)

leaflet() %>%
  addProviderTiles("Esri.WorldImagery") %>%
  leafem:::addGeoRaster(
    file = "/home/timpanse/software/data/srtm_39_03.tif"
    , group = "test"
    , layerId = "testid"
  )

Currently, only supports single layer file (tested) and url (not tested). stars/raster support will basically be either passing on the file path or write to tempfile and use that.

This is only a POC, but it seems reasonably performant, so I have hopes. Much of the performance seems to depend on the resolution which is currently hardcoded to be 96 px.

Feel free to try yourselves with your own files. I think they need to be in 4326 but I am not sure (haven't tested others).

EDIT: Here's the commit r-spatial/leafem@7ae2b86

[lbusett]

This seems very promising! Thanks!

Concerning projection, 4326 does not seem mandatory. At least UTM projection works:

f <- system.file("external/test.grd", package="raster")
a <- stars::read_stars(f) %>% stars::st_warp(crs = 32632, cellsize = 2)
tmpfile <- tempfile(fileext = ".tif")
stars::write_stars(a, tmpfile)

leaflet() %>%
  addProviderTiles("Esri.WorldImagery") %>%
  leafem:::addGeoRaster(
    file = tmpfile,
    , group = "test"
    , layerId = "testid"
  )

However, this doesn't:

f <- system.file("external/test.grd", package="raster")
leaflet() %>%
  addProviderTiles("Esri.WorldImagery") %>%
  leafem:::addGeoRaster(
    file =tmpfile,
    , group = "test"
    , layerId = "testid"
  )

, though it may be related to the projection of the original file (strereographic something)

Concernig the issue about stars proxy objects / downsampling, I'll let you decide: I could come up with something quick to fill the void while you develop the "new" technique, or I can stay like this. Just let me know.

[tim-salabim]

Thanks! UTM seems to be nativley supported, see https://github.com/GeoTIFF/georaster-layer-for-leaflet/blob/master/utils/utm.js

Re stars downsampling, I would be grateful for a quick solution as I don't know how long I will need to get addGeoRaster to a stable, usable state. I am a bit pressed for time as I need to get both leafem and mapview updated on CRAN before mid August and there are still a few kinks to be ironed out. Hence, a PR with would be appreciated.

[lbusett]

Ok, I'll see what I can do !

[tim-salabim]

I'll try to implement a more robust POC for addGeoRaster tonight, so we can play a bit more with it, especially passing options like opacity, resolution, etc...

One thing that may turn into a can of worms is the fact that we need/are able to handle coloring on the javascript side. This can be very powerful, as it is not as static as creating color vectors in R and passing them. On the other hand, we need to carefully think about what we need to pass from R and how flexible that can be (e.g. palettes created with colorRampPalette).

The noce thing about being able to handle all sorts of things on the JS side is that we can also do band calculations on the fly. See https://github.com/GeoTIFF/georaster-layer-for-leaflet/blob/master/tests/ndvi.html#L37-L82 for an example. But again, this diverts far enough from what I've implemented in any method so far, that I cannot foresee how much work this will entail.

But I am quite excited about all the possibilities :-)

[edzer]

Is this going towards some kind of in-browser google earth engine?

[tim-salabim]

I honestly don't know where this will lead. Though I have no intentions to reinvent existing wheels. This is just an attempt to adress the long standing issue of lacking capabilities to interactively view large raster data from R. Given that without shiny we have no way of communicating anything back to the R session once things are rendered I think it's good that calculations like downsampling happen when they need to happen, i.e. on zoom/pan changes. And this is what leafem::addGeoRaster is aiming at.

However, given that with georaster-layer-for-leaflet we have a tool to do this, I guess it's only natural to try out other calculations too. But I would not overload any of these as my feeling is that this will degrade performance quickly if calculations become too involved. NDVI is easy and I think that's why it's included as an example, for more computationally heavy stuff, we'd need to resort to something like GPU.js (at which stage I guess one could argue that we're headed towards a google earth engine type thing).

We'll see how this evolves. But now that you have me thinking about this, I guess it's not a bad thing to at least try and implement something along the way of reactable's ability to pass javascript functions for cell value rendering. If a user knows JS, why not give them the tools to make use of it.

[tim-salabim]

But if we wanna go this way, that's the next step I guess: https://github.com/GeoTIFF/geoblaze

[lbusett]

Concerning this:

Thanks! UTM seems to be nativley supported, see https://github.com/GeoTIFF/georaster-layer-for-leaflet/blob/master/utils/utm.js

Re stars downsampling, I would be grateful for a quick solution as I don't know how long I will need to get addGeoRaster to a stable, usable state. I am a bit pressed for time as I need to get both leafem and mapview updated on CRAN before mid August and there are still a few kinks to be ironed out. Hence, a PR with would be appreciated.

1. I can put out something very quickly to allow rendering/downsampling of stars proxy.

2. Regarding downsampling of standard stars objects, as @mtennekes I think that having st_downsample exported would greatly facilitate things.

@tim-salabim If you agree, I would make a PR to implement 1). Then I'd wait for @edzer to export st_downsample to implement the rest. I do not think that you wish to rely on a github version, so I think this could/should wait that the "exported st_downsample " hits CRAN. If however you are ok on relying on a github version for this (we could make it "conditional" on pkg version), or also about using a non exported function via ::: just tell me.

[tim-salabim]

Thanks @lbusett! I think approach 1) i.e. only implement for proxy for now would be best. I agree that there is no rush, so we should wait until st_downsample is available on CRAN.

On that note, I've just used st_downsample to play around a bit more with addGeoRaster. Had to downsample a global elevation grid from

dim(tst)
    x     y 
32727 15474 

to

dim(tst_ds)
    x     y 
21818 10316

to be able to render it. But the experience once rendered is great. I've exposed resolution and opacity for people to play with. Feed back very welcome! Next is color handling, but for that I need some proper programming time (not in the pub :-) )

[lbusett]

Ok, PR submitted - sorry for some mess with indentation.

Now we get for example on proxy objects:

f <- system.file("external/test.grd", package="raster")
r <- stars::read_stars(f, proxy = T)
mapview(r, maxpixels = 1000)

Number of pixels is above 1000.Only about 1000 pixels will be shown.
You can increase the value of `maxpixels` to 9200 to avoid this.

On non proxy ones:

f <- system.file("external/test.grd", package="raster")
r <- stars::read_stars(f, proxy = F)
mapview(r, maxpixels = 1000)

Number of pixels is above 1000. Automatic downsampling of `stars` object is not yet implemented, so rendering may be slow.
You can pass a `stars` proxy object to mapview to get automatic downsampling.

@mtennekes : I basically replicated your approach for computation of the downsamplig factor - hope you don't mind!

[tim-salabim]

Continuing the addGeoRaster stuff here:

r-spatial/leafem#25

Closing this, as the solution proposed by @lbusett has now been merged into develop.