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?

5
  • 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, 2020 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, 2020 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 gis.stackexchange.com/questions/204387/… in another context.
    – Houska
    Feb 17, 2020 at 18:34
  • You could use a model for that.
    – Babel
    Dec 28, 2020 at 16:46
  • Can you expand on this? What do you mean with model?
    – Matteo
    Jan 5, 2021 at 12:31

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.