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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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