7

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.

  • It might be worth looking closer to find a faster way in R/raster, there are a few ways to do it and the details matter. Happy to do that perhaps outside this q if you want – mdsumner Mar 8 '14 at 6:19
  • Sure, I'm always up for leveling up in R when it comes to raster, especially when it comes to processing CONUS scale stacks of rasters! – Prophet60091 Mar 8 '14 at 12:57
  • 3
    Another approach might be to use python/numpy to do this. Could throw in the pandas module if you want to do something more complicated in terms of stats. – Roland Mar 10 '14 at 15:25
  • you could also have a look at OTB – radouxju Aug 19 '14 at 12:43
  • 1
    I second python+numpy. Once you have your 'raster stack' which is a numpy array n by m by i (i being the number of rasters). A quick numpy.std or numpy.mean will give some stats. But be sure to apply it to the correct axis (2 for this problem, I think) so it returns a n by m result. You can use numpy to load .asc text raster really easily (ignore headers and parse by spaces). – user1269942 Sep 22 '14 at 4:58
7

pardon any syntax errors...hope this helps. [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)
  • +1 This is a really good solution. Please consider adding a few lines of text describing what is going on. – Aaron Sep 22 '14 at 6:17
4

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

  • The script is incomplete. It only calculates the average pixel value of all input rasters (treated as individual rasters vs. a raster stack). – Prophet60091 Sep 3 '14 at 18:29
  • the script is for a list : the equation with a raster stack is given in my second code snippets, which can be extended. I've never had problem with the lenght of expression with otbcli_bandmath. If you really want an adaptive per pixel calculation with OTB, you need to create your own filter. You can do a lot more with the library than with the applications). I'll propose to the OTB team to add this in the next release (you need to read an otbVectorImage as input, then you write a filter that process each pixel as a vector). – radouxju Sep 3 '14 at 19:34
1

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)
  • I just tried our rios-1.3.1, but couldn't get it to work. The script runs without any (obvious) error, but there's no output. testrios.py passes all tests. – Prophet60091 Aug 30 '14 at 2:41
  • That's weird - does it print status about percent completion. I've uploaded a copy of the script here: gist.github.com/danclewley/b769b81947ea24807a04 if you want to test this incase it's an indentation problem – danclewley Sep 3 '14 at 6:01
  • Oops. Looks like the issue was that I installed riosRoot into the original untarred rios-1.3.1 folder (e.g., /home/mine/rios-1.3.1), instead of a separate directory (e.g. /home/mine/rios). Rerunning got me: "Warning 6: Driver GTiff does not support COMPRESSED creation option Warning 6: Driver GTiff does not support IGNOREUTM creation option". It did work though, and it was pretty fast! Definite +1 (i'd +2 if possible) – Prophet60091 Sep 3 '14 at 14:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.