以编程方式减去不同维度的光栅

Question

问题：

如何以编程方式返回两个（大小不同）红色带差的光栅？在

即

gdal_calc.py -A 'WARPED.tif'  -B 'DSC_1636.tif' --outfile = 'dif.tif' --calc = "A-B"

QGIS光栅计算器很好地完成了这个功能。但是，前面的代码返回以下错误。在

Exception: Error! Dimensions of file DSC_1636.tif (7380, 4928) are different from other files (7743, 5507). Cannot proceed

我目前的印象是，我应该使用一个定义的范围，通过找到下面显示的重叠来创建光栅，但我仍然不能使这个工作。在

^{pr2}$

warped image size is (7743, 5507)
initial image (image 2) size is (7380, 4928)

# Use GDAL to get information about images
def get_extent(fn):
    '''Returns min_x, max_y, max_x, min_y'''
    ds = gdal.Open(fn)
    gt = ds.GetGeoTransform()
    return (gt[0], gt[3], gt[0] + gt[1] * ds.RasterXSize,
        gt[3] + gt[5] * ds.RasterYSize)

print('extent of warped.tif is %s' % str(get_extent('DSC_0934-warped.tif')))
print('extent of 1636.png is %s' % str(get_extent('DSC_1636.png')))

extent of warped.tif is (-375.3831214210602, 692.5167764068751, 7991.3588371542955, -5258.102875649754)
extent of 1636.png is (0.0, 0.0, 7380.0, 4928.0)

r1 = get_extent('DSC_0934-warped.tif')
r2 = get_extent('DSC_1636.png')

# Get left, top, right, bottom of dataset's bounds in pixel coordinates
intersection = [max(r1[0], r2[0]), 
                min(r1[1], r2[1]), 
                min(r1[2], r2[2]), 
                max(r1[3], r2[3])]

print('checking for overlap')
if (intersection[2] < intersection[0]) or (intersection[1] > intersection[3]):
    intersection = None
    print('no overlap')
else:
    print('intersection overlaps at: %s' % intersection)

^{4磅}$

Answer 1

最直截了当的回答是将图像作为一个定义维度的数组读入。在

无需重新发布上面用于检查重叠位置的代码，就可以通过以下添加来获得解决方案。（谢谢@Val）

# Get the data
ds1_src = gdal.Open( "DSC_1636.png" )
ds2_src = gdal.Open( "DSC_0934-warped.tif")

ds1_bnd = ds1_src.GetRasterBand(1).ReadAsArray(xoff=0, yoff=0, win_xsize=7380, win_ysize=4928)
ds2_bnd = ds2_src.GetRasterBand(1).ReadAsArray(xoff=0, yoff=0, win_xsize=7380, win_ysize=4928)

# Do the maths...
data_out = ds2_bnd - ds1_bnd

#Write the out file
driver = gdal.GetDriverByName("GTiff")
dsOut = driver.Create("out.tiff", 7380, 4928, 1, GDT_Byte)
CopyDatasetInfo(ds1_src,dsOut)
bandOut=dsOut.GetRasterBand(1)
BandWriteArray(bandOut, data_out)

#Close the datasets
ds1_src = None
ds2_src = None
ds1_bnd = None
ds2_bnd = None
bandOut = None
dsOut = None