Support for `stars` format for raster layers?

[cboettig]

Apologies if there are good reasons this doesn't make sense, but I love the way tmap makes it easy to make thematic maps that mix raster and vector data (including support for sf) and very intrigued by a lot of the functionality of stars in working with rasters. Is this in scope for support in tmap?

Thanks for considering!

[mtennekes]

The first stars are supported!

library(tmap)
library(stars)
#> Loading required package: abind
#> Loading required package: sf
#> Linking to GEOS 3.6.2, GDAL 2.2.3, PROJ 4.9.3
tif = system.file("tif/L7_ETMs.tif", package = "stars")
x = read_stars(tif)
qtm(x)

tm_shape(x) + tm_rgb()

^{Created on 2020-01-09 by the reprex package (v0.3.0)}

[mtennekes]

This is a very lazy solution, basically done with as(x, "Raster"). It does not support curvilinear grids yet.

The next (bigger) step will be to support stars natively. The question is how. I think there are two options:

• Add stars processing alongside the existing raster (class) processing.
• Replace raster processing with stars processing. I.e. raster objects are transformed to stars.

There is also terra (https://github.com/rspatial/terra), which will be the replacement of raster. It is not an option to transform stars to terra objects, since terra will probably not support curvilinear grids.

Any thoughts? @Robinlovelace @Nowosad @edzer

[tim-salabim]

My experience is that stars processing is faster than Raster processing, so from that point of view option 2 is something to consider.

[edzer]

I think the stars data model includes that of Raster. terra is not on CRAN, and I haven't seen any ETA for CRAN submission; the GH page mentions an alpha release. AFAIK its data structure will resemble that of Raster (i.e., raster maps or a stack of them).

[edzer]

@mtennekes is this in a branch? After installing tmap from master I see

> library(tmap)
> library(stars)
Loading required package: abind
Loading required package: sf
Linking to GEOS 3.8.0, GDAL 3.0.2, PROJ 6.2.1
> #> Loading required package: abind
> #> Loading required package: sf
> #> Linking to GEOS 3.6.2, GDAL 2.2.3, PROJ 4.9.3
> tif = system.file("tif/L7_ETMs.tif", package = "stars")
> x = read_stars(tif)
> qtm(x)
Error in get_projection(mshp_raw, output = "crs") : 
  shp is neither a sf, sp, nor a raster object
> sessionInfo()
R version 3.6.2 (2019-12-12)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 18.04.3 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/atlas/libblas.so.3.10.3
LAPACK: /usr/lib/x86_64-linux-gnu/atlas/liblapack.so.3.10.3

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] stars_0.4-1 sf_0.8-0    abind_1.4-5 tmap_2.3-2 

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.3         compiler_3.6.2     later_1.0.0        RColorBrewer_1.1-2
 [5] class_7.3-15       tools_3.6.2        digest_0.6.23      viridisLite_0.3.0 
 [9] lattice_0.20-38    rlang_0.4.2        shiny_1.4.0        DBI_1.1.0         
[13] crosstalk_1.0.0    parallel_3.6.2     rgdal_1.4-8        fastmap_1.0.1     
[17] e1071_1.7-3        raster_3.0-7       htmlwidgets_1.5.1  rgeos_0.5-2       
[21] classInt_0.4-2     leaflet_2.0.3      grid_3.6.2         R6_2.4.1          
[25] XML_3.98-1.20      sp_1.3-4           magrittr_1.5       promises_1.1.0    
[29] tmaptools_2.0-2    codetools_0.2-16   leafsync_0.1.0     htmltools_0.4.0   
[33] units_0.6-5        dichromat_2.0-0    mime_0.8           xtable_1.8-4      
[37] httpuv_1.5.2       KernSmooth_2.23-15 lwgeom_0.1-7      

[mtennekes]

Sorry, forgot to mention that the GitHub version of tmaptools is also required.

When I will to the big raster to stars transition I'll start a new branch.

[edzer]

Confirmed!

More urgently I need your review of https://github.com/r-spatial/sf/issues/1225. Notably, tmap hard-codes crs objects, making it for sf impossible to change their structure, which is needed. E.g. in tmaptools, create_crs should use sf::st_crs always and never

get_proj4_code.R:        structure(list(epsg = as.integer(NA), proj4string = x), class = "crs")

More in general tmap seems to hinge entirely on proj4strings, which are now legacy, and no longer the recommended way to describe a CRS (or transformation).

[Nowosad]

@mtennekes, great news!

Some of my comments:

1. CRC issue is probably the priority now.
2. stars objects could represent both raster and vector data (cubes). Of course, it would be great to start somewhere (e.g., simple regular rasters), but it is crucial to have it at the back of the head that there is more to it than simple rasters...
3. Rasters in stars could have many forms (regular, rotated, sheared, rectilinear, curvilinear rasters)
4. One of the great features in stars is described at https://r-spatial.github.io/stars/articles/stars2.html#stars-proxy-objects - "only those pixels are read that can be seen on the plot". It makes visualizations for large objects faster and less memory intensive. It would be great to have something similar in tmap.
5. Another aspect of raster data cubes is its dimensionality. stars objects could have many attributes, where each attribute can have several bands plus timestamps, etc. It would be essential to decide what is the default action for stars objects and how to specify other actions (e.g., one band only, one attribute only, etc.)

I hope that my comments are useful here. Best, J.

[edzer]

Thanks @Nowosad ! Yes, two areas where tmap could shine that would benefit from what stars supports but raster doesn't are:

• (time) sequences of vector maps, e.g. shown at https://r-spatial.github.io/stars/index.html
• two-dimensional layouts of raster or vector maps; for raster e.g. a time series (horizontal) of different spectral bands (vertical)

Another one, pretty big, is that of plotting massive raster images at screen resolution. This turned out to be pretty much impossible by design to do with ggplot2, and it would be great if tmap could do this. See e.g. https://r-spatial.github.io/stars/articles/stars2.html#stars-proxy-objects Package starsdata is installed by

install.packages("starsdata", repos = "http://gis-bigdata.uni-muenster.de/pebesma/", type = "source")

Let me know if you need examples and/or code snippets.

[tim-salabim]

Reagrding @edzer 's first point from an interactive point of view (i.e. in tmap_mode("view")), @timelyportfolio 's leaftime package is now on CRAN. I've been meaning to look into it for supporting space-time data in mapview, but got side-tracked integrating flatgeobuf first. This is next on my list, though!

[mtennekes]

Thanks for all inputs.

More and more stars can be plotted with tmap now: see the demo script https://github.com/mtennekes/tmap/blob/stars/sandbox/stars.R , which will be extended step by step.

• Forms of the grid: sheared/rotated/rectilinear/curvilinear. They work will in plot mode, but in view mode, there are still problems because leaflet expects epsg 3857.
• Dimensionality: still working on that. Selecting attributes or bands works (e.g. the land object is now a multi-attribute stars object). I still have to think about how to deal with timestamps. However, I think this should be relatively easy to implement.
• The final big step is how to plot massive stars. Currently, the non-regular stars are cast to sf objects, simply because grid::rasterGrob expects regular rasters.

If any of you has examples that illustrate the features of stars (other than mentioned in the stars-vignettes), please let me know, or add to https://github.com/mtennekes/tmap/blob/stars/sandbox/stars.R.

[edzer]

I get

> qtm(s)
Error in if (!tmaptools::is_projected(shp)) { : 
  missing value where TRUE/FALSE needed

you need to use isTRUE() around st_is_longlat if you want a logical result.

[edzer]

Ah, installing stars branch of tmap resolved that.

[edzer]

Great progress, @mtennekes ! About the massive raster data, the approach would be probably to create a raster tile server on-the-fly, and have that serve the pixel values. That should do an st_warp() to web mercator, first, for this purpose, I guess, or even do that on the fly.

[tim-salabim]

Heads-up, I've implemented that once in mapview using gdal2tiles.py via system call, but it was very slow. Too slow for on-the-fly creation for the scope of mapview. The code for this can be found here:

https://github.com/r-spatial/mapview/blob/c0cf7d2071aafc405bbec3cfc2bc572c1c07caa8/R/extensions.R#L1075-L1148

Note, that for serving data from a tile folder, we already have leafem::addTiles().

I've played around a little bit with libvips which seems to be sufficiently fast for on-the-fly creation, but I don't think there's an R biding for it + plus I have yet to figure out how to properly render the resulting tiles on a leaflet map.

I'd be very interested in realising something that is fast enough for on-the-fly creation, as it's been nagging me for some time now that we still don't have proper ways to display large rasters quickly.

[edzer]

I was thinking about creating tiles on-the-fly at request time, your solution seems to compute them prior to serving them is that right?

[tim-salabim]

Yes, it first writes all tiles to disk and then invokes the rendering engine leaflet. Are you aware of any (JavaScript) tools to create tiles on-the-fly at request time?

[tim-salabim]

For reference https://sharp.pixelplumbing.com/api-output#tile

[mtennekes]

@edzer : is it possible to warp a rotated stars object? It doesn't work for a small object (see below). @tim-salabim : are curvilinear stars object supported in leafem::addStarsImage? It runs, but doesn't seem to work as expected (see below).

library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.2, PROJ 6.2.1
library(stars)
#> Loading required package: abind

# create regular stars object
m = matrix(1:20, nrow = 5, ncol = 4)
dim(m) = c(x = 5, y = 4) 
(s = st_as_stars(m))
#> stars object with 2 dimensions and 1 attribute
#> attribute(s):
#>       A1        
#>  Min.   : 1.00  
#>  1st Qu.: 5.75  
#>  Median :10.50  
#>  Mean   :10.50  
#>  3rd Qu.:15.25  
#>  Max.   :20.00  
#> dimension(s):
#>   from to offset delta refsys point values    
#> x    1  5      0     1     NA FALSE   NULL [x]
#> y    1  4      0     1     NA FALSE   NULL [y]
attr(s, "dimensions")[[2]]$delta = -1
plot(s)

# change rotation
s2 = s
attr(attr(s2, "dimensions"), "raster")$affine = c(0.1, 0.1)
plot(s2)

# set crs to 4326
(s3 = st_set_crs(s2, 4326))
#> stars object with 2 dimensions and 1 attribute
#> attribute(s):
#>       A1        
#>  Min.   : 1.00  
#>  1st Qu.: 5.75  
#>  Median :10.50  
#>  Mean   :10.50  
#>  3rd Qu.:15.25  
#>  Max.   :20.00  
#> dimension(s):
#>   from to offset delta                       refsys point values    
#> x    1  5      0     1 +proj=longlat +datum=WGS8... FALSE   NULL [x]
#> y    1  4      0    -1 +proj=longlat +datum=WGS8... FALSE   NULL [y]
#> sheared raster with parameters: 0.1 0.1

# warp to 3857
(s4 = st_warp(s3, crs = 3857)) # can you warp a rotated object?
#> Error: cannot allocate vector of size 1103.7 Gb

# transform to 3857
(s5 = st_transform(s3, crs = 3857))
#> stars object with 2 dimensions and 1 attribute
#> attribute(s):
#>       A1        
#>  Min.   : 1.00  
#>  1st Qu.: 5.75  
#>  Median :10.50  
#>  Mean   :10.50  
#>  3rd Qu.:15.25  
#>  Max.   :20.00  
#> dimension(s):
#>   from to offset delta                       refsys point
#> x    1  5     NA    NA +proj=merc +a=6378137 +b=... FALSE
#> y    1  4     NA    NA +proj=merc +a=6378137 +b=... FALSE
#>                       values    
#> x   [5x4] 61225.7,...,539900 [x]
#> y [5x4] -384285,...,-5565.98 [y]
#> curvilinear grid

plot(s5) # looks alright except that it seems cropped?

^{Created on 2020-01-27 by the reprex package (v0.3.0)}

library(leaflet)
library(leafem)

leaflet() %>% 
    addTiles() %>% 
    addStarsImage(s5) # why is it not rotated?

I could have opened issues in the stars and leafem repo; sorry for my laziness.

[tim-salabim]

@mtennekes currently we simply transform the value matrix of the layer to colors and then convert that to a raw, write to png and base64encode the png data. See

https://github.com/r-spatial/leafem/blob/master/R/addStarsImage.R#L107-L142

I pretty much copied that code from the leaflet::addRasterImage function.

If I am not mistaken, plot.stars converts a curvilinear grid to polygons and then plots those.

https://github.com/r-spatial/stars/blob/master/R/plot.R#L210

If we follow the same path, then the grid is rotated.

leaflet() %>% 
  addTiles() %>% 
  addPolygons(data = st_as_sf(s3)) # is now rotated!

[mtennekes]

Thanks! I have followed the same approach in the plot mode, so this should be easy. And I notice that it is not cropped anymore.

[mtennekes]

@tim-salabim Where shall implement this? It makes sense to do this in tmap/mapview, but then it would be good to have a warning or error when an user uses addStarsImage for a non-regular stars object.

[tim-salabim]

Thanks! I have followed the same approach in the plot mode, so this should be easy. And I notice that it is not cropped anymore.

I am just a bit concerned about size... Sure, these are simple polygons but there will likely be very any of them?

Where shall implement this? It makes sense to do this in tmap/mapview, but then it would be good to have a warning or error when an user uses addStarsImage for a non-regular stars object.

I think that makes sense. Let me add a warning to addStarsImage. Is there a dedicated function in stars to test this @edzer?

[tim-salabim]

@mtennekes we now have

leaflet() %>% 
   addTiles() %>% 
   addStarsImage(rev(s5), group = "strs1") %>%
   addStarsImage(s3, group = "strs2", project = TRUE) %>%
   addLayersControl(overlayGroups = c("strs1", "strs2"))

Warning messages:
1: cannot handle curvilinear or sheared stars images. Rendering regular gird. 
2: cannot handle curvilinear or sheared stars images. Rendering regular gird. 

[tim-salabim]

And I just fixed the typo... :-)

[edzer]

@tim there are no exported functions to figure out what kind of data cube we have, but I see that this makes sense. Proposal at https://github.com/r-spatial/stars/issues/248.

[mtennekes]

FYI: I placed the stars helper functions that I use in tmap in one file: https://github.com/mtennekes/tmap/blob/stars/R/stars_misc.R Many of them are copy pasted from stars, including functions like is_regular_grid and is_curvilinear. Not sure if it make sens to export all of them in stars, but at least for me it would be helpful.

Oh, and the last function in that script, cut_world_edges is a quick fix for the 4326 to 3857 warp problem (see https://github.com/r-spatial/mapview/issues/256#issuecomment-577337931).

[Nowosad]

I stared using the stars branch lately and found my first issue. The code below works fine with tmap < 3.0:

library(tmap)
#> Warning: replacing previous import 'sf::st_make_valid' by
#> 'lwgeom::st_make_valid' when loading 'tmap'
library(spDataLarge)
library(raster)
#> Loading required package: sp

tm_shape(nlcd) + 
  tm_raster(col = "levels")
#> Error: Invalid color specification. The available raster variables are: "layer".

^{Created on 2020-02-15 by the reprex package (v0.3.0)}

[mtennekes]

FYI, I rebranched tmap:

• I merged the 3.0/stars branch into master
• The old master (tmap 2.x) is renamed to the branch 'v2'.
• The same is done for tmaptools

[mtennekes]

I stared using the stars branch lately and found my first issue. The code below works fine with tmap < 3.0:

library(tmap)
#> Warning: replacing previous import 'sf::st_make_valid' by
#> 'lwgeom::st_make_valid' when loading 'tmap'
library(spDataLarge)
library(raster)
#> Loading required package: sp

tm_shape(nlcd) + 
  tm_raster(col = "levels")
#> Error: Invalid color specification. The available raster variables are: "layer".

Created on 2020-02-15 by the reprex package (v0.3.0)

Strange that it worked before, since levels is not a name of the raster layer (see names(nlcd)), but rather a generic column name in the factor attributes table.

[Nowosad]

I think it was implemented in relation to this issue - https://github.com/mtennekes/tmap/issues/204. I cannot find the commit though.

A related topic is that stars do not store raster levels. Small example from https://nassgeodata.gmu.edu/CropScape/ .

[edzer]

I'm not sure what you meant by "raster levels", @Nowosad , but stars now reads color tables through GDAL. stars can handle factor arrays properly, but I'm not sure whether GDAL has a mechanism for reading in / passing on legend entries (factor levels). Do you have an example file with colors and legend entries?

library(stars)
r = read_stars("CDL_2019_clip_20200301104844_1946824056.tif")
plot(r, key.pos = NULL)

[Nowosad]

NLCD data stores GDALRasterAttributeTable. You can see it, for example, in a small dataset for Puerto Rico - https://s3-us-west-2.amazonaws.com/mrlc/PR_landcover_wimperv_10-28-08_se5.zip.

[edzer]

Thanks! With

library(stars)
r = read_stars("pr_landcover_wimperv_10-28-08_se5.img", RAT = "Land Cover Class")
r = droplevels(r)
plot(r, key.pos = 4, key.width = lcm(5))

we now get

Do raster or terra read and plot such labels? or tmap or maview?

[tim-salabim]

For mapview it's a nope currently...

[Nowosad]

raster does read the labels:

library(raster)
#> Loading required package: sp
tf = tempfile(fileext = ".tif")
download.file("https://s3-us-west-2.amazonaws.com/mrlc/PR_landcover_wimperv_10-28-08_se5.zip", tf)
dir.create("tmp")
unzip(tf, exdir = "tmp")

x = raster("tmp/pr_landcover_wimperv_10-28-08_se5.img")
levels(x)[[1]][c(10:15), ]
#>    ID   COUNT Red Green Blue Opacity   Land.Cover.Class
#> 10  9       0   0     0    0     255                   
#> 11 10       0   0     0    0     255                   
#> 12 11 2699704  71   107  161     255         Open Water
#> 13 12       0 209   222  250     255 Perennial Snow/Ice
#> 14 13       0   0     0    0     255                   
#> 15 14       0   0     0    0     255
plot(x)

^{Created on 2020-03-17 by the reprex package (v0.3.0)}

[Nowosad]

@mtennekes Close?