How to center continuous raster layer by subtracting mean and dividing by standard deviation using QGIS raster calculator?

Question

I'm trying to extrapolate a spatial model to a map in QGIS using the raster calculator. The model predicts the value of a pixel given the environmental conditions at this location (e.g. elevation, slope, ect.). Hence, each explanatory variable is a raster layer. For the model I centered my covariates by subtracting the mean and dividing by the standard deviation. Now I want to do the same to my raster layers. As trivial as this seems, I didn't manage to do so. I tried (among similar attempts - like changing upper case, etc.):

(my.raster-mean(my.raster))/sd(my.raster)

The calculator tells me the expression is valid, however, mean() and sd() don't work.

And I could not find a site with the built-in functions for the raster calculator.

Answer 1

Raster calculator is a wrong tool to do that because it is for map algebra. The mean(my.raster) or sd(my.raster) expressions return my_raster (not the constants that you are hoping). For this reason, your expression (my.raster-mean(my.raster))/sd(my.raster) will be evaluated to 0 for each pixel.

At the below image you can see that any "function" it will be "valid" and it will return the same raster.

As it was pointed out by Michael, you need to evaluate mean and standard deviation constants first by other methods. You can do that with GDAL-python, PyQGIS or R (with rgdal o raster libraries).

Next, one example with PyQGIS:

For the complete process, you need to adapt this solution:

Qgis raster calculator not giving output

Editing Note 1:

The complete solution is this:

from qgis.analysis import QgsRasterCalculator, QgsRasterCalculatorEntry

layer = iface.activeLayer()
entries = []

# Define band1
band1 = QgsRasterCalculatorEntry()
band1.ref = 'band@1'
band1.raster = layer
band1.bandNumber = 1
entries.append( band1 )

renderer = layer.renderer()
provider = layer.dataProvider()
extent = layer.extent()

stats = provider.bandStatistics(1, QgsRasterBandStats.All,extent, 0)

myMean = stats.mean
myStdDev = stats.stdDev

myFormula = '(band@1 -' + str(myMean) +')/' + str(myStdDev)

print "mean = ", myMean 

print "stdev = ", myStdDev

# Process calculation with input extent and resolution
calc = QgsRasterCalculator( myFormula, 
                            '/home/zeito/pyqgis_data/outputfile.tif', 
                            'GTiff',
                            layer.extent(), 
                            layer.width(), 
                            layer.height(), 
                            entries )

calc.processCalculation()

It works nicely.

Editing Note 2:

This is the complete code with R (Benedikt Gehr approach and RobertH approach). It is really simple and the results are the same (see image below with Value Tool Plugin).

setwd('pyqgis_data') #where is my_raster
library(raster)
my_raster<-raster('utah_demUTM2.tif')
#Benedikt Gehr approach
centered_raster<-(my_raster-cellStats(my_raster,mean))/(cellStats(my_raster,sd))
#RobertH approach
centered_raster2 <- scale(my_raster)
writeRaster(centered_raster, 'centered_raster.tif', format='GTiff',overwrite=TRUE)
writeRaster(centered_raster2, 'centered_raster2.tif', format='GTiff',overwrite=TRUE)

However, these approaches require the rgdal library installed. Sometimes, Windows users have problems with its installation.