Stars: Add support for subsetting stars proxy objects by non-rectangular regions

[kendonB]

This may still be on the to-do list but I couldn't find an issue related to it. It would be great to be able to subset a stars_proxy object by a set of 1D indices or (equivalently) an arbitrary set of points in the datacube. Is there a way to do this currently that does this in a minimal way in the backend?

I had a go at this for the raster package a while back:

https://github.com/rspatial/raster/blob/e2a8861f4612b7234f454f85faa0e86e1ad453a6/R/readCellsGDAL.R#L56-114

However, that patch is no longer actually exported. The approach I took breaks the IO stages into contiguous blocks in the 1D cell direction and iterates through those. I found that it had massive speed benefits when subsetting small parts of a raster as the previous version would iterate through and read entire rows of the raster then subset the rows. Is something like this (likely better than this) on the cards for stars?

[kendonB]

A nice interface might be to pass a data.frame of indices:

library(stars); library(tidyverse)
#> 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")
x1 = read_stars(tif, proxy = TRUE)
indices = tribble(
  ~x, ~y, ~band,
   1,  2,     5,
   1,  3,     2,
   1,  3,     3,
)
x1 = x1[indices]
x1 %>% st_as_stars()

[kendonB]

Or perhaps a data.frame as the second argument:

x1[, indices]

[edzer]

I can see that would come in handy, but I don't see how the output of such a selection can still be a data cube (rather than a data.frame).

[mdsumner]

(I didn't realise quite how that raster patch was working until you explained it here. The problem I saw in raster was its row-iteration worked against the native tiling within a raster - if it had one. ).

You seem to be asking two things, arbitrary cell access (and no longer cube as Edzer says) AND block wise access for efficiency. I guess the first is supported by the second if you classify the input query into native blocks?

[edzer]

Currently using the [,rows,cols,band] selector, you do get a sub-cube with gaps in the dimensions:

# example(read_stars)
> x1[,c(1,3),c(2,5),c(2,4)]
stars object with 3 dimensions and 1 attribute
attribute(s):
  L7_ETMs.tif   
 Min.   :45.00  
 1st Qu.:53.25  
 Median :62.00  
 Mean   :64.00  
 3rd Qu.:76.50  
 Max.   :84.00  
dimension(s):
     from to offset delta                       refsys point
x       1  2     NA    NA +proj=utm +zone=25 +south... FALSE
y       1  2     NA    NA +proj=utm +zone=25 +south... FALSE
band    1  2     NA    NA                           NA    NA
                                   values    
x        [288791,288819), [288848,288876) [x]
y    [9120718,9120690), [9120633,9120604) [y]
band                         [2,3), [4,5)    

but I guess that is not what you wanted, right?

[kendonB]

Currently using the [,rows,cols,band] selector, you do get a sub-cube with gaps in the dimensions

Did not know this but this is good to know. Does it still read these in in the subcubes or does it subset after reading in?

I don't think that one is quite the solution I was thinking though. My real-life example is about reading about 1000 .nc files of NZ stored in latlong. Even in the nice projections NZ is very diagonal but in latlong it is really diagonal:

https://ggplot2.tidyverse.org/reference/coord_map-1.png

So reading the whole cube reads in loads of NAs - about 6 times as many NAs as actual data. The IO stage is quite slow.

The request I was thinking of would, as you suggested, not extract a cube from the on-disk data. However, there's no reason why the in-memory object couldn't be a cube initially. Later you could have an option to make this a sparse representation of the array. You could initially allocate the cube in memory as NAs then fill in only the indices requested in the extraction.

The approach in the link I sent still fills in the data in a traditional 2D matrix in-memory, it just asks GDAL to only read the data from the indices requested rather than subsetting entire rows of the raster in memory after the IO stage. I'm not quite sure if you get big time savings for these diagonal extractions when you're pulling out most of the extent of the underlying raster.

I can see that would come in handy, but I don't see how the output of such a selection can still be a data cube (rather than a data.frame).

The output can be a data cube where the indices both in the bounds of the cube and not requested by the extraction are set to NA. Again, there might be good use cases for having an option to store these sparse like a data.frame.

You seem to be asking two things, arbitrary cell access (and no longer cube as Edzer says) AND block wise access for efficiency. I guess the first is supported by the second if you classify the input query into native blocks?

I don't have a lot of familiarity with how native tiling works in the common formats (.tif, .nc) but I presume it stores subregions of the object next to each other on disk. You could add an outer loop that loops through tiles in that case.

[edzer]

It's still unclear to me whether you are missing functionality or are looking for a particular access strategy.

[kendonB]

It is both. The user-facing functionality that's missing is extraction by cell index. The access strategy that's missing (and enables the performance gain for such user-facing functionality) is to extract those cell indices via a minimal number of contiguous blocks. For tiled data, the identification of the contiguous blocks would be done within each tile.

[edzer]

Closing, by lack of activity.