I two rasters with differen't dimensions, but they overlap. I want to change the extent of the larger raster so that it is 'clipped' to the smaller raster so I can stack them.

I found out how to determine where the intersection from this post, How to add different sized rasters in GDAL so the result is only in the intersected region

but I think I am getting something wrong in the final step or two.

So far I have this:

#read in the files
larger_raster = r'C:\large.tif'
larger_raster = gdal.Open(larger_raster)

smaller_raster = r'C:\small.tif'
smaller_raster = gdal.Open(smaller_raster)

#get the geotransformation information
gt = smaller_raster.GetGeoTransform()
larger_raster = lc.GetGeoTransform()

#get the xmin, xmax, ymin, ymax of each
r1 = [gt[0], gt[3], gt[0] + (gt[1] * smaller_raster.RasterXSize), gt[3] + (gt[5] * smaller_raster.RasterYSize)]

r2 = [gt2[0], gt2[3], gt2[0] + (gt2[1] * larger_raster.RasterXSize), gt2[3] + (gt2[5] * larger_raster.RasterYSize)]

#find where they intersect
intersection = [max(r1[0], r2[0]), min(r1[1], r2[1]), min(r1[2], r2[2]), max(r1[3], r2[3])]

#only select the locations where they intersect from the larger raster
new_extent = larger_raster.ReadAsArray(intersection)

this last line is where I go wrong though and this error is returned:

TypeError: in method 'BandRasterIONumPy', argument 3 of type 'double'

I have also tried this (picking up after finding the intersection):

xoffset = int(abs(intersection[0] - gt2[0]) / gt2[1]) 
yoffset = int(abs(intersection[1] - gt2[3]) / gt2[5])

 new_extent = lc.ReadAsArray(xoffset, yoffset, int(gt2[1]), int(gt2[5]))

but this returns:

ValueError: negative dimensions are not allowed


gt2 returns:

(-11119505.203637687, 463.3127136230469, 0.0, 8895604.157333, 0.0, -463.3127136230469)

intersection returns:

[-8895604.157333, 7783653.637667, -7783653.637667, 6671703.118]

xoffset returns:


and yoffset returns:


The code runs when I change the new extent calculation to this:

new_extent = large_raster.ReadAsArray(xoffset, yoffset, 1, 1)

but the output shape is (1L, 1L) whereas the shape of smaller_raster is (2400L, 2400L).

  • You are using georeference coordinates where the ReadAsArray are in rows and columns of the raster being read. Find the intersection area of r1 and r2 then calculate the offset, rows and columns to read. Jul 6, 2017 at 21:05
  • something like this (I have negative x values so I use absoute value) xoffset = abs(intersection[0] - gt2[0]) and yoffset = abs(intersection[1] - gt2[3]) and then new_extent = larger_raster.ReadAsArray(xoffset, yoffset, gt[1], gt[5]) ? Jul 6, 2017 at 21:23

1 Answer 1

from osgeo import gdal

# given that a geotransform object is in this form:
#adfGeoTransform[0] /* top left x */
#adfGeoTransform[1] /* w-e pixel resolution */
#adfGeoTransform[2] /* 0 */
#adfGeoTransform[3] /* top left y */
#adfGeoTransform[4] /* 0 */
#adfGeoTransform[5] /* n-s pixel resolution (negative value) */

smaller_raster = gdal.Open(r'C:\small.tif')
larger_raster  = gdal.Open(r'C:\large.tif')

gt = smaller_raster.GetGeoTransform()
lt = larger_raster.GetGeoTransform()

SmlMaxX = gt[0] + (gt[1] * smaller_raster.RasterXSize)
SmlMinY = gt[3] + (gt[5] * smaller_raster.RasterYSize)
Xoff = int((gt[0] - lt[0])/lt[1]) # cols to skip
Yoff = int((gt[3] - lt[3])/lt[5]) # rows to skip
Cols = int((SmlMaxX - gt[0])/lt[1])
Rows = int((SmlMinY - gt[3])/lt[5])

# print out some numbers so you can check manually
print("X offset {}, Y offset {}".format(Xoff,Yoff))
print("Xmax {}, Ymin {}".format(SmlMaxX,SmlMinY))
print("Reading {} cols, {} rows".format(Cols,Rows))

band = larger_raster.GetRasterBand(1)
data = larger_raster.ReadAsArray(Xoff,Yoff,Cols,Rows) # read the larger raster

It's not good to reuse your variables. In the start of your script larger_raster is a string then it's a gdal dataset then it's a geotransform array..

When you read an extent from a raster you need to do some checks to make sure that you're not addressing lower than 0 in either rows or columns and not exceeding the number of rows or columns in the raster you're reading. I have not added such checks as I believe it would confuse the intent of your question.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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