Comparing two rasters, based on a complex set of rules

Question

I am trying to compare two rasters, say A.tif and B.tif. Both correspond to river bed level changes under two different flow conditions. Now I want to produce a comparison raster having values based on certain rules. The rules for the value of pixel are following

Pixel Value	Rule
0	Both A and B does not cause any significant bed level change (less than 10% of Min (Range_A, Range_B)
1	Both A and B are positive and are equal (the difference is less than 10% of Min(Range_A, Range_B)
2	Both A and B are negative and are equal (difference less than 10% of Min (Range_A, Range_B)
3	A is positive and B is negative (A > 0 and B < 0)
4	A is negative and B is positive (A < 0 and B > 0)
5	Both A and B are positive. Positive change due to A is more than that of B (i.e., A > B)
6	Both A and B are positive. Positive change due to B is more than that of A (i.e., A < B)
7	Both A and B are negative. Negative change due to A is more than that of B (i.e., A < B)
8	Both A and B are negative. Negative change due to B is more than that of A (i.e., A > B)

The raster calculator in the current version of QGIS (3.22.2) supports conditional if statements. Although, it is possible to do it in the raster calculator using a complex and lengthy expression. Probably by looping many if statements in each other. The problem with using Raster Calculator is that it cannot programmatically get the range of a certain raster. I have to take maximum and minimum values by myself (from the properties of the raster file) and calculate the range myself for each comparison I have to do. I have tried other open source Raster Calculators. I have tried Raster calculators for SAGA, GRASS and gvSIG. But this limitation is also there. The Raster calculator's expression (without following proper syntax) approximately takes the following form

A  = Raster layer A
B = Raster Layer B
Max value in A = MAX (A)
Min value in A = MIN (A)
R_A = Range of Cell Values in Layer A = (MAX (A) - MIN (A))
Max value in B = MAX (B)
Min value in B = MIN (B)
R_B = Range of Cell Values in Layer B = (MAX (B) - MIN (B))

If0 ( 
    ((A) < 0.1 * Min(R_A, R_B) AND (B) < 0.1 * Min(R_A, R_B))
    AND 
    ((A) > -0.1 * Min(R_A, R_B)) AND ((B) > -0.1 * Min(R_A, R_B)), 0, 
    If1 ( 
        ((A > 0) and (B > 0)) 
        AND
        (A < 0.1 * Min_Range + B) AND (A > 0.1 * Min_Range - B)), 1,
        If2 (((A < 0) AND (B < 0)) 
             AND
            (A < 0.1 * Min_Range + B) AND (A > 0.1 * Min_Range - B), 2
        ...

With this, the script becomes verbose and convoluted, and secondly, there are many more similar comparisons therefore I am searching for some scripting (either in Python or in R). So that, the method may be easily repeatable and reproducible. I have a beginner level familiarity with scripting (Python and R). I am wondering if there is an easy way to do it in the scripting?

Answer 1

Here is a little documented R sketch using the rgdal and raster package.

# Load the raster driver kit
require(rgdal)
require(raster)
Set the work dir
WORK_DIR <- getwd()
DATA_DIR <- WORK_DIR
Create the file names
file.a <- sprintf('%s/%s', DATA_DIR,'A.tif')
file.b <- sprintf('%s/%s', DATA_DIR,'B.tif')
file.c <- sprintf('%s/%s', DATA_DIR,'C.tif')
Copy the file A to get the metadata for the result
like resolution, CRS etc. (quick & dirty variant)
file.copy(file.a, file.c)
Show what in the input
GDALinfo(file.a)
GDALinfo(file.b)
Open the source files
A = GDAL.open(file.a)
B = GDAL.open(file.b)
Open the result file writeable
C = GDAL.open(file.c, read.only = FALSE)
Plot the source files
displayDataset(A)
displayDataset(B)
Extract the raster data
mx.A <- getRasterData(A)
mx.B <- getRasterData(B)
mx.C <- getRasterData(C)
Build the difference
mx.D <- mx.A-mx.B
Show the differences
plot(raster(mx.D))
hist(mx.D)
Get the range of the differences to build the threshold
range.dif <- range(mx.D, na.rm = TRUE)
Not sure if the right threshold is applied here
thresh    <- diff(range(mx.D, na.rm=TRUE))/10
Clear all in data in MX C
mx.C[,] <- 0
Apply the rules ..not all implemented ..example only
You could use percentile based rules here.
mx.C[  (mx.D > -thresh) & (mx.D < thresh)   ]         <- 10
mx.C[ (mx.A > 0) & (mx.B > 0) & (  mx.D < thresh )  ] <- 20
mx.C[ (mx.A < 0) & (mx.B < 0) & ( -mx.D < thresh )  ] <- 30
mx.C[ (mx.A > 0) & (mx.B < 0) ]                       <- 40
mx.C[ (mx.A < 0) & (mx.B > 0) ]                       <- 50
mx.C[ (mx.A > 0) & (mx.B > 0) & (mx.A > mx.B) ]       <- 60
mx.C[ (mx.A > 0) & (mx.B > 0) & (mx.A < mx.B) ]       <- 70
etc. ... many more
Store the raster back to the geotiff
putRasterData(C, mx.C)
Show the geotiff
displayDataset(C)
Store the geotiff
saveDataset(C, file.c)

Source A

Source B

Difference

Result

Answer 2

You can use numpy, example:

import numpy as np
from osgeo import gdal
a = QgsProject.instance().mapLayersByName('A')[0]
b = QgsProject.instance().mapLayersByName('B')[0]
def givearray(rasterlayer):
    ds = gdal.Open(rasterlayer.source())
    myarray = np.array(ds.GetRasterBand(1).ReadAsArray())
    return myarray
aarr = givearray(a)
barr = givearray(b)
rule3 = np.where((aarr>=0) & (barr<0), 3, None)
#Other rules
#Then combine the arrays and write to raster: https://gis.stackexchange.com/questions/37238/writing-numpy-array-to-raster-file

Answer 3

Both the answers, by @huckfinn and @BERA are worthy and helped me to reach the final code. For the sake of completeness, I am writing the final code for helping anyone else reaching here.

import numpy as np
from osgeo import gdal
import rasterio
def givearray(rasterlayer): #Checked
    ds = gdal.Open(rasterlayer)
    myarray = np.array(ds.GetRasterBand(1).ReadAsArray())
    return myarray
def range(rasterlayer):
    # Calculates the range of a raster layer
    ds = gdal.Open(rasterlayer)
    ds_min = ds.GetRasterBand(1).GetStatistics(0,1)[0]
    ds_max = ds.GetRasterBand(1).GetStatistics(0,1)[1]
    range = ds_max - ds_min
    return range
def threshold(rasterlayer1, rasterlayer2):
    # calculates the threshold value from two given raster layers
    threshold = 0.1 * min(range(rasterlayer1), range(rasterlayer2))
    return threshold
def applyRules (rasterlayer1, rasterlayer2):
    A = givearray(rasterlayer1)
    B = givearray(rasterlayer2)
    thld = threshold(rasterlayer1, rasterlayer2)
    rule1 = np.where((((A < thld) & (A >-thld)) & ((B < thld) & (B >-thld))), 1, 0)
    rule2 = np.where(((A > 0) & (B > 0)) & ((A < (B + thld)) & (A > (B - thld))), 2, 0)
    rule3 = np.where(((A < 0) & (B < 0)) & ((A < (B + thld)) & (A > (B - thld))), 3, 0)
    rule4 = np.where(((A > 0) & (B < 0)), 4, 0)
    rule5 = np.where(((A < 0) & (B > 0)), 5, 0)
    rule6 = np.where(((A > 0) & (B > 0)) & ((A >=  (B + thld))), 6, 0)
    rule7 = np.where(((A > 0) & (B > 0)) & ((A <=  (B - thld))), 7, 0)
    rule8 = np.where(((A < 0) & (B < 0)) & ((A <=  (B - thld))), 8, 0)
    rule9 = np.where(((A < 0) & (B < 0)) & ((A >=  (B + thld))), 9, 0)
Final = rule1 + rule2 + rule3 + rule4 + rule5 + rule6 + rule7 + rule8 + rule9 
return Final


def saveResult(rasterlayer1, Final, Output_File):
    raster_file = rasterio.open(rasterlayer1)
    raster_array = raster_file.read()
    driver = "GTiff"
    dim = raster_array.shape #For rows and columns of the raster array
    height = dim[1]
    width = dim[2]
    count = 1
    crs = raster_file.crs
    gdal_transform = raster_file.get_transform() #Returns a GDAL geotransform in its native form.
    # For (west, north, x-resolution, y-resolution)
    transform = rasterio.transform.from_origin(gdal_transform[0],gdal_transform[3],gdal_transform[1],-gdal_transform[5])
    data = Final.reshape(1, Final.shape[0], Final.shape[1]) #Reshaping the matrix so that its shape matches rasterio's requirements
    with rasterio.open(Output_File, "w", 
                        driver=driver, width=width, height=height, 
                        count=1, dtype="int32", crs=crs, transform=transform) as dst:
        dst.write(data)
def saveResult_alternative(rasterlayer1, Final, Output_File):
    # register all of the GDAL drivers
    gdal.AllRegister()
    reference_raster_data = gdal.Open(rasterlayer1)
    rows = reference_raster_data.RasterYSize
    cols = reference_raster_data.RasterXSize
    driver = reference_raster_data.GetDriver()
    raw_File = driver.Create(Output_File, cols, rows, 1, gdal.GDT_Int32)
    raw_File_band = raw_File.GetRasterBand(1)
    # Write the data
    raw_File_band.WriteArray(Final,0,0)
#Write the data to the disk
raw_File_band.FlushCache()

# Georeferencing and Projection
raw_File.SetGeoTransform(reference_raster_data.GetGeoTransform())
raw_File.SetProjection(reference_raster_data.GetProjection())


def totalComparison (rasterlayer1, rasterlayer2, Output_File):
    Final = applyRules(rasterlayer1,rasterlayer2)
    saveResult(rasterlayer1, Final, Output_File)
a = "A.tif" #Assuming the the raster files are in the current working directory
b = "B.tif"
Output_File = "Result_Final.tif"
totalComparison(a,b,Output_File)