1

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?

Matteo
  • 157
  • 4
  • 1
    (Thinking out loud) Could you load your csv into your project as a Delimited text layer with no geometry, and then access it natively without pandas? You would then build your "FORMULA" parameter in your python code, referencing the proper "feature" in the csv layer to pull the coefficient values? – Houska Feb 14 '20 at 15:12
  • Great idea, I've tried it but I can't still get the desired result. I've update the question with the new code that I've attempted – Matteo Feb 17 '20 at 17:34
  • Hmm. I'm not that experienced with applying the gdal calculator specifically. Ideas a) get it to debug-output formula_y to really check it's what you would use manually (you may have already done that), b) make sure your NO_DATA values are right (if applicable), c) give up on GDAL raster calculator and do a numpy calculation directly, see for instance https://gis.stackexchange.com/questions/204387/gdal-calc-raster-calculator-syntax-for-logical-operators-and-other-functions in another context. – Houska Feb 17 '20 at 18:34
  • You could use a model for that. – Babel Dec 28 '20 at 16:46
  • Can you expand on this? What do you mean with model? – Matteo Jan 05 '21 at 12:31

0 Answers0