pythonrasterio

How to substract raster1 and raster2 in rasterio


I have tried to subtract the raster1.tif and raster2.tif using rasterio in python. The raster1.tif file is larger and completely overlaps the raster2.tif.

I need to subtract raster1.tif - raster2.tif in such a way so that the output raster extent will be equal to raster1.tif.

I have written the following code,

import rasterio 

with rasterio.open('raster1.tif') as src1, rasterio.open('raster2.tif') as src2:
    data1 = src1.read(1)
    data2 = src2.read(1)
    data = data1 - data2

    with rasterio.open("output.tif", 'w', **src1.meta) as dst:
        dst.write(data,1)

It produces the following error,

ValueError: operands could not be broadcast together with shapes (790,1554) (57,67) 

I know this error is generated since both my raster has different bounds. Can anyone help me to fix this issue? is there any way to reshape the data2 so I can subtract easily?

PS: I have reproduced the required thing with ArcGIS Pro using a raster calculator with the condition below,

Con(IsNull('raster2.tif'), 'raster1.tif', 'raster1.tif' - 'raster2.tif')

I am looking for the equivalent code in rasterio.

Image showing raster1 and raster2


Solution

  • I successfully implemented the solution code as below after an hours of trial and errors. Thank you everyone,

    with rasterio.open('raster1.tif') as src1, rasterio.open('raster2.tif') as src2:
        data1 = src1.read(1)
    
        # extract metadata from the first raster
        meta = src1.meta.copy()
    
        # read the window of the second raster with the same extent as the first raster
        window = src2.window(*src1.bounds)
    
        # read the data from the second raster with the same window as first raster
        data2 = src2.read(1, window=window, boundless=True, fill_value=0)
        data2 = np.where(data2 == src2.nodata, 0, data2)
    
        # calculate the difference
        data = data1 - data2
    
        # write the result to a new raster
        with rasterio.open("output.tif", 'w', **meta) as dst:
            dst.write(data,1)