How to add different sized rasters in GDAL so the result is only in the intersected region

Question

I'm writing a Python method that adds two rasters and generates a single raster output. For reasons beyond my control, the extents of the input rasters are different, but they do overlap.

Is it possible to operate exclusively on the input raster pixels that are overlapped in the 2 overlapped regions to generate my ouput such that the output raster extent is only the intersecting region of the two inputs?

Answer 1

The first thing to do is to determine the overlapping rectangle in geospatial coordinates. To do this, you get the geotransform for each source image:

gt1 = ds1.GetGeoTransform()
# r1 has left, top, right, bottom of dataset's bounds in geospatial coordinates.
r1 = [gt1[0], gt1[3], gt1[0] + (gt1[1] * ds1.RasterXSize), gt1[3] + (gt1[5] * ds1.RasterYSize)]

# Do the same for dataset 2...

intersection = [max(r1[0], r2[0]), min(r1[1], r2[1]), min(r1[2], r1[2]), max(r1[3], r2[3])]

Then convert that rectangle into pixels for each image by subtracting the top and left coordinates and dividing by the pixel size, rounding up.

From here you can call ReadRaster() on each image, giving it the pixel extents you've just calculated:

band.ReadRaster(px1[0], px1[1], px1[2] - px1[0], px1[3] - px1[1], px1[2] - px1[0], px1[3] - px1[1], <band's datatype here>)

I'm a bit tired, so if this doesn't make much sense, let me know!

Answer 2

Noticed a typo. Third element of intersection should be min(r1[2], r2[2]):

intersection = [max(r1[0], r2[0]), min(r1[1], r2[1]), min(r1[2], r2[2]), max(r1[3], r2[3])]

Also, I would recommend some logic to check that the datasets actually intersect. This is one approach:

if (intersection[2] < intersection[0]) or (intersection[1] < intersection[3]):
    intersection = None