Программно вычитать растры разных размеров - PullRequest
Новая
       35

Программно вычитать растры разных размеров

0 голосов
/ 16 октября 2018

Вопрос:

Как программно вернуть растр, который отличается от двух (разных по размеру) красных полос?

т.е. 

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

Растровый калькулятор QGIS выполняет эту функцию просто отлично.Однако предыдущий код возвращает следующую ошибку.

Исключение: ошибка!Размеры файла DSC_1636.tif (7380, 4928) отличаются от других файлов (7743, 5507).Невозможно продолжить

В настоящее время у меня сложилось впечатление, что я должен читать растры, используя определенный экстент, созданный путем нахождения перекрытия, как показано ниже, но я все еще не могу выполнить эту работу. 

# Subtract two rasters of different dimensions
# Pixel coordinates define overlap

import os, sys
from PIL import Image
from osgeo import gdal, ogr, osr
gdal.UseExceptions()

    # Use PIL to get information from images
im1 = Image.open('DSC_0934-warped.tif')
print('warped image size is %s ' % str(im1.size))

im2 = Image.open('DSC_1636.png')
print('initial image (image 2) size is %s' % str(im2.size))

размер деформированного изображения составляет (7743, 5507)
исходное изображение (изображение 2), размер составляет (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')))

Степень warped.tif составляет (-375,3831214210602, 692,5167764068751, 7991,3588371542955, -5258.102875649754)
степень 1636.png составляет (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)
проверка перекрытия
пересечения пересечений в: [0,0, 0,0, 7380,0, 4928,0]

1 Ответ

0 голосов
/ 17 октября 2018

Самый простой ответ - читать изображения как массив заданных размеров.

Без повторного размещения приведенного выше кода, используемого для проверки того, где находится перекрытие, решение может быть найдено со следующими дополнениями.(Спасибо, @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
...