2

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

Also,

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:

4800

and yoffset returns:

2400

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. – Michael Stimson Jul 6 '17 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]) ? – Stefano Potter Jul 6 '17 at 21:23
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, privacy policy and cookie policy

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