Per-pixel (statistical) calculations on a raster stack using GDAL

Question

In R, it's relatively trivial to perform per pixel calculations based on a raster stack (e.g., get std.dev for each pixel on a 12 layer GeoTIFF). Unfortunately, the speed is less than desirable when working with large rasters.

Recently, I've been trying to use GDAL for this purpose, but I can't seem to get the expression down. I've tried the following, but it doesn't seem to work as advertised (or how I think it should work):

gdal_calc.py -A tmean_monthly_1980.tif --calc="std(A)" --outfile=tstdev_1980.tif --allBands=A --calc="std(A)" --NoDataValue=-9999

In the above example tmean_monthly_1980.tif is a raster composed of 12 layers (Jan-Dec) of mean monthly temperatures in 1980 (stacked using gdal_merge). Although the GDAL command produces an output, the results are definitely not the standard deviation per pixel.

Another option might be to use single band files, but gdal_calc.py appears limited to "only" 26 inputs (A-Z; i.e., "-A file -B file -C file etc."), and the resulting command line/script would be fairly cumbersome when working with large, multi-layered rasters (e.g., time series).

Is there a specific syntax for performing per-pixel raster calculations on a multi-band image using GDAL? If only individual files/bands can be specified using gdal_calc.py, please provide a scripted example using Bash, Python, etc.

Specifically, my goal is to calculate the standard deviation and cv for each pixel in a multi-layer raster stack, or using individual files that make up the same raster stack.

Note: I'm holding out for a GDAL specific answer, but other FOSS solutions are welcomed and will be upvoted, especially if they are relatively compact and adaptable to large number of inputs.

Answer 1

Edit: adding some more stuff to make clearer, as per suggestion

#import the numpy and gdal libraries
import numpy as np
from osgeo import gdal

#an empty array/vector in which to store the different bands
layers = []

#open raster
ds = gdal.Open('raster.tif')

#loop thru bands of raster and append each band of data to 'layers'
#note that 'ReadAsArray()' returns a numpy array
for i in range(1, ds.RasterCount+1):
    layers.append(ds.GetRasterBand(i).ReadAsArray())

#dstack will take a number of n by m in tuple or list and stack them
#in the 3rd dimension so you end up with raster_stack being n by m by i, 
#where i is the number of bands
raster_stack = np.dstack(layers)

#call built in numpy functions std and mean, with a specified axis. if   
#no axis is set then it will return a number (scaler) but specifying
#axis=2 means it will calculate along the 'depth' axis, per pixel.
#with the return being n by m, the shape of each band.
std_raster = np.std(raster_stack, axis=2)
mean_raster = np.mean(raster_stack, axis=2)

Answer 2

You can use OTB for this. You can build a vrt to stack all your bands, then run OTB bandmath for the per pixel calculation.

otbcli_BandMath -il image_stack.vrt -out out.tif -exp "your expression"

Image file(s) are referenced in the expression in the order they appear after the -il switch (e.g., im1b1 refers to the first image im1 and first band b1. In the example below, im1b1 refers to band 1 of image_stack.vrt. If using a multi-band file, one just specifies the same image (im1), while using b1, b2, etc. to indicate the different bands.

To determine the cv with OTB, you begin by computing the mean:

otbcli_BandMath -il image_stack.vrt -out mean.tif -exp "avg(im1b1,im1b2,im1b3)"

You can then reuse the mean in order to compute the variance :

otbcli_BandMath -il image_stack.vrt mean.tif -out variance.tif -exp "( (im1b1-im2b1)^2+ (im1b2-im2b1)^2+(im1b3-im2b1)^2 )/2"

Once you have the mean and the variance, you can compute cv

 otbcli_BandMath -il variance.tif mean.tif -out cv.tif -exp "( (im1b1^0.5)/im2b1 )"

Another trick is to compute those values recursively (credit to J. Inglada from the OTB forum). You can of course adapt this bash expression for gdal_calc.py if you prefer.

RESULT_IMAGE = /tmp/tmp.tif
  │   NB_IMAGES = $(ls -l $INPUT_IMAGE_DIR/*.tif | wc -l)
  │   # Create an empty image with the appropriate size
  │   otbcli_BandMath -il $INPUT_IMAGE_DIR/image_name.tif -exp "im1b1*0" -out $RESULT_IMAGE
  │
  │   for IMAGE in $(find $INPUT_IMAGE_DIR/ -name *.tif)
  │   do
  │       otbcli_BandMath -il $IMAGE $RESULT_IMAGE -exp "im1b1+im2b1" -out $RESULT_IMAGE
  │   done
  │
  │   # Divide by the number of images
  │   otbcli_BandMath -il $RESULT_IMAGE -exp "im1b1/$NB_IMAGES" -out $RESULT_IMAGE

Note: you can use the same type of loop for variance computation based on the recursive variance calculation

Answer 3

As the question allows for other packages I'd like to propose a solution using RIOS (https://bitbucket.org/chchrsc/rios/). It is build on top of the GDAL Python bindings but provides a simpler interface, taking care of the actual raster I/O.

RIOS will provide the pixel values as a NumPy array, so you could use any of the inbuit stats functions or utilise other Python modules (e.g., pandas as suggested by Roland) for more advanced statistics. The data is read block at a time so it will work well for large rasters (which is why I've suggested it over the GDAL ReadAsArray function, which would work in a similar manner if you don't mind loading the entire dataset to memory).

An example Python script using RIOS is:

import sys
from rios import applier
from rios import cuiprogress
import numpy

def calcstats(info, inputs, outputs):
    # Calc standard deviation of bands for each pixel
    stdev_pixels = numpy.std(inputs.inimage, axis=0)
    # Reshape to make 3-dimensional array (required for output)
    outputs.outimage = stdev_pixels.reshape((1,stdev_pixels.shape[0],stdev_pixels.shape[1]))

# Set up input and output files
infiles = applier.FilenameAssociations()
infiles.inimage = 'tmean_monthly_1980.tif'

outfiles = applier.FilenameAssociations()
outfiles.outimage = 'tstdev_1980.tif'

# Set up options for output file
controls = applier.ApplierControls()
controls.setOutputDriverName("GTiff") # Output
controls.setCalcStats(True)           # Create stats and overviews for output file
controls.progress = cuiprogress.CUIProgressBar() # Show progress

applier.apply(calcstats, infiles, outfiles, controls=controls)