# Comparing two rasters, based on a complex set of rules

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?

• I guess some example data would be helpful. Also, since the resulting map will probably be very complex, from a cartographic point of view it'd be preferrable if you somehow reduced the categories per map (and thus reduce the complexity of the calculation. Addtional thought: you can chain expressions using `+`. Wouldn't that work out for you?
– Erik
Dec 22, 2021 at 11:25
• Thank you @Erik, The example data is [here] (github.com/datakeen/VariablesnRatios/blob/master/…). I am trying to work with your suggestion of using ‘+’ method. Dec 22, 2021 at 13:27
• Besides, can you suggest to get the range of a raster programmatically? That is in the Raster Calculator. Dec 22, 2021 at 13:36

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)

``````

## Result

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())
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
``````

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)
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)
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)
``````