How to calculate GRASS covariance matrices in QGIS with r.covar?

Question

I would like to calculate covariance matrices between different bands of the same raster in QGIS (2 or 3). I have found only one built-in processing tool capable of doing this, named r.covar from GRASS.

Unfortunately, in QGIS this tool is structured to only accept "entire" raster layers as input, while in GRASS it is indeed possible to calculate covariance matrices between bands of the same raster dataset.

I have tried to trick PyQGIS to accepting GDAL's bands as inputs, but it only accepts raster layers. Does anyone have a solution to this problem or perhaps a user script that does creates covariance and correlation matrices?

I would strongly prefer not to save each band as a separate file strictly for this problem.

Sample code:

import processing,gdal

rasterLayer = iface.activeLayer()
ds = gdal.Open("C:\\Multispectral_Clipped.tif")

outputRaster= "C:\\Output.html"
extent =rasterLayer.extent()

xmin = extent.xMinimum()
xmax = extent.xMaximum()
ymin = extent.yMinimum()
ymax = extent.yMaximum()

#Not working option
output=processing.runalg('grass7:r.covar', [ds.GetRasterBand(1),ds.GetRasterBand(2), ds.GetRasterBand(3)],True,"%f,%f,%f,%f" %(xmin, xmax, ymin, ymax),outputRaster,None)

#Working option
output=processing.runalg('grass7:r.covar', [rasterLayer],True,"%f,%f,%f,%f" %(xmin, xmax, ymin, ymax),outputRaster,None)

EDIT I suspect my request above is not possible. So the question is then just how to run r.covar with multiple inputs in PyQGIS (QGIS 2.18)?

Even with having the line:

output=processing.runalg('grass7:r.covar', [rasterList],True,"%f,%f,%f,%f" %(xmin, xmax, ymin, ymax),outputHtml,outputRaw)

And the rasterList containing a Python list like so: ['C:\\band1.tiff', 'C:\\band2.tiff', 'C:\\band3.tiff']

I am getting a Wrong parameter value error. I've also tried loading them as QgsRasterLayer, like done here with the StringToRaster function, but got the same error.

How did you define outputHtml and outputRaw? Does it work if you replace these with None? — Joseph, Mar 29 '18 at 11:28
These are simply defined as paths, e.g. outputHtml = "C:\\correlation_output.html" and outputRaw = "C:\\correlation_output.txt". I tried replacing them with None, it still doesn't work. Furthermore, the error says Wrong parameter value: ['C:\\band1.tiff', 'C:\\band2.tiff', 'C:\\band3.tiff'], pointing to the input error there. — 15Step, Mar 29 '18 at 11:36
They are .tiff files, saved from an original composite 'tif' file using this little loop: for i in range (1,ds.RasterCount+1): srcband = ds.GetRasterBand(i) out_ds = gdal.Translate(out_path + 'band' + str(i) + '.tiff', ds, format='GTiff', bandList=[i]) out_ds=None I tried saving them as .tif, but got the same error. Also, r.covar works in the interface with the rasters, both as .tif or tiff. — 15Step, Mar 29 '18 at 11:55

score 1 · Accepted Answer · answered Apr 26 '18 at 11:29

It appears the problem was in the order of the extent variables, as well as the rasterList variable type and format. Here is the corrected code:

def getRasterExtent(raster_path):
    src = gdal.Open(raster_path)
    ulx, xres, xskew, uly, yskew, yres  = src.GetGeoTransform()
    lrx = ulx + (src.RasterXSize * xres)
    lry = uly + (src.RasterYSize * yres)
    return [ulx,lrx,lry,uly] #xmin, xmax,ymin,ymax


rasterList=[]
extent=[]
for filename in os.listdir(input_dir): #Loop over all individual single band rasters in a folder and create correlation matrix
    if filename.endswith(".tif"): 
        rasterPath = os.path.join(input_dir, filename)
        rasterList.append(rasterPath)
        extent =getRasterExtent(rasterPath)


rasterString =str( ";".join(rasterList)) #convert the list of raster paths to a *semicolon* seperated STRING (the way GDAL expects it).
xmin,xmax,ymax,ymin = extent[0],extent[1],extent[2],extent[3]

output=processing.runalg('grass7:r.covar', rasterString,True,"%f,%f,%f,%f" %(xmin, xmax, ymax, ymin),None,None)

So this works and creates a correlation matrix from a given directory of single-band, same extent, rasters.

How to calculate GRASS covariance matrices in QGIS with r.covar?

1 Answers1