QGIS from console: raster algebra

Question

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)

Answer 1

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

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)

Answer 2

You find an excerpt of a script of mine. I am using GDAL there and the calculation is done with arrays. Important lines for you are the RGBI_DS - here you put your path and filename. and the line dataout - here you can replace the 2.0 by your constant. I hope this gets you started.

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

# Raster Dataset
RGBI_DS = "D:/Rohdaten/RGBI_DS.img" 

# 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 * 2.0)

# Output
outFile = "D:/NDVI.tif"

# 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