Sign up ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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?

share|improve this question

2 Answers 2

up vote 12 down vote accepted

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!

share|improve this answer
It makes perfect sense, thank you! –  Rich Nov 13 '11 at 18:14
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 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
share|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

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