Calculating difference between two raster layers in Python

Question

I have two raster layers(same extent) that represent land cover classes. So, I want to extract only class 1 and calculate the difference between class 1.

The expected output should be like the figure above.

I provide rasters example here

So far I just load raster using rasterio.open

import rasterio

class2019 = rasterio.open('classification_base.tif')
class2021 = rasterio.open('classification_target.tif')

Can you provide any ideas or example code?

Answer 1

You can read them as numpy arrays, subtract, and write the result to file.

I dont know what you want to do where the result is -1, or where the classes are 2 or 3. You can modify the code to do what you want.

import rasterio
import numpy as np

class2019 = rasterio.open("/home/bera/Downloads/classification_base.tif-20230721T095311Z-001/classification_base.tif")
a2019 = class2019.read(1) #Read first band as a numpy array
class2021 = rasterio.open("/home/bera/Downloads/classification_target.tif-20230721T095949Z-001/classification_target.tif")
a2021 = class2021.read(1)

a2019[np.isin(a2019, [2,3])]=np.nan #Where the value is 2 or 3, set it to nan so the result will be nan when subtracting.
a2021[np.isin(a2021, [2,3])]=np.nan

adiff = a2021-a2019 #Subtract 2021 from 2019
adiff[np.isnan(adiff)]=99 #Where either a2021 or a2019 is nan, the result is nan. Set it to 99.

#Save the result to file
output = "/home/bera/Downloads/target_diff_base.tif"

with rasterio.open(
    output,
    'w',
    driver='GTiff',
    height=adiff.shape[0],
    width=adiff.shape[1],
    count=1,
    dtype=adiff.dtype,
    crs=class2019.crs,
    transform=class2019.transform,
    compress="lzw") as dst:
    dst.write(adiff, 1)

class2019.close()
class2021.close()