In the source for gdal_calc.py, the calculation is made directly using eval:
myResult = eval(opts.calc, global_namespace, local_namespace)
That would suggest that any well-formed expression that also evaluates on the command line will work. According to the documentation, you may use gdalnumeric syntax with +-/*, and/or numpy functions. You can test your functions using small dummy arrays in the interactive shell, then use the same calls in gdal_calc.
Keep in mind that chaining together multiple numpy functions is likely to produce temporary in-memory arrays that can substantially increase memory usage, especially when dealing with large images.
You can look at the numpy documentation for a list of all the functions: routines. The ones you are after are likely here: math or here: routines.logic.
This is where functions like minimum are coming from, it's just that the namespace is already imported. Really, it's numpy.minimum, etc
Here is an approach that imputes NA values based on a local polynomial regression (loess). The default smoothing parameter (s) is fairly stable but, if the data is stochastic (say, climate data), is also something that you may want to test.
I would also point out that, in the temporal analysis literature, it is common practice to smooth a time-series using a loess regression. I added an argument "smooth.data" that allows for the function to return not only imputed NA values but also the smoothed time-series. If smooth.data = FALSE only the NA values are replaced else all of the values are replace with the smoothed loess estimate. If the data is smoothed attention needs to be paid to the smoothing parameter (s) to ensure that the data is being adequately smoothed without also changing the properties of the original distribution. In a recent monthly 20-year LAI trend analysis I found that the default s=0.80 was adequate for imputation of NA values however, was too aggressive in smoothing the entire time-series. I had to specify s=0.30 to achieve the desired smoothing.
#### Loess impute function
# y Vector with NA values or to be smoothed
# x.length length of resulting vector
# s Smoothing parameter (see loess span argument)
# smooth.data (FALSE/TRUE) smooth all of the data
impute.loess <- function(y, x.length = NULL, s = 0.80,
smooth.data = FALSE, ...) {
if(is.null(x.length)) { x.length = length(y) }
options(warn=-1)
x <- 1:x.length
p <- loess(y ~ x, span = s, data.frame(x=x, y=y))
if(smooth.data == TRUE) {
y <- predict(p, x)
} else {
na.idx <- which( is.na(y) )
if( length(na.idx) > 1 ) {
y[na.idx] <- predict(p, data.frame(x=na.idx))
}
}
return(y)
}
Here we apply the impute.loess function. I have had numerous problems using calc to apply this function so, I am using a for loop to iterate through the rows of the raster stack with getValues. An alternative approach is to coerce into a SpatialPixelsDataFrame and operate on the @data slot using apply, the same as below. This is more direct and does not require the back-and-forth between raster and sp. In this example, I make a copy of the raster stack and populate it with NA values. To be memory safe you could use writeStart and writeValues to create and write to a raster on disk.
ndvi.new <- ndvi
ndvi.new[] <- NA
for( rl in 1:nrow(ndvi) ) {
v <- getValues(ndvi, rl, 1)
ndvi.new[rl,] <- as.matrix(t(apply(v, MARGIN=1, FUN=impute.loess)))
}
Since apply can take additional function arguments you can change the s and smooth.data agruments by supplying them after the loess.smooth call in apply eg.,
as.matrix(t(apply(v, MARGIN=1, FUN=impute.loess, smooth.data = TRUE, s = 0.20)))
Occasionally, this function will miss NA values on the absolute tails of the time-series because it is impossible to converge on an end value. This is a statistical limitation and can only really be solved by taking the next/previous value and assigning to the missing data. You could add an if/else statement in the function to account for this or just deal with it in your analysis.
Best Answer
you must use maximum instead of max