Jak dodać rastry różnej wielkości w GDAL, aby wynik był tylko w przecinającym się regionie

11

Piszę metodę Python, która dodaje dwa rastry i generuje jeden wynik rastrowy. Z przyczyn niezależnych ode mnie zakresy rastrów wejściowych są różne, ale się pokrywają.

Czy możliwe jest działanie wyłącznie na wejściowych pikselach rastrowych, które nakładają się na 2 zachodzące na siebie regiony, w celu wygenerowania mojego wyjścia tak, że wyjściowy zasięg rastra jest tylko obszarem przecinającym się na dwóch wejściach?

python raster gdal Bogaty
źródło

Odpowiedzi:

12

Pierwszą rzeczą do zrobienia jest określenie nakładającego się prostokąta we współrzędnych geoprzestrzennych. Aby to zrobić, otrzymasz geotransformę dla każdego obrazu źródłowego:

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])]

Następnie przekonwertuj ten prostokąt na piksele dla każdego obrazu, odejmując górne i lewe współrzędne i dzieląc przez rozmiar pikseli, zaokrąglając w górę.

Stąd możesz wywoływać ReadRaster()każdy obraz, nadając mu zakres pikseli, który właśnie obliczyłeś:

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>
)

Jestem trochę zmęczony, więc jeśli to nie ma sensu, daj mi znać!

MerseyViking
źródło
Czy działa to również w przypadku rastrów o różnych rzutach (aka systemy odniesienia za pomocą współrzędnych / systemy odniesienia przestrzennego)? I nawet jeśli projekcje są takie same: Czy to działa również wtedy, gdy gt1[1]i gt2[1](lub gt1[5]i gt2[5]) mają przeciwne znaki? (Jak sądzę, który przewróciłby jeden z rastrów w pionie lub w poziomie). Albo jeśli abs(gt1[2])i abs(gt1[4])są większe niż abs(gt1[1])i abs(gt1[5])ale abs(gt2[2])i abs(gt2[4])są mniejsze niż abs(gt2[1])i abs(gt2[5])(co prawdopodobnie przewróciłoby jeden z rastrów po przekątnej)?
das-g
6

Trzecim elementem przecięcia powinno być min (r1 [2], r2 [2]):

intersection = [max(r1[0], r2[0]), min(r1[1], r2[1]), min(r1[2], r2[2]), max(r1[3], r2[3])]

Ponadto zaleciłbym pewną logikę, aby sprawdzić, czy zbiory danych faktycznie się przecinają.

To jest jedno podejście:

if (intersection[2] < intersection[0]) or (intersection[1] < intersection[3]):
    intersection = None
David Shean
źródło