Python – Comparing Two Rasters Based on Complex Rules with Python

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