3

In the sequence of my previous question (Can't use Raster Calculator with two different rasters), I am stuck again in my project.

With my source rasters, I created 3 different outputs using the following code:

calc100 = QgsRasterCalculator('aa@1*(1/28) + ao@1*(2/28) + am@1*(1/28) + ct@1*(6/28) + ln@1*(6/28) + ea@1*(2/28) + gl@1*(6/28) + gm@1*(2/28) + rs@1*(2/28)', 'D:/maps/z100.tif', 'GTiff', aco_ant.extent(), aco_ant.width(), aco_ant.height(), entries)
calc100.processCalculation()

calc100to5000 = QgsRasterCalculator('ln@1 * (6/28+0.665*(d@1-100)/4900) + ct@1 * (6/28-0.665*(d@1-100)*6/(4900*22)) + gl@1 * (6/28-0.665*(d@1-100)*6/(4900*22)) + gm@1 * (2/28-0.665*(d@1-100)*2/(4900*22)) + rs@1 * (2/28-0.665*(d@1-100)*2/(4900*22)) + ao@1 * (2/28-0.665*(d@1-100)*2/(4900*22)) + ea@1 * (2/28-0.665*(d@1-100)*2/(4900*22)) + am@1 * (1/28-0.665*(d@1-100)/(4900*22)) + aa@1 * (1/28-0.665*(d@1-100)/(4900*22))', 'D:/maps/z100to5000.tif', 'GTiff', aco_ant.extent(), aco_ant.width(), aco_ant.height(), entries)
calc100to5000.processCalculation()

calc5000 = QgsRasterCalculator('ln@1 * 0.879 + ct@1 * 0.033 + gl@1 * 0.033 + gm@1 * 0.011 + rs@1 * 0.011 + ao@1 * 0.011 + ea@1 * 0.011 + am@1 * 0.005 + aa@1 * 0.005', 'D:/maps/z5000.tif', 'GTiff', aco_ant.extent(), aco_ant.width(), aco_ant.height(), entries)
calc5000.processCalculation()

However, I just want to create one single output raster. I have one raster with the distance to shoreline at each point, and I was thinking if it's possible to insert an if statement to identify the value in that raster and use the respective formula to the pixel in question, something like this:

if d@1 < 100:
calc100 = QgsRasterCalculator('aa@1*(1/28) + ao@1*(2/28) + am@1*(1/28) + ct@1*(6/28) + ln@1*(6/28) + ea@1*(2/28) + gl@1*(6/28) + gm@1*(2/28) + rs@1*(2/28)', 'D:/maps/z100.tif', 'GTiff', aco_ant.extent(), aco_ant.width(), aco_ant.height(), entries)
calc100.processCalculation()

if d@1 >= 100 and d@1 <= 5000:    
calc100to5000 = QgsRasterCalculator('ln@1 * (6/28+0.665*(d@1-100)/4900) + ct@1 * (6/28-0.665*(d@1-100)*6/(4900*22)) + gl@1 * (6/28-0.665*(d@1-100)*6/(4900*22)) + gm@1 * (2/28-0.665*(d@1-100)*2/(4900*22)) + rs@1 * (2/28-0.665*(d@1-100)*2/(4900*22)) + ao@1 * (2/28-0.665*(d@1-100)*2/(4900*22)) + ea@1 * (2/28-0.665*(d@1-100)*2/(4900*22)) + am@1 * (1/28-0.665*(d@1-100)/(4900*22)) + aa@1 * (1/28-0.665*(d@1-100)/(4900*22))', 'D:/maps/z100to5000.tif', 'GTiff', aco_ant.extent(), aco_ant.width(), aco_ant.height(), entries)
calc100to5000.processCalculation()

if d@1 > 5000:    
calc5000 = QgsRasterCalculator('ln@1 * 0.879 + ct@1 * 0.033 + gl@1 * 0.033 + gm@1 * 0.011 + rs@1 * 0.011 + ao@1 * 0.011 + ea@1 * 0.011 + am@1 * 0.005 + aa@1 * 0.005', 'D:/maps/z5000.tif', 'GTiff', aco_ant.extent(), aco_ant.width(), aco_ant.height(), entries)
calc5000.processCalculation()

If it's not possible, could you recommend a viable way?

2

I found an alternative answer to my problem, so I think I should share it with anyone interested.

So, instead of trying to find a "if" function that worked, I created masks for each function, in order to only give results for the areas that I wanted to and 0 for the other zones.

calc100 = QgsRasterCalculator('(d@1 < 100)*(aa@1*(1/28) + ao@1*(2/28) + am@1*(1/28) + ct@1*(6/28) + ln@1*(6/28) + ea@1*(2/28) + gl@1*(6/28) + gm@1*(2/28) + rs@1*(2/28))', 'D:/maps/z100.tif', 'GTiff', aco_ant.extent(), aco_ant.width(), aco_ant.height(), entries)
calc100.processCalculation()

calc100to5000 = QgsRasterCalculator('(d@1 >= 100 AND d@1<=5000)*(ln@1 * (6/28+0.665*(d@1-100)/4900) + ct@1 * (6/28-0.665*(d@1-100)*6/(4900*22)) + gl@1 * (6/28-0.665*(d@1-100)*6/(4900*22)) + gm@1 * (2/28-0.665*(d@1-100)*2/(4900*22)) + rs@1 * (2/28-0.665*(d@1-100)*2/(4900*22)) + ao@1 * (2/28-0.665*(d@1-100)*2/(4900*22)) + ea@1 * (2/28-0.665*(d@1-100)*2/(4900*22)) + am@1 * (1/28-0.665*(d@1-100)/(4900*22)) + aa@1 * (1/28-0.665*(d@1-100)/(4900*22)))', 'D:/maps/z100to5000.tif', 'GTiff', aco_ant.extent(), aco_ant.width(), aco_ant.height(), entries)
calc100to5000.processCalculation()

calc5000 = QgsRasterCalculator('(d@1 > 5000)*(ln@1 * 0.879 + ct@1 * 0.033 + gl@1 * 0.033 + gm@1 * 0.011 + rs@1 * 0.011 + ao@1 * 0.011 + ea@1 * 0.011 + am@1 * 0.005 + aa@1 * 0.005)', 'D:/maps/z5000.tif', 'GTiff', aco_ant.extent(), aco_ant.width(), aco_ant.height(), entries)
calc5000.processCalculation()

Then, I only had to sum all the maps using the QgsRasterCalculator in the same way as above, since the maps with the mask don't overlap each others:

vuln = QgsRasterCalculator('v100@1 + v100to5000@1 + v5000@1', 'D:/maps/vuln.tif', 'GTiff', z100.extent(), z100.width(), z100.height(), entries)
vuln.processCalculation()

Thanks to all, Pedro

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.