Algorytmy danych geoprzestrzennych
Przetwarzanie danych rastrowych

Krzysztof Dyba

Rozpocznijmy od wczytania danych rastrowych z pliku DEM.tif i wyświetlenia jego podstawowych metadanych.

import os
from qgis.core import QgsRasterLayer

# wczytanie rastra
filepath = os.path.join("dane", "DEM.tif")
raster = QgsRasterLayer(filepath, "DEM")
print(raster.isValid())
True
# wyświetlenie metadanych
print("Liczba kolumn (Szerokość):", raster.width())
print("Liczba wierszy (Wysokość):", raster.height())
print("Liczba komórek:", raster.height() * raster.height())
print("Liczba kanałów:", raster.bandCount())
print("Zakres:", raster.extent().toString(precision = 2))
print("CRS:", raster.crs().authid())
Liczba kolumn (Szerokość): 533
Liczba wierszy (Wysokość): 608
Liczba komórek: 369664
Liczba kanałów: 1
Zakres: 253698.33,353734.36 : 520058.23,657570.76
CRS: EPSG:2180

Odczyt wartości

Wartości rastra można wczytać na dwa sposoby – definiując punkt (QgsPointXY) lub blok (QgsRectangle). W drugim przypadku, jeśli wskażemy wszystkie kolumny oraz wiersze, wtedy do pamięci zostanie wczytany cały raster.

Punkt

W pierwszym kroku należy utworzyć obiekt QgsPointXY reprezentujący punkt, dla którego zostanie pozyskana wartość komórki. Konieczne jest określenie długości (X) oraz szerokości (Y) geograficznej punktu w układzie odniesienia warstwy rastrowej, tj. układ punktu oraz rastra muszą być identyczne! Następnie należy wykorzystać metodę sample() z klasy QgsRasterDataProvider i wskazać punkt oraz indeks kanału zaczynający się od 1. W przypadku gdy punkt znajduje się poza zasięgiem warstwy lub podano nieprawidłowy indeks kanału, to zostanie zwrócona wartość NaN. Dla większej liczby kanałów i/lub punktów należy zastosować pętlę.

from qgis.core import QgsPointXY

point = QgsPointXY(389900, 507600) # długość, szerokość
value = raster.dataProvider().sample(point, 1)
print(value)

# punkt jest poza zakresem rastra
point = QgsPointXY(0, 0)
value = raster.dataProvider().sample(point, 1)
print(value)
(141.6381072998047, True)
(nan, False)

Blok

Zamiast pojedynczego punktu możemy wykorzystać blok, który umożliwi pozyskanie wartości z zakresu prostokąta zdefiniowanego przez QgsRectangle. Do jego utworzenia niezbędne jest określenie minimalnych i maksymalnych współrzędnych (cztery wartości). Następnie użyjemy metodę block() również z klasy QgsRasterDataProvider, która wymaga wskazania:

  • indeksu kanału,
  • obiektu QgsRectangle,
  • liczby kolumn oraz wierszy.
from qgis.core import QgsRectangle

provider = raster.dataProvider()
rect = QgsRectangle(346950, 454028, 393355, 493871) # xmin, ymin, xmax, ymax

# określenie liczby kolumn i wierszy
width = int(rect.width() / raster.rasterUnitsPerPixelX())
height = int(rect.height() / raster.rasterUnitsPerPixelY())

# pobranie wartości pikseli dla prostokąta
block = provider.block(1, rect, width, height)
print(block.isValid())
print("Liczba komórek:", block.width() * block.height())
True
Liczba komórek: 7268

W rezultacie otrzymamy QgsRasterBlock, dla którego możemy wyświetlić numera kolumny, wiersza i odpowiadającą im wartość komórki.

for x in range(0, 3): # block.width()
    for y in range(0, 3): # block.height()
        print(x, y, round(block.value(x, y), 2))
0 0 79.54
0 1 83.73
0 2 85.62
1 0 78.35
1 1 78.12
1 2 80.94
2 0 77.76
2 1 78.12
2 2 80.94

Docinanie rastra

Wczytany blok w rzeczywistości jest niczym innym jak rastrem dociętym przez zakres przestrzenny. Wynik tej operacji możemy zapisać na dysku.

from qgis.core import QgsRasterFileWriter, QgsRasterPipe, QgsCoordinateTransformContext

pipe = QgsRasterPipe()
provider = raster.dataProvider()
pipe.set(provider.clone())

writer = QgsRasterFileWriter("DEM_crop.tif")
writer.setCreateOptions(["COMPRESS=LZW"])
writer.setOutputFormat("GTiff")

status = writer.writeRaster(
    pipe,
    width,
    height,
    rect,
    raster.crs(),
    QgsCoordinateTransformContext()
)

if status != 0:
    print("Błąd zapisu")
else: 
    print("OK")

Raster

Jeśli posiadamy wystarczającą ilość pamięci RAM, to możemy wczytać cały raster do pamięci. Wygląda to podobnie jak w przypadku wczytywania za pomocą bloku, z tą różnicą iż podajemy cały zakres rastra oraz wszystkie kolumny i wiersze. Możliwe jest dokonanie konwersji z QgsRasterBlock do macierzy numpy używając metody frombuffer() oraz reshape(). Zwróć uwagę, że w numpy pierwszy wymiar macierzy to wiersze, a drugi to kolumny (w QGIS jest na odwrót). Domyślnie dane nie są maskowane (wymagane jest zastosowanie modułu masked arrays).

import numpy as np

extent = raster.extent()
cols = raster.width()
rows = raster.height()
block = provider.block(1, extent, cols, rows)
array = np.frombuffer(block.data(), dtype = np.float32)
array = array.reshape(rows, cols)
print(array.shape) # wiersze, kolumny
print(array[100:104, 100:104])
(608, 533)
[[ 99.89686  117.18553  122.53293  132.92046 ]
 [104.782074 121.91743  122.99693  122.71936 ]
 [122.32443  120.56837   98.11017  105.36948 ]
 [100.51226  113.48083  120.141685 105.36948 ]]

W QGIS 3.40 dodano nową metodę as_numpy(), która upraszcza konwersję do macierzy.

Wyświetlenie

Mając dane w postaci macierzy, możemy je wyświetlić używając biblioteki matplotlib w następujący sposób:

import matplotlib.pyplot as plt

ext = [extent.xMinimum(), extent.xMaximum(), extent.yMinimum(), extent.yMaximum()]
nadata = raster.dataProvider().sourceNoDataValue(1)
data = np.ma.masked_values(array, nadata) # maskowanie wartości

plt.figure(figsize = (4, 4))
plt.imshow(data, extent = ext, cmap = "terrain")
plt.colorbar()
plt.title("Wysokość terenu [m n.p.m.]")
plt.xlabel("Długość geograficzna [m]")
plt.ylabel("Szerokość geograficzna [m]")
plt.show()

Algebra rastrów

Klasa QgsRasterCalculator umożliwia dostęp do kalkulatora rastrów w QGIS i wykorzystanie różnych operatorów uwzględniając matematyczne, trygonometryczne, porównania czy logiczne (dzięki czemu możliwe jest zamaskowanie czy reklasyfikacja wartości rastra). Potencjalne zastosowania obejmują:

  • uwypuklenie form rzeźby terenu wykorzystując współczynnik przewyższenia,
  • wyznaczenie dolin rzecznych wykorzystując maskowanie,
  • wyznaczenie klas terenu używając prostych reguł klasyfikacyjnych.

Pamiętaj, że w przypadku obliczeń uwzględniających kilka rastrów, to bezwzględnie muszą posiadać jednakowe zasięgi, rozdzielczość oraz układy współrzędnych.

from qgis.analysis import QgsRasterCalculator, QgsRasterCalculatorEntry

a = QgsRasterCalculatorEntry()
a.raster = raster
a.bandNumber = 1
a.ref = "DEM@1"

calculator = QgsRasterCalculator(
    "'DEM@1' * 3",                   # formulaString
    "output.tif",                    # outputFile
    "GTiff",                         # outputFormat
    raster.extent(),                 # outputExtent
    raster.crs(),                    # outputCrs
    raster.width(),                  # nOutputColumns
    raster.height(),                 # nOutputRows
    [a],                             # rasterEntries
    QgsCoordinateTransformContext()  # transformContext
)

status = calculator.processCalculation()

if status != 0:
    print("Błąd zapisu")
else: 
    print("OK")

Maskowanie

Aby stworzyć maskę wartości należy wykorzystać operator logiczny i ustalony warunek, np. wysokość terenu jest mniejsza niż 80 m n.p.m. W ten sposób otrzymamy raster binarny reprezentujący obszary, które spełniają zadany warunek (wartość 1; prawda) oraz nie spełniają (wartość 0; fałsz). Następnie, jeśli przemnożymy maskę przez wartości rastra, to otrzymamy oryginalne wartości dla obszarów spełniających warunek, w przeciwnym razie iloczyn operacji wynosi 0 (czyli komórka zostaje zamaskowana).

# maskowanie
formulaString = "('DEM@1' < 80) * 'DEM@1'"

Reklasyfikacja

Podobnie do maskowania działa reklasyfikacja wartości rastra. Jeśli chcemy zreklasyfikować raster na trzy klasy wysokości (np. niziny, wyżyny, góry), to możemy zrobić to w następujący sposób:

# reklasyfikacja
formulaString = "('DEM@1' < 94) * 1 + ('DEM@1' >= 94 AND 'DEM@1' <= 125) * 2 + ('DEM@1' > 125) * 3"

Obliczenie histogramu

Histogram umożliwia analizę rozkładu częstości wartości komórek poprzez zgrupowanie i zliczenie ich w określonych przedziałach (interwałach). Do stworzenia histogramu musimy zdefiniować liczbę przedziałów (wyznaczają w jaki sposób zostaną zgrupowane wartości) oraz kanał, dla którego ma zostać obliczony. Wymienione parametry przekazujemy do metody histogram() z dataProvider(). Wartości histogramu przechowywane są w atrybucie histogramVector.

from qgis.core import QgsRasterHistogram

provider = raster.dataProvider()

# oblicz histogram
bin_count = 100  # liczba przedziałów w histogramie
histogram = provider.histogram(1, bin_count)
vals = histogram.histogramVector

Do wizualizacji wykorzystamy wykres słupkowy, gdzie na osi Y znajdzie się częstość występowania wartości, a na osi X wysokość terenu. Aby wyświetlić prawidłowe etykiety na osi X musimy je wcześniej opracować, tj. najpierw podzielić zakres wartości na równe części (step), a następnie przypisać każdemu przedziałowi od 1 do 100 odpowiednią etykietę (labels).

# wyliczenie etykiet dla osi X
elev_max = provider.bandStatistics(1).maximumValue
elev_min = provider.bandStatistics(1).minimumValue
step = (elev_max - elev_min) / (bin_count - 1)
labels = [elev_min + i * step for i in range(bin_count)]
import matplotlib.pyplot as plt

plt.figure(figsize = (4, 2))
plt.bar(labels, vals, width = step)
plt.xlabel("Wysokość terenu [m n.p.m]")
plt.ylabel("Częstość")
plt.show()

Przepróbkowanie

Do przepróbkowania (resampling) rastra można wykorzystać dwa zaawansowane programy pochodzące z GDAL, tj. gdalwarp (używany do transformacji rastrów między różnymi układami odniesienia) oraz gdal_translate (używany do konwersji formatu i typów danych). Aby je zastosować należy najpierw zainicjować QGIS (bez tego nie będzie możliwe znalezienie modułów przetwarzania).

# https://docs.qgis.org/latest/en/docs/user_manual/processing_algs/gdal/rasterprojections.html#warp-reproject
from qgis import processing

params = {
    "INPUT": filepath,          # nazwa obiektu lub ścieżka do pliku
    "TARGET_CRS": "EPSG:2180",  # reprojekcja
    "TARGET_RESOLUTION": 1000,  # rozdzielczość
    "RESAMPLING": 0,            # numer metody
    "OUTPUT": "DEM_1000.tif"
}

processing.run("gdal:warpreproject", params)

Alternatywnie można użyć odwołań systemowych poprzez moduł Pythona subprocess do qgis_process.

import subprocess

params = [
    "--INPUT=" + filepath,
    "--TARGET_CRS=EPSG:2180",
    "--TARGET_RESOLUTION=1000",
    "--RESAMPLING=0",
    "--OUTPUT=DEM_1000.tif"
]
result = subprocess.run(["qgis_process", "run", "gdal:warpreproject"] + params,
                        capture_output = True, text = True)
print(result.stdout)

Tworzenie rastra

Tworzenie nowego rastra składa się z trzech zasadniczych kroków:

  1. Zainicjowanie pustego rastra na dysku (jeden lub wiele kanałów).
  2. Wczytanie wartości i ich zastąpienie w QgsRasterBlock.
  3. Aktualizacja (nadpisanie) wartości rastra w pliku.

Na początku stwórzmy macierz z wartościami od 1 do 100, którą zapiszemy na dysku.

data = np.arange(1, 101).reshape(10, 10)
print(data)
[[  1   2   3   4   5   6   7   8   9  10]
 [ 11  12  13  14  15  16  17  18  19  20]
 [ 21  22  23  24  25  26  27  28  29  30]
 [ 31  32  33  34  35  36  37  38  39  40]
 [ 41  42  43  44  45  46  47  48  49  50]
 [ 51  52  53  54  55  56  57  58  59  60]
 [ 61  62  63  64  65  66  67  68  69  70]
 [ 71  72  73  74  75  76  77  78  79  80]
 [ 81  82  83  84  85  86  87  88  89  90]
 [ 91  92  93  94  95  96  97  98  99 100]]

Do zainicjowania pustego rastra na dysku służy klasa QgsRasterFileWriter i metoda createOneBandRaster(), która wymaga wskazania typu danych, liczby kolumn i wierszy, zakresu przestrzennego oraz układu współrzędnych. W drugim kroku stosujemy metodę block, która zasięgiem obejmuje cały raster i iterujemy po wszystkich komórkach rastra nadając im odpowiednie wartości z wejściowej macierzy (metoda setValue()). Na końcu zapisujemy wartości oraz wyłączamy tryb edycji warstwy.

ext = QgsRectangle(347000, 454000, 393000, 494000)
crs = QgsCoordinateReferenceSystem("EPSG:2180")

# (1)
writer = QgsRasterFileWriter("test.tif")
provider = writer.createOneBandRaster(
    dataType = Qgis.DataType.Float32,
    width = 10,
    height = 10,
    extent = ext,
    crs = crs
)

# (2)
block = provider.block(1, provider.extent(), provider.xSize(), provider.ySize())
for x in range(block.width()):
    for y in range(block.height()):
        value = data[x][y]
        block.setValue(x, y, value)

# (3)
provider.writeBlock(block, band = 1)
provider.setEditable(False)

Zauważ, że dla danych rastrowych typ danych określamy z klasy Qgis.DataType, a nie QMetaType.Type, jak w przypadku danych wektorowych.

Łączenie rastrów

Do łączenia kilku jednokanałowych rastrów w jeden wielokanałowy można wykorzystać podejście zaprezentowane w poprzedniej sekcji dotyczącej tworzenia nowego rastra z tą różnicą, iż zamiast metody createOneBandRaster należy użyć createMultiBandRaster, w której każdy raster zostanie zapisany w pętli jako osobny kanał. Pamiętaj, że rastry muszą posiadać tę samą liczbę wierszy, kolumn, zakres oraz układ przestrzenny. W przeciwnym razie wymagane jest ich wcześniejsze dopasowanie (przepróbkowanie).

Uproszczone podejście opiera się o wykorzystanie dwóch narzędzi GDAL:

  1. gdal_merge – fizycznie łączy dane z wielu rastrów w jeden plik.
  2. gdalbuildvrt – tworzy wirtualny raster (plik z rozszerzeniem .vrt oparty o XML), który wyłącznie odwołuje się do wejściowych plików rastrowych bez ich kopiowania, tj. działa jako wskaźnik do oryginalnych plików.
# https://docs.qgis.org/latest/en/docs/user_manual/processing_algs/gdal/rastermiscellaneous.html#merge

params = {
    "INPUT": [filepath, filepath],
    "SEPARATE": True,
    "NODATA_OUTPUT": 9999,
    "OUTPUT": "Merged.tif"
}

processing.run("gdal:merge", params)
# https://docs.qgis.org/latest/en/docs/user_manual/processing_algs/gdal/rastermiscellaneous.html#build-virtual-raster

params = {
    "INPUT": [filepath, filepath],
    "SEPARATE": True,
    "OUTPUT": "Merged.vrt"
}

processing.run("gdal:buildvirtualraster", params)

Zadania

  1. W katalogu landsat znajdziesz kanał panchromatyczny (B08) w rozdzielczości 15 m oraz kompozycję RGB w rozdzielczości 30 m. Napisz funkcję służącą do pansharpeningu, tj. zwiększenia rozdzielczości przestrzennej kanałów spektralnych na podstawie kanału panchromatycznego. Do obliczeń wykorzystaj średnią arytmetyczną: wynik = (kanał + kanał panchromatyczny) / 2. Wynik operacji zapisz na dysku jako jeden raster składający się z trzech kanałów w wyższej rozdzielczości w formacie .tif. Wyniki pośrednie zapisz jako pliki tymczasowe używając modułu tempfile. Dla chętnych można zaimplementować inną wybraną metodę pansharpeningu, np. Brovey’a.