I need to calculate an expression over 25 different rasters. The function is the same for all rasters, but the coefficients are raster-specific. I would like to tell PyQGIS to select the relevant coefficients and then feed them into GDAL raster calculator formula.
To be more specific, I have rasters named as nl_y, where y is a numeric identifier, and coefficients relative to each y stored in a table (my_file_csv).
What I have tried to do - following Houska's advise - is to: (i) import the csv table as a layer without geometry; (ii) create a string reporting the expression to be applied in the raster calculator; (iii) feed this string into GDAL raster calculator.
Below is the code doing this:
## (i) Import csv with coefficients
uri = "file:///my_file_csv?delimiter=,"
L_coeff = QgsVectorLayer(uri, "csvLayer", "delimitedtext")
QgsProject.instance().addMapLayer(L_coeff)
for y in year:
## (ii) Select coefficients and compute formula
for f in L_coeff.getFeatures():
if f["Var name"] == "Var 0":
c0_y = f["Column " + y]
if f["Var name"] == "Var 1":
c1_y = f["Column " + y]
if f["Var name"] == "Var 2":
c2_y = f["Column " + y]
formula_y = "(" + str(c0_y) + ")* (A > 0) + (" + str(c1_y) + ")* A + (" + str(c2_y) + ")* (A)^2"
## (iii) Raster calculation
Fname_y = "nl_" + y + ".tif"
processing.run("gdal:rastercalculator", {
"INPUT_A": input_dir + "\\" + Fname_y,
"BAND_A": 1,
"FORMULA": formula_y,
"OUTPUT": output_dir + "\\" + Fname_y})
I get output rasters but they have no-sense values - i.e. between 1.79769e+308 and -1.79769e+308, instead of having 0-63 range. My guess is that GDAL raster calculator does not correctly read the formula, even though when I print it, it is exactly identical to what I would write manually in GDAL algorithm.
Any idea on how I can apply raster calculation (not necessarily on GDAL!) in an automated manner using Python console in QGIS?
"FORMULA"
parameter in your python code, referencing the proper "feature" in the csv layer to pull the coefficient values?formula_y
to really check it's what you would use manually (you may have already done that), b) make sure yourNO_DATA
values are right (if applicable), c) give up on GDAL raster calculator and do anumpy
calculation directly, see for instance gis.stackexchange.com/questions/204387/… in another context.