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

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?

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], r2[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!

• Does this also work for rasters with different projections (a.k.a. coordinate reference systems / spatial reference systems)? And even if the projections are the same: Does this also work if `gt1[1]` and `gt2[1]` (or `gt1[5]` and `gt2[5]`) have opposing signs? (Which would flip one of the rasters vertically or horizontally, I think.) Or if `abs(gt1[2])` and `abs(gt1[4])` are greater than `abs(gt1[1])` and `abs(gt1[5])` but `abs(gt2[2])` and `abs(gt2[4])` are smaller than `abs(gt2[1])` and `abs(gt2[5])` (which probably would flip one of the rasters diagonally)? – das-g Oct 6 '15 at 7:20

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