Первое, что нужно сделать, это определить перекрывающийся прямоугольник в геопространственных координатах. Для этого вы получите геотрансформацию для каждого исходного изображения:
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], r2[2]), max(r1[3], r2[3])]
Затем преобразуйте этот прямоугольник в пиксели для каждого изображения, вычтя верхнюю и левую координаты и разделив на размер пикселя, округлив в большую сторону.
Отсюда вы можете вызвать ReadRaster()
каждое изображение, указав только что рассчитанные пиксельные экстенты:
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>
)
Я немного устал, поэтому, если это не имеет особого смысла, дайте мне знать!
gt1[1]
иgt2[1]
(илиgt1[5]
иgt2[5]
) имеют противоположные признаки? (Какой бы переворачивать один из растров по вертикали или по горизонтали, я думаю.) Или , еслиabs(gt1[2])
иabs(gt1[4])
больше , чемabs(gt1[1])
и ,abs(gt1[5])
ноabs(gt2[2])
иabs(gt2[4])
меньше ,abs(gt2[1])
иabs(gt2[5])
(что , вероятно, переверните один из растров по диагонали)?Третий элемент пересечения должен быть min (r1 [2], r2 [2]):
Кроме того, я бы порекомендовал некоторую логику, чтобы убедиться, что наборы данных фактически пересекаются.
Это один из подходов:
источник