5

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?

3
  • 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.
    – datakeen
    Dec 22, 2021 at 13:27
  • Besides, can you suggest to get the range of a raster programmatically? That is in the Raster Calculator.
    – datakeen
    Dec 22, 2021 at 13:36

3 Answers 3

4

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

enter image description here

Source B

enter image description here

Difference

enter image description here

Result

enter image description here

3

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
0

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)

Your Answer

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

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