Indexing in more than 2 dimensions in R?

[marciolegal]

Hi,
I'm looking to use Tile DB as a replacement for NetCDF to store tridimensional (row, col, time) data in R. I'm currently facing three issues in R.

The first issue is trying to write data from a tileDB array using 3-dimensional R arrays.
For example, I changed Example 1 to include an extra dimension "time", using the same attributes:

dom <- tiledb_domain(dims = list(
  tiledb_dim("rows", c(1L, 10L), 10L, "INT32"),
  tiledb_dim("cols", c(1L, 5L), 5L, "INT32"),
  tiledb_dim("time", c(1L, 2L), 2L, "INT32")))
schema <- tiledb_array_schema(dom, attrs = c(tiledb_attr("a", type = "INT32"),
                                             tiledb_attr("b", type = "FLOAT64"),
                                             tiledb_attr("c", type = "CHAR", ncells = NA_integer_)))

I can't just pass a list of tridimentional R arrays to the tileDB array. It only recognizes the first 2 dimensions of the R array:

data <- list(a=array(seq(1:100), dim = c(10,5, 2)),
             b=array(as.double(seq(101,by=0.5,length=100)), dim = c(10,5,2)),
             c=array(rep(c(letters[1:26], "brs", "asdf", LETTERS[1:22]), 2), dim = c(10,5,2)))
A <- tiledb_array(uri = uri, as.data.frame = TRUE, query_layout = "ROW_MAJOR")
A[] <- data

The above throws:

Error in libtiledb_query_submit(qryptr) :
[TileDB::Writer] Error: Buffer sizes check failed; Invalid number of cells given for attribute 'c' (50 != 100)

If I instead just pass a list of two-dimensional arrays with 100 values, it works.

Is it not currently possible to feed the n-dimentional tileDB array with a n-dimensional R array?

The second issue is that I also can't access the values using tridimensional indexing. After feeding my A array with a list of 2-dimensional R arrays as said above, A[1:2, 1:2, 1] gives:

"rows" "cols" "time" "a" "b" "c"
1 1 1 1 101 "a"
1 1 2 51 126 "a"
1 2 1 11 106 "k"
1 2 2 61 131 "k"
2 1 1 2 101.5 "b"
2 1 2 52 126.5 "b"
2 2 1 12 106.5 "l"
2 2 2 62 131.5 "l"

When I expected to see only data from time 1. Instead it is returning all times. In the R docs it says the third argument is ignored. So how do we query multidimensional arrays in R?

The third problem arises when trying to write the data using the query:

qry <- tiledb_query_set_buffer(qry, "c", data[["c"]])

throws:

Error in libtiledb_query_set_buffer(query@ptr, attr, buffer) :
Invalid attribute buffer type for attribute 'c': STRSXP

It seems that character arrays are not supported in this way. Works fine when using A[]<-data

[eddelbuettel]

Thanks for the bug report. I transferred it from the repository for TileDB Embedded, our underlying C++ library, to the repository for the R package.

I will have to take a closer look at the dense case, and three dimensions.

One issue is that because R imposes some constraints on the [ operator (namely the signature with i, j and three dots), I do not think we can generally support triple indexing directly in the operator. I opted to generalize it with a general accessor function selected_ranges() which can receive a (named) list with entries for each dimension (or a subset) where each list element is a two-column matrix where each row contains a pair of 'from' and 'to' selections. (So to select element 1:10 and 100:110, we would pass two rows with value 1 and 10, and 100 and 110, respectively.)

Below is a small variation on your example where I use a sparse array, created from a data frame which contains your input data arrays 'flattened' into vectors forming a data frame. I used sampling for the row, col and time indices here; we could equally use sequences. One advantage of using a sparse array is that we can then mix index type -- your time column could actually be an R Date or Datetime object (as well as one with a different time resolution as TileDB supports more date types).

library(tiledb)
set.seed(123)  # make the RNG draws reproducible
moddata <- data.frame(rows=sample(10, 100, TRUE),
                      cols=sample(5, 100, TRUE),
                      time=sample(2, 100, TRUE),
                      a=seq(1:100),
                      b=as.double(seq(101,by=0.5,length=100)),
                      c=rep(c(letters[1:26], "brs", "asdf", LETTERS[1:22]), 2))
uri <- tempfile()
fromDataFrame(moddata, uri, sparse=TRUE, col_index=c("rows", "cols", "time"))

## read back
arr <- tiledb_array(uri, as.data.frame=TRUE, extended=TRUE) 
selected_ranges(arr) <- list(rows=cbind(1,2),       # rows 1 to 2
                             cols=cbind(1,2),       # cols 1 to 2
                             time=cbind(1,1))       # time equal to 1 
res <- arr[]
res

This generates (with the given seed) a four-row selection:

> library(tiledb)
> set.seed(123)  # make the RNG draws reproducible
> moddata <- data.frame(rows=sample(10, 100, TRUE),
+ cols=sample(5, 100, TRUE),
+ time=sample(2, 100, TRUE),
+ a=seq(1:100),
+ b=as.double(seq(101,by=0.5,length=100)),
+ c=rep(c(letters[1:26], "brs", "asdf", LETTERS[1:22]), 2))
> uri <- tempfile()
> fromDataFrame(moddata, uri, sparse=TRUE, col_index=c("rows", "cols", "time"))
> arr <- tiledb_array(uri, as.data.frame=TRUE, extended=TRUE)
> selected_ranges(arr) <- list(rows=cbind(1,2), 		# rows 1 to 2
+ cols=cbind(1,2),       # cols 1 to 2
+ time=cbind(1,1))       # time equal to 1
> res <- arr[]
> res
  rows cols time  a     b c
1    1    2    1 42 121.5 N
2    1    2    1 55 128.0 e
3    1    2    1 58 129.5 h
4    2    2    1 84 142.5 F
> 

If we, say, also selected rows 7 to 8, we get twelve rows back:

> selected_ranges(arr) <- list(rows=cbind(c(1,7), c(2,8)),
+ cols=cbind(1,2),       # cols 1 to 2
+ time=cbind(1,1))       # time equal to 1
> res <- arr[]
> res
   rows cols time  a     b c
1     1    2    1 42 121.5 N
2     1    2    1 55 128.0 e
3     1    2    1 58 129.5 h
4     2    2    1 84 142.5 F
5     7    1    1 34 117.5 F
6     7    1    1 54 127.5 d
7     7    1    1 74 137.5 x
8     7    1    1 19 110.0 s
9     7    2    1 98 149.5 T
10    8    1    1 71 136.0 u
11    8    2    1 51 126.0 a
12    8    2    1 95 148.0 Q
> 

[marciolegal]

Thank you for your response and code. selected_ranges works with my example so it's fine with dense arrays too. Quick question: does a tiedb array support parallel writes?

[eddelbuettel]

Great -- and yes, parallel reads and writes are supported on all filesystems that support it.

(I took the liberty of editing your reply down as it was shown at GitHub will all lines from the email. You could cut the quoted lines when you reply by email.)

[marciolegal]

Thanks.

I noticed that unfortunately select_ranges doesn't seem to work for writing values though.
For example, I can use: `

A <- tiledb_array(uri = uri, as.data.frame = TRUE, 
                  selected_ranges = list(rows = cbind(1,2), cols = cbind(1,2), time = cbind(1,1)),
                  query_type = "WRITE")

which has the same affect as creting the array and then change the ranges right?

However, if I try to set the values for these indices, it asks for the entire set:

f = data.frame(a = c(1:4), b = c(1:4), c = LETTERS[1:4])
A[] <- f

this throws:

Error in libtiledb_query_submit(qryptr) :
[TileDB::Writer] Error: Buffer sizes check failed; Invalid number of cells given for attribute 'c' (4 != 100)

but:

A[]
rows cols time a b c
1 1 1 1 NA NaN €
2 1 2 1 NA NaN €
3 2 1 1 NA NaN €
4 2 2 1 NA NaN €

So my hope was that I would only write these four values of a-c,

[eddelbuettel]

Your schema has six columns, but your data.frame f only has three. When you use a sparse array, you must supply its index columns, here rows, cols and time, as well. If you do that you can write incrementally.

One caveat is that the original domain and extent values may limit what you write later, and one solution is to explicitly create the sparse array with 'wider' limits via the tile_domain. Here is a short example from one of the unit tests where I write out the penguins data set (whicg has three known 'species' (NB: we have no domain restriction on character columns) and observed years 2007 to 2009. So I set year explicitly to something wider and can then, as a test, add two more random values:

library(palmerpenguins)
uri <- tempfile()
fromDataFrame(penguins, uri, sparse = TRUE,
              col_index = c("species", "year"),
              tile_domain=list(year=c(1966L, 2021L)))
arr <- tiledb_array(uri, as.data.frame=TRUE)
## new data
newdf <- penguins[1:2,]                           # quicker than assigning all columns...
newdf$species <- c("Fred", "Ginger")
newdf$island <- c("Manhattan", "Staten Island")
newdf$year <- c(1966L, 1969L)                     # int is important
arr[] <- newdf

[eddelbuettel]

Hi @marciolegal -- your original issue of writing the three-dimensional dense array can be addressed as well (as we have done in this branch which should get merged soon as well).

However, given that arrays with more than two dimensions are somewhat infrequently used, we currently do not have a convenient way to retrieve a arrays with dimensions of, say, c(10, 5, 2) as here. Would it help your use case if we created such an option? tiledb_array can be asked to return a data.frame, or a matrix; we could generalize this to a dimension setting. As the size of returned objects may vary, maybe that is too restrictive (though maybe a count of unique elements per dimension may help) ?

[marciolegal]

Hi Dirk Yes that would be great for my application. I believe spatio-temporal indexing would be interesting in general for the remote sensing community using R. See for example the gdal_cubes and stars project in R. They currently write using gdal and it appears the preferred format for storage is now netcDF for spatio-temporal arrays, which doesn't allow for parallel writes and is harder to scale compared to tileDB.

[eddelbuettel]

The correction (and additional test) from the branch I mentioned were merged to the mainline a little earlier.

[eddelbuettel]

@marciolegal -- yes, it would relatively straightforward to do something like the following:

> uri <- "someArray"        # adjust as needed, written with code above and package updated after #260
> library(tiledb)                                   
> arr <- tiledb_array(uri)   
> domain(schema(arr))       # we can access the domain of the schema programmatically too
> d <- c(10, 5, 2)          # this could be a function of the retrieved data too
> res <- arr[]              # retrieve data, here all
> array(res$a, dim=d)
, , 1

      [,1] [,2] [,3] [,4] [,5]
 [1,]    1   11   21   31   41
 [2,]    2   12   22   32   42
 [3,]    3   13   23   33   43
 [4,]    4   14   24   34   44
 [5,]    5   15   25   35   45
 [6,]    6   16   26   36   46
 [7,]    7   17   27   37   47
 [8,]    8   18   28   38   48
 [9,]    9   19   29   39   49
[10,]   10   20   30   40   50

, , 2

      [,1] [,2] [,3] [,4] [,5]
 [1,]   51   61   71   81   91
 [2,]   52   62   72   82   92
 [3,]   53   63   73   83   93
 [4,]   54   64   74   84   94
 [5,]   55   65   75   85   95
 [6,]   56   66   76   86   96
 [7,]   57   67   77   87   97
 [8,]   58   68   78   88   98
 [9,]   59   69   79   89   99
[10,]   60   70   80   90  100

and we could of course place the array creation inside the res <- arr[] accessor as we do for matix and data.frame.

Otherwise, would you have a suitable example data set used with with gdalcubes and/or stars? We do have a number of examples in Python, often relying on gdal binaries that are TileDB-aware, to create disk arrays from tiff or other files. We could work on R importers, but it might be helpful to tie it to a particular data set example if you have one.

[marciolegal]

Looks good. Do you think something like that could be used in the [ and [<- implementations in TileDB?

Here is a netCDF file with some variables covering Sao Felix - BR Municipality and some code using the stars package:

s = read_stars("testencX.nc", proxy = TRUE) # proxy avoids loading everything into memory. used here only to illustrate
r = s[1:2, 1:1, 1:1, 1:3] %>% as.data.frame 

the quadruple index is for variable, row, column and time. Stars proxy objects with multiple attributes can be converted to data frames directly, but I could get a single attribute by converting the result to a stars object (memory) then refering the variable with $:

r = s[1:2, 1:2, 1:2, 1:3] %>% st_as_stars %>% .$forestcover

, , 1

     [,1] [,2]
[1,]    1    1
[2,]    1    1

, , 2

     [,1] [,2]
[1,]    1    1
[2,]    1    1

, , 3

     [,1] [,2]
[1,]    1    1
[2,]    1    1

However, this is much slower than using the ncdf4 package to extract the same data: 1.14s in stars against 0.0s in ncdf4 in my laptop.

[eddelbuettel]

Looks good. Do you think something like that could be used in the [ and [<- implementations in TileDB?

The [ accessor could probably be extended to return arrays as sketched above if the index information is complete (as in our example here).

The [<- accessor for writes cannot, I think, write that way for dense arrays created this way. We can write using subarray settings and suitable tiles and extents. But in general you will have more flexibility with sparse arrays.

Thank you a lot for the nc file which I was able to download. I can run both your examples and will take a closer look. For me, running system.time({s = read_stars("testencX.nc", proxy = TRUE) ; as.data.frame( s[1:2, 1:1, 1:1, 1:3] ) }) takes only about 0.4 seconds. Is it that step that takes 1.14s on your machine?

TileDB could be an alternative for you to the current NetCDF format. You would have to re-ingest the data, we have for example Docker containers that can help read from other (geospatial) formats.

[marciolegal]

Good to know you you were able to run the example.
the following took 1.14s in my 1.8GHZ intel processor:
r = s[1:2, 1:2, 1:2, 1:3] %>% st_as_stars %>% .$forestcover
there might be some overhead there somewhere

Data ingestion is ok, we are being able to do it using gdal to read geotiffs into R then writing it to netcdf (faster than using gdal virtual raster to read/write data afterwards). I'll try to use now tileDB with our big files (~70GB) and check speeds.

Do you think sparse arrays reading / writing speeds would be comparable to netcdf's? Because then we could use multiple cores to go faster - that's the idea.

[eddelbuettel]

Quicker here for me (on an i7 desktop with decent SSD disks):

> system.time({ r = s[1:2, 1:2, 1:2, 1:3] %>% st_as_stars %>% .$forestcover })
   user  system elapsed 
  0.294   0.021   0.314 
There were 19 warnings (use warnings() to see them)
> 

Going via gdal is what we do too -- see for example this Docker container which has TileDB backings. I would expect TileDB to be competitive in terms of timing and read/write performance. There are always different parameters to tune and we would be happy to help further once you have a first example.