Array transformation: RasterStacks of time-series into Array of matrices

I have 9 RasterStacks with 65 layers containing time-series of satellite imagery. Each stack contains a time-serie of each Sentinel-2 band (B2, B3, B4, B5, B6, B7, B8, B11 and B12). The shape of the stacks is 256 x 256 x 65 (X pixels, Y pixels, Time). I need to create a 256 x 256 array containing a matrix, where the Y axis of the matrix is each band, and the X axis is the pixel values in each time. This looks like this:

Any help or code example where to approach this? I need help to start developing this transformation.

• A matrix can only be 2-dimensional - anything 3 or more dimensional in R is an "array" and can be constructed with `array()`. I'm not certain what you want to end up with but is it a 256x256x585 sized 3-d array? Where 585 is 65 (layers) times 9 (stacks)? Do you care about the ordering of time within stacks in the third dimension? Or do you want a four-dimensional 256x256x9x65 array? The elements of a matrix (or an array) can only be scalar values. Commented Aug 23, 2022 at 13:51
• I need a 4d array of (256x256x9x65). To better understand why I need this. I would like to train a Convolutional Neural Network with each of these 9x65 matrices. This is my approach to applying a CNN to a time series to generate a pixel-based classification. Commented Aug 23, 2022 at 15:02

I've been working on an array library for Scheme: https://srfi.schemers.org/srfi-231/srfi-231.html

In this library, the code would be a one-liner:

``````(array-curry (array-stack 2 (list B2 B3 B4 B5 B6 B7 B8 B11 B12)) 2)
``````

Many of the routines in my library are motivated by similar routines in Python, etc. `array-stack` is very much like `numpy.stack`: https://numpy.org/doc/stable/reference/generated/numpy.stack.html which you can use here similarly.

Once you have your 256 x 256 x 9 x 65 array in Python, there's slicing notation to choose subarrays, e.g., `a[i,j,:,:]` should give you the 9 x 65 array associated with the `[i,j]` pixel.

• Thank you very much for your effort. I'm going to check your library. Commented Nov 3, 2022 at 12:03

@Spacedman, for the moment I was able to do this using R. Below is a simple example with 4 RasterStack of 5x5 pixels and 12 dates. So in this example, my goal is to create an array of (5x5x4x12).

``````### Create toy data:

# First we create matrices containing the values of the time-series.
# Each matrix corresponds to a layer in the RasterStack containing the time-series.
toyMatrix1 <- matrix(c(1:25), nrow = 5, ncol =5)
toyMatrix2 <- toyMatrix1 + 1; toyMatrix3 <- toyMatrix1 + 2
toyMatrix4 <- toyMatrix1 + 3; toyMatrix5 <- toyMatrix1 + 4
toyMatrix6 <- toyMatrix1 + 5; toyMatrix7 <- toyMatrix1 + 6
toyMatrix8 <- toyMatrix1 + 7; toyMatrix9 <- toyMatrix1 + 8
toyMatrix10 <- toyMatrix1 + 9; toyMatrix11 <- toyMatrix1 + 10
toyMatrix12 <- toyMatrix1 + 11

# We stack the matrices to create a 3D array.
toyArray <- simplify2array(list(toyMatrix1, toyMatrix2, toyMatrix3, toyMatrix4, toyMatrix5, toyMatrix6,
toyMatrix7, toyMatrix8, toyMatrix9, toyMatrix10, toyMatrix11, toyMatrix12))

# Then we create a RasterStack object.
toyRaster <- raster::stack(raster::brick(toyArray))
plot(toyRaster)

# Now we create the time-series of the Sentinel-2 bands. We are going to create only the 2,3,4 and 8 toy bands
# and we are going to assume that all RasterStacks are equals.

tsB02 <- tsB03 <- tsB04 <- tsB08 <- toyRaster
``````

At this point we have a dataset comparable to what I need, but simpler. My approach is:

``````# We convert the RasterStack into 3D arrays
a02 <- as.array(tsB02) #! We use this first Array as reference for the loops
a03 <- as.array(tsB03)
a04 <- as.array(tsB04)
a08 <- as.array(tsB08)

# All the arrays should have the same shape: dim(a02) == dim(a003) == dim(a04) == dim(a08)
dim(a02) == dim(a03); dim(a02) == dim(a04); dim(a02) == dim(a08)

matVec <- list() # Empty list to fill with the matrices (Bands x Time)

c = 0

for(i in 1:dim(a02)[1]){
for(j in 1:dim(a02)[2]){
c = c+1
# Create empty vectors to fill
v_b02 <- list()
v_b03 <- list()
v_b04 <- list()
v_b08 <- list()
for(t in 1:dim(a02)[3]){
# Looping pixel values in the -t- dimension
pix02 <- a02[i,j,t]; pix03 <- a03[i,j,t]; pix04 <- a04[i,j,t];pix08 <- a08[i,j,t]
# Append pixel values in the -t- dimension to create temporal vectors
v_b02[[t]] <- pix02
v_b03[[t]] <- pix03
v_b04[[t]] <- pix04
v_b08[[t]] <- pix08
}
# Create matrix for the (i,j) pixel
mat <- matrix(c(unlist(v_b02),
unlist(v_b03),
unlist(v_b04),
unlist(v_b08)),
nrow=length(v_b02), ncol=4)
mat <- t(mat)

# Append (i, j) pixel matrices to matVec
matVec[[c]] <- array(mat, dim=c(nrow(mat), ncol(mat), 1))
#matVec <- append(matVec, array(mat, dim=c(nrow(mat), ncol(mat), 1)))
}
}
``````

Here the matVec list contains my (4x12) arrays. Now I need to reshape this elements into a (5x5x4x12) array. I tryed this:

``````finalMatrix <- matrix(matVec, ncol = ncol(a02), nrow = nrow(a02))

> finalMatrix
[,1]       [,2]       [,3]       [,4]       [,5]
[1,] numeric,48 numeric,48 numeric,48 numeric,48 numeric,48
[2,] numeric,48 numeric,48 numeric,48 numeric,48 numeric,48
[3,] numeric,48 numeric,48 numeric,48 numeric,48 numeric,48
[4,] numeric,48 numeric,48 numeric,48 numeric,48 numeric,48
[5,] numeric,48 numeric,48 numeric,48 numeric,48 numeric,48
``````

However each element of `finalMatrix` contains a list of one element instead of an array.

``````> finalMatrix[1,1]
[[1]]
, , 1

[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
[1,]    1    2    3    4    5    6    7    8    9    10    11    12
[2,]    1    2    3    4    5    6    7    8    9    10    11    12
[3,]    1    2    3    4    5    6    7    8    9    10    11    12
[4,]    1    2    3    4    5    6    7    8    9    10    11    12
``````