R package gapfill: How to convert R raster stack to 4-dimensional array and then back to raster stack?

Question

I'm currently analyzing MODIS EVI time series (MOD13Q1) for a small region between 2000 and 2015. Now, I would like to fill missing values using the new R package gapfill.

What I have done so far:

• Downloaded and pre-processed data using R MODIS package. This gave me 377 GeoTiff-files that are named MOD13Q1.A2000081.250m_16_days_EVI.tif and so on.
• Used the QA layers from MOD13Q1 to remove pixels with clouds.

As a result, I now have 377 GeoTiffs with NA values for clouds:

> stack(evi_without_clouds_file_paths)
class       : RasterStack 
dimensions  : 121, 122, 14762, 377  (nrow, ncol, ncell, nlayers)
resolution  : 257.2911, 257.3028  (x, y)
extent      : 448840.3, 480229.9, 7766995, 7798129  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=38 +south +datum=WGS84 +units=m +no_defs     +ellps=WGS84 +towgs84=0,0,0 
names       : MOD13Q1.A//6_days_EVI, MOD13Q1.A//6_days_EVI, MOD13Q1.A//6_days_EVI, MOD13Q1.A//6_days_EVI, MOD13Q1.A//6_days_EVI, MOD13Q1.A//6_days_EVI, MOD13Q1.A//6_days_EVI, MOD13Q1.A//6_days_EVI, MOD13Q1.A//6_days_EVI, MOD13Q1.A//6_days_EVI, MOD13Q1.A//6_days_EVI, MOD13Q1.A//6_days_EVI, MOD13Q1.A//6_days_EVI, MOD13Q1.A//6_days_EVI, MOD13Q1.A//6_days_EVI, ... 
min values  :                 -1195,                  -846,                  -220,                  -185,                  -321,                  -344,                  -128,                  -159,                  -309,                    99,                  -281,                  -268,                  -200,                   218,                  -290, ... 
max values  :                  9205,                  9902,                  8466,                  7933,                  7696,                  6700,                  7019,                  6825,                  5578,                  5163,                  4934,                  5028,                  4561,                  5691,                  4262, ... 

However, the input data for gapfill has to be in the following format: "Numeric array with four dimensions. [...]. the data should have the dimensions: x coordinate, y coordinate, seasonal index (e.g., day of the year), and year."

I'm familiar with the raster package, but I have no experience whatsoever working with multidimensional arrays. As a result, I spent the whole day trying to convert my data, but had no success so far. So my questions would be:

1. How to convert raster stacks to 4-dimensional named arrays?
2. How to convert 4-dimensional named arrays back to either a brick with 377 layers or 377 single GeoTiff files after having gaps filled with gapfill?

I apologize for not being able to provide sample data. I tried to simulate similar data, but didn't manage to do so.

I know that it is difficult to answer the question without sample data, but I would also be very happy about some pointers to the right direction.

Here is the closest I have gotten so far:

require(tidyr)
evi.stack <- stack(evi_without_clouds_file_paths)
x <- as.data.frame(evi.stack, xy = TRUE)
x <- gather(x, key = "scene", "value", -x, -y)
x$scene <- extract.id(x$scene)
x$day.of.year <- as.numeric(str_sub(x$scene, 6, 8))
x$year <- str_sub(x$scene, 2, 5)

evi.dim.long <- unique(x$x)
evi.dim.lat <- unique(x$y)
evi.dim.day.of.year <- unique(x$day.of.year)
evi.dim.year <- unique(x$year)

x.vec <- x$value
dim(x.vec) <- c(length(evi.dim.long),
                length(evi.dim.lat),
                length(evi.dim.day.of.year),
                length(evi.dim.year))
dimnames(x.vec) <- list(evi.dim.long,
                        evi.dim.lat,
                        evi.dim.day.of.year,
                        evi.dim.year)

If I only import data from a single year, this seems to kind of work (although the resulting image was "upside down"). I was able to fill gaps with Gapfill, but then had no idea how to convert data back to a format that could be further processed with the raster package. If I try to import data from more than one year, the dim command fails.

Answer 1

See the example below:

library("gapfill")
library("raster")
## create raster data for demo
## extent: x=10, y=10, month=3, years=2
data <- array(runif(600), c(10,10,6))
data[c(1,5,54,76,150,450,556)] <- NA
input_stack <- stack(brick(data))
plot(input_stack)

## create array and predict missing values with gapfill()
tmp <- array(input_stack, dim=c(10,10,3,2))
input_array <- aperm(tmp, c(2,1,3,4))
Image(input_array, col=rev(terrain.colors(100)), asRaster=TRUE)
output <- Gapfill(data=input_array)
output_array <- output$fill
Image(output_array, col=rev(terrain.colors(100)))

## convert back
output_stack <- stack(brick(array(output_array, c(10,10,6))))
plot(output_stack)

If the dataset is too big to be loaded into the memory, one strategy is to divide the dataset into smaller parts, gapfill them, and join the gapfilled parts afterwards. Gapfill() predicts each missing value based on a subset of the data. Depending on the chosen subset strategy different divisions of the dataset make sense. For the default subset strategy you could, e.g., divide the data into smaller spatial regions. Make sure that the smaller parts have enough overlap.

For more information see http://dx.doi.org/10.1109/TGRS.2017.2785240, http://user.math.uzh.ch/furrer/download/tgrs/supplementary_material.pdf, the help pages: ?Gapfill; ?Extend; ?Subset-Predict.