Take the 2-minute tour ×
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 11 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], 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!

share|improve this answer
    
It makes perfect sense, thank you! –  Rich Nov 13 '11 at 18:14

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
share|improve this answer

Your Answer

 
discard

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.