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, gt, gt + (gt * smaller_raster.RasterXSize), gt + (gt * smaller_raster.RasterYSize)] r2 = [gt2, gt2, gt2 + (gt2 * larger_raster.RasterXSize), gt2 + (gt2 * larger_raster.RasterYSize)] #find where they intersect intersection = [max(r1, r2), min(r1, r2), min(r1, r2), max(r1, r2)] #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 - gt2) / gt2) yoffset = int(abs(intersection - gt2) / gt2) new_extent = lc.ReadAsArray(xoffset, yoffset, int(gt2), int(gt2))
but this returns:
ValueError: negative dimensions are not allowed
[-8895604.157333, 7783653.637667, -7783653.637667, 6671703.118]
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