Tmap: tm_raster - a number of pixels plotted

[Nowosad]

tm_raster() is fast for small and medium-sized rasters (see the last example).
However, it is not so good for large ones (>few millions of pixels) - then it takes a long time (and an enormous part of RAM) to make a plot.
Is it possible to add a new argument similar to maxpixels in rasterVis?

# prepare data ------------------------------------------------------------
library(raster)
r = raster(nrow = 1000, ncol = 1000, vals = sampleInt(10, 1e+6, replace = TRUE))

# rasterVis ---------------------------------------------------------------
library(rasterVis)
system.time(print(levelplot(r, maxpixels = 1e6)))

#>    user  system elapsed 
#>   5.744   0.103   5.937
system.time(print(levelplot(r, maxpixels = 1e4)))

#>    user  system elapsed 
#>   0.506   0.005   0.517
system.time(print(levelplot(r, maxpixels = 1e1)))

#>    user  system elapsed 
#>   0.093   0.000   0.094

# tmap --------------------------------------------------------------------
library(tmap)
system.time(print(tm_shape(r) + tm_raster()))
#> Linking to GEOS 3.6.1, GDAL 2.2.4, proj.4 4.9.3

#>    user  system elapsed 
#>   0.499   0.028   0.530

Created on 2018-07-19 by the reprex package (v0.2.0).

[mtennekes]

Definitely a great idea!

The implementation is not as straightforward as it may sound. Currently the processing of raster objects is:

1. determine type of raster data
2. reproject if necessary
3. extract and preprocess raster data (e.g. use color tables, and fix factor levels if step 2 is applied)
4. process raster data (apply class intervals etc.)
5. plot the raster object

It would make sense to apply the aggregation between 2 and 3. However, this would influence the result of step 4, i.e. the class intervals may be different.

[tim-salabim]

Different from what? The original that was passed in, or the one you see in the plot. I think it is more important that things are coherent within the plot. A suitable warning that you are messing with the data should be enough for the user (this is something that caught me multiple times when working with raster plots - as there is no warning).

Just to provide some perspective, in mapview we do it between 1. and 2. using this function. Whether it is done before reprojection does not matte I think, but it should definitely happen before any color creation etc. Additionally, reprojection takes a lot longer for large rasters, that's why we opted to sample first...

[mtennekes]

Thx Tim! I followed your suggestion.

max.raster is added as tmap option:

> tmap_options("max.raster")
$max.raster
 plot  view 
1e+07 1e+06 

The default values are chosen arbitrarily. According to these values, the maximum raster in plot mode is around 3000 x 3000, and for interactive mode 1000 x 1000. Hopefully, it makes sense.

[Nowosad]

Awesome Martijn!

Small benchmark of this option:

# prepare data ------------------------------------------------------------
library(raster)
#> Loading required package: sp
r = raster(nrow = 3000, ncol = 3000, vals = sampleInt(10, 9e+06, replace = TRUE))

# rasterVis ---------------------------------------------------------------
library(tmap)
tmap_options("max.raster")
#> $max.raster
#>  plot  view 
#> 1e+07 1e+06
system.time(print(tm_shape(r) + tm_raster()))
#> Linking to GEOS 3.6.1, GDAL 2.2.4, proj.4 4.9.3

#>    user  system elapsed 
#>   5.311   1.133   6.463

tmap_options(max.raster = c(plot = 1e+04))
system.time(print(tm_shape(r) + tm_raster()))
#> Raster object has 9e+06 (3000 by 3000) cells, which is larger than 10000, the maximum size determined by the option max.raster. Therefore, the raster will be shown at a decreased resolution of 10000 cells. Set tmap_options(max.raster = c(plot = 9e+06, view = 9e+06)) to show the whole raster.

#>    user  system elapsed 
#>   2.831   0.100   2.942

tmap_options(max.raster = c(plot = 1e+02))
system.time(print(tm_shape(r) + tm_raster()))
#> Raster object has 9e+06 (3000 by 3000) cells, which is larger than 100, the maximum size determined by the option max.raster. Therefore, the raster will be shown at a decreased resolution of 100 cells. Set tmap_options(max.raster = c(plot = 9e+06, view = 9e+06)) to show the whole raster.

#>    user  system elapsed 
#>   0.255   0.052   0.308

Created on 2018-07-31 by the reprex package (v0.2.0).