[GIS] QGIS from console: raster algebra

map-algebrapyqgis

I'm newbie to programming in pyqgis.

I need to perform a simple "GeoTiff * A" type of calculation(where A is a constant), but I am struggling to do so.

I tried processing.runalg('saga:rastercalculator') but failed 'coz my constant is NOT a raster.

I suppose I should be using 'qgismodelertools:calculator' or possibly something like to what they use here http://geoexamples.blogspot.it/2012/12/raster-calculations-with-gdal-and-numpy.html

Any clues?

Gene's answer iss succinct and logical but rasterio does not come out of the QGIS 2.8 box, which I'm using from OSGEO4W.

In the meantime, I used EikeMike's useful code by turning it into a working, callable function:

def raster_calc(RGBI_DS,outFile,constant):

  import processing, os, time, subprocess, array
  from osgeo import gdal  
  from osgeo.gdalnumeric import *  
  from osgeo.gdalconst import *

  # Read Raster dataset
  GDALRgbiDS = gdal.Open(RGBI_DS, GA_ReadOnly )
  band1rot = GDALRgbiDS.GetRasterBand(1)  

  # Convert band to Array
  Rdata = BandReadAsArray(band1rot) * 1.0

  # calculate raster * constant
  dataOut = (Rdata * constant)

  # Export result as Geotiff
  driver = gdal.GetDriverByName("GTiff")  
  dsOut = driver.Create(outFile, GDALRgbiDS.RasterXSize, GDALRgbiDS.RasterYSize, 1, 
   gdal.GDT_Float32)  
  CopyDatasetInfo(GDALRgbiDS,dsOut)  
  bandOut=dsOut.GetRasterBand(1)  
  BandWriteArray(bandOut, dataOut)

  #Deletion of variables, this actually triggers the writing of the output file
  dsOut =None
  bandOut = None
  driver =None

Later, I found myself a much easier and shorter solution.

Here it is:

Stick to 'saga:rastercalculator' (or similar algorithm). Simply build the mathematical formula therein as a new variable where your constant is itself a variable. In other words:

 import processing

 formula = ('a * %s' % (str(constant))) # where 'a' is going to point to yourraster   and 'constant' is the numeral you want to multiply your raster by
 processing.runalg('saga:rastercalculator',yourraster,None,formula,outFile)

Best Answer

Unfortunately the question is about PyQGIS and not GDAL (osgeo) (use rasterio, easier, see below in 3) and there are many examples in GIS Stack Exchange as How to evaluate raster calculator expressions from the console?

processing.runalg('saga:rastercalculator') is the Saga raster Calculator processing.runalg('grass:r.mapcalculator') is the GRASS GIS raster calculator

To use the QGIS raster calculator in the console, I detail here the process using the example of Cómo usar el Raster Calculator desde la Python Console de José Guerrero (in Spanish)

1) The raster calculator dialog

enter image description here

2) The script in the console

Characteristics of the raster layer

layer=iface.activeLayer()
# compare with the values in the dialog
print layer.extent().xMinimum(), layer.extent().xMaximum()
202086.577 205625.414407
print layer.extent().yMinimum(), layer.extent().yMaximum()
88411.048 90534.3504441
print layer.width(), layer.height()
100 60
name = layer.name()
print name
test
bands = layer.bandCount()
print bands
1

To use the raster calculator with the only band, the expression is therefore test@1 as in the dialog 

from qgis.analysis import QgsRasterCalculator, QgsRasterCalculatorEntry
resultinglayer = QgsRasterCalculatorEntry()
resultinglayer.ref = name + "@1"
resultinglayer.raster = layer
resultinglayer.bandNumber = 1
entries = [resultinglayer]
# perform a simple "GeoTiff * A" type of calculation(where A is a constant)
calc = QgsRasterCalculator('(test@1)*3', 
                          '/Users/Shared/resulting.tif', 
                          'GTiff', 
                          layer.extent(), 
                          layer.width(), 
                          layer.height(), 
                          entries )
calc.processCalculation()
# load the resulting raster layer.
rlayer = QgsRasterLayer('/Users/Shared/resulting.tif', 'resulting')
QgsMapLayerRegistry.instance().addMapLayer(rlayer)

3) EikeMike gives you a solution with GDAL, compare with a rasterio solution

import rasterio
with rasterio.open('/Users/Shared/test.tif') as original:
    layer = original.read() # or read_band(1), but deprecated
    resultinglayer = layer*3
    kwargs = original.meta
    with rasterio.open("/Users/Shared/resulting.tif", 'w', **kwargs) as output:
         output.write(resultinglayer)
Related Question