raster fails rendering if not already in 3857 projection

[lbusett]

Hi,

working on Windows, with tmap and tmaptools at v3.0, sf 0.9-2 installed this morning from github, when tmap_mode is "view", tmap fails rendering raster data unless they are already in 3857 projection. In practice:

this works

s = stars::read_stars(system.file("tif/L7_ETMs.tif", package = "stars"))[,,,1]
tmap_mode("plot")
tm_shape(s) + tm_raster()

this fails - keeps running endlessly

tmap_mode("view")
tm_shape(s) + tm_raster()

this works

s <- raster::raster(system.file("tif/L7_ETMs.tif", package = "stars"))
s <- raster::projectRaster(s, crs = "+init=epsg:3857", res = c(28.5, 28.5))
tm_shape(s) + tm_raster()

My guess is that you may need to set "project = TRUE" in this call to "leafem::addStarsImage":

https://github.com/mtennekes/tmap/blob/dfea5be4d1ba76b57daabaf497b30a0aab4c0080/R/view_tmap.R#L543

, or add an automatic reprojection beforehand ?

HTH,

Lorenzo

PS: as always, thanks for your work on tmap

[lbusett]

UPDATE

I had a better look, and saw that the culprit seems to be here:

https://github.com/mtennekes/tmap/blob/3b76d2e544d512d343d83719e1d041fea38663c6/R/preprocess_shapes.R#L85

, which returns a curvilinear-grid stars object.

Modifying the line to:

shp <- stars::st_warp(shp, crs = 3857) #transform to 3857

seems to solve the issue. I could do a PR if you wish. (probably also line

https://github.com/mtennekes/tmap/blob/3b76d2e544d512d343d83719e1d041fea38663c6/R/preprocess_shapes.R#L89

should be changed accordingly ?

[mtennekes]

Thanks Lorenzo!

It is so indeed very slow because it creates a curvilinear grid, where grid cells are converted to polygons (> 100,000 in this case).

We could indeed do a warp, which returns a regular raster, but with a new grid. For some applications (like this) this is preferable, but for other applications a transformed grid may be better. E.g. smaller rasters or applications where it is important to plot without resampling (at the cost of being slow).

It's up to the user which method to use, but the question for us is how. What we could do is add an argument in tm_shape which defines whether a grid is warped or transformed, with a default setting depending on the size of the raster. Currently, there is a tmap option max.raster, which defines above which resolution a raster is downscaled. I just found out that this option is not implemented yet in tmap 3.0. Maybe there is less need, since users can use starts proxy objects. We could add an option warp.raster that defines above which raster size it is warped by default?

What do you think?

[lbusett]

Hi @mtennekes

thanks for replying.

I think giving the "choice" to the user could be a good idea. To simplify things, I would propose to just use st_warp as default, giving the user the choice to eventually "activate" the "transform" route. In my opinion, due to the fact that the "view" mode of tmap is generally used for "quick" interactive exploration of datasets, it is important that the default is reasonably fast. At the moment, even quite small raster datasets based on Sentinel-2 imagery (10 m resolution), that used to render reasonably fast on tmap 2.3, lead to apparent unresponsiveness. However, also having a default based on dimension would be ok for me, but I'd suggest the "maximum size" to be kept quite small.

Concerning the "max.raster" option: wasn't that already available on older "tmap" versions? Was something changed in 3.0 ?

Lorenzo

[mtennekes]

Sounds reasonable. @Nowosad @zross @edzer Any preferences from your side? warp or transform?

Concerning the "max.raster" option: wasn't that already available on older "tmap" versions? Was something changed in 3.0 ?

Yes, under the hood, tmap switched from raster to stars (the whole reason of 2.x -> 3.x).

[lbusett]

@mtennekes

Yes, under the hood, tmap switched from raster to stars (the whole reason of 2.x -> 3.x).

Yes, I noticed that. So, the max.raster option is currently ignored? Can I ask what tmap is currently doing for "large" rasters? Does it use automatic subsampling based on stars proxy (as in https://cran.r-project.org/web/packages/stars/vignettes/stars2.html) ?

[Nowosad]

@mtennekes, I think we have a few issues here:

1. Data type - raster vs vector; I assume we are focussing on the former here
2. Raster type - regular grids, rotated and sheared grids, rectilinear grids, curvilinear grids; the question is which of those are (will be) supported by tmap?
3. Input type - e.g. stars vs stars.proxy - should the behaviour of tmap be identical for both or should it be different?
4. Raster size - what is the threshold (number of cols and rows) when we want to decrease a resolution of the input raster?

[mtennekes]

It should be improved now (commits https://github.com/mtennekes/tmap/commit/a656d63660174617a62b18bb9e72d873db2d750a and https://github.com/mtennekes/tmap/commit/eb78a1ddbcbbc86f60532f0c15859bd019d6a594).

See also the test script: https://github.com/mtennekes/tmap/blob/master/sandbox/stars.R

• After testing different stars, I noticed that a warp is (almost) always preferable over a transform. So I followed @lbusett suggestion. Question: do we still want users to be able to transform stars after reprojection? If so, I will add an argument to tm_shape to specify this.
• max.raster is working again. I have set it so 1e6, which is similar to a 1000 by 1000 raster. This is much more responsive than the original value 1e7 (for plot mode), which is about 3200 by 3200. Do you agree? Should we add an argument that specifies whether or not to downsample the stars object? If so, I will add an arugment to tm_shape for this as well.

@Nowosad:

1. Yes, indeed. I haven't thought about an upper limit for the number of sf features. Could be worthwhile to think about...
2. All should be working, if not let me know! See the test script (link above).
3. Yes, they should be identical. If it is a stars proxy object, it will resize to max.raster cells. If it is a stars object, it will downscale to max.raster cells when it has more cells.
4. Currently, max.raster is set to 1e6 for both plot and view mode. As said, if the total number of cells is higher it will downscale. It takes the aspect ratio into account. For instance a 2:1 raster will be downscaled to 1415 by 707.

Suggestions and feedback are welcome.

[lbusett]

@mtennekes

seems good to me! Thanks ! I reply inline on some points:

After testing different stars, I noticed that a warp is (almost) always preferable over a transform. So I followed @lbusett suggestion. Question: do we still want users to be able to transform stars after reprojection? If so, I will add an argument to tm_shape to specify this.

Do you mean if it should be possible to be able to switch to the "transform" route? I'd guess so, if it is not much of a problem. Showing the true "footprints" of original raster pixels is in some cases necessary (in particular for coarse resolution grids), in my opinion. Regarding the slowness previosly experienced, I was wondering if maybe it could be reduced by switching to the new "leafgl" feature of mapview.

Should we add an argument that specifies whether or not to downsample the stars object? If so, I will add an arugment to tm_shape for this as well.

I think it could be a good idea. Also, sending a message to the user when downsampling occurs could be useful.

All should be working, if not let me know! See the test script (link above).

Most examples of the test script are working. I see however strange results at line 74:

[edzer]

That is the result from

st_warp(st_set_crs(r, 4326), crs = 3857)

which assumes WGS84 for r, then creates a regular grid from a rectilinear grid in web Mercator. I think this case illustrates that it would be good to be able to control (some of) the parameters that go into st_warp from tmap; this for instance looks better:

a = sum(st_area(st_set_crs(r, 4326))[[1]]) %>% units::drop_units()
b = st_warp(st_set_crs(r, 4326), crs = 3857, cellsize = sqrt(a)/100)
tm_shape(b) +
	tm_raster(style = "cat", palette = "cat") +
tm_shape(rsf) +
	tm_borders()

[mtennekes]

raster.downsample and raster.warp added to tm_shape. Some (messy) examples can be found in https://github.com/mtennekes/tmap/blob/master/sandbox/stars.R

I think this case illustrates that it would be good to be able to control (some of) the parameters that go into st_warp from tmap; this for instance looks better:

The grid warp only seems to go wrong for small stars, like this one. Extra arguments, or arguments that can be passed on to st_warp could be nice to have too, but 99% of the users just want good out-of-the-box maps. I'm thinking about an (automatic) method to improve the quality of warped grids for small stars. @edzer do you have any tricks how to do this? For instance, something like: if the number of cells is less than 100, set cellsize of st_warp to ... (e.g. your example: sqrt area / 100). I also saw there is an argument segments, but I'm not sure what is does, and how it relates to cellsize.