Praca z rastrami w Pythonie

Materiał przygotowany z pomocą AI

Wymagania wstępne

Wpis zakłada znajomość podstaw Pythona i NumPy w stopniu podstawowym — patrz Python - podstawy. Skupiamy się tu na danych rastrowych; pojęcia GIS spoza tego obszaru (wektory, układy współrzędnych) są wzmiankowane tylko tam, gdzie są niezbędne.

Środowisko do uruchomienia

Bloki kodu zostały zweryfikowane w Pythonie 3.12.12 z numpy 2.4.4, rasterio 1.5.0, xarray 2026.4.0, rioxarray 0.22.0, matplotlib 3.10.9. Pełna instrukcja odtworzenia środowiska: Środowisko Python dla przykładów.

Rastry to jeden z dwóch fundamentalnych modeli danych w GIS. Numeryczny model terenu (NMT), zdjęcie satelitarne, mapa pokrycia terenu, mapa opadu interpolowana ze stacji meteorologicznych — to wszystko rastry. W hydrologii i hydraulice rastry są wszechobecne: NMT jest punktem wyjścia do wyznaczania grnic zlewni, mapa pokrycia terenu wchodzi do modelu odpływu, a wyniki symulacji 2D też często są rastrami (głębokość, prędkość, zasięg).

Ten wpis przedstawia podstawy: czym raster jest, jakie są jego cechy, w jakich formatach się go zapisuje, jakie typy danych są w użyciu i jak wykonać podstawowe operacje w Pythonie. Przykłady budujemy wokół jednego pliku — schematycznego DEM zlewni 10×10, na którego bazie wyliczamy mapę pokrycia terenu, mapę nachyleń i maskę cieku.

Czym jest raster

Raster to dwuwymiarowa siatka (wiersze i kolumny) komórek (pikseli) o stałym rozmiarze, gdzie każda komórka przechowuje jedną lub więcej wartości. W modelu wektorowym świat opisujemy obiektami (punkt, linia, poligon) z atrybutami; w modelu rastrowym dzielimy obszar na prostokąty i każdemu przypisujemy wartość liczbową.

Najczęstsze przykłady:

NMT (Numeryczny Model Terenu)/DEM (Digital Elevation Model) — pojedyncza wartość wysokości na komórkę
Ortofotomapa — trzy lub cztery pasma (R, G, B, ewentualnie NIR)
Pokrycie terenu — kod klasy (np. 1=woda, 2=las, 3=zabudowa…)
Maska — wartość binarna (prawda/fałsz, 1/0)
Wynik symulacji 2D — głębokość, prędkość, itp.

Anatomia rastra

Każdy raster opisany jest przez kilka cech:

Wymiary — liczba wierszy (height) i kolumn (width); razem dają liczbę komórek
Pasma (bands) — pojedyncze pasmo (DEM) lub wielopasmowy (RGB ma 3, multispektralny Sentinel-2 ma ich kilkanaście). Konwencja indeksowania w GDAL/rasterio: pasma numerowane od 1 (nie od 0)
Rozdzielczość przestrzenna (pixel size) — ile metrów (lub stopni, minut, sekund) odpowiada jednej komórce; np. NMT generowany z LiDAR-u może mieć 1 m, a raster z akumulacją spływu z HydroSHEDS może mieć 15’
Zasięg (extent, bounding box) — prostokąt geograficzny, który raster pokrywa: (xmin, ymin, xmax, ymax)
Układ współrzędnych (CRS, Coordinate Reference System) — informacja, w jakim układzie odniesienia leżą współrzędne; w Polsce typowo EPSG:2180 (geodezyjny PUWG 1992) lub EPSG:4326 (geograficzny WGS84)
Oznaczenie braku danych — specjalna wartość oznaczająca “brak danych” (poza zlewnią, poza zasięgiem zdjęcia, brakujący odczyt). Ważne, żeby NoData nie kolidował z prawdziwymi wartościami w danych; typowo używa się -9999, 255, NaN (w zależności od typu danych, tj. int, uit, float oraz ich wielkości)

Raster vs wektor

Cecha	Raster	Wektor
Reprezentacja	siatka komórek o stałym rozmiarze	obiekty (punkty, linie, poligony)
Dyskretyzacja	regularna, narzucona z góry	dowolna, dopasowana do obiektu
Powierzchnia ciągła	naturalnie (NMT, opad)	trudniej (TIN, izolinie)
Topologia	brak (sąsiedztwo to indeksy macierzy)	jawna (sąsiedztwo, zawieranie)
Rozmiar pliku	rośnie z rozdzielczością i obszarem	rośnie z liczbą obiektów

W praktyce większość projektów GIS używa obu modeli — np. DEM (raster) + granica zlewni (wektor, poligon) + stacje pomiarowe (wektor, punkt) razem w jednym workflow.

Popularne formaty rastrów

Format pliku zapisuje wszystkie powyższe informacje (dane + metadane) na dysku. Wybór formatu ma znaczenie: wpływa na rozmiar pliku, szybkość odczytu, wsparcie kompresji i pasm, kompatybilność z narzędziami.

Format	Rozszerzenie	Zastosowanie	Cechy
GeoTIFF	`.tif`, `.tiff`	uniwersalny standard de facto	TIFF + georeferencja, wsparcie wszystkiego
NetCDF	`.nc`	klimatologia, meteorologia, oceanografia	wielowymiarowy (czas, głębokość, …), CF
HDF5	`.h5`, `.hdf5`	dane satelitarne	hierarchiczny, samoopisujący
ASCII Grid	`.asc`	wymiana z prostymi narzędziami	tekstowy — łatwy do podejrzenia, ogromny
JPEG2000	`.jp2`	duże zdjęcia satelitarne (Sentinel-2 L1C. )	silna kompresja, wolniejszy odczyt
PNG/JPG	`.png`/`.jpg`	tylko tam, gdzie nie potrzeba georeferencji	kompresja, brak metadanych przestrzennych

GeoTIFF — domyślny wybór

GeoTIFF to rozszerzenie formatu TIFF o tagi geograficzne (CRS, transform). To najczęściej spotykany format dla pojedynczego rastra: NMT, ortofoto, mapy pochodne. Wspiera dowolne typy danych, dowolną liczbę pasm, kompresję, NoData, piramidy.

Kompresje typowo używane w GeoTIFF:

Kompresja	Stratność	Charakterystyka
brak	—	szybki odczyt, ogromne pliki
LZW	bezstratna	dobry kompromis dla danych całkowitych
Deflate	bezstratna	popularna, podobna do LZW
ZSTD	bezstratna	nowsza, szybsza dekompresja, gęstsza
PackBits	bezstratna	prosta RLE, dobra dla masek
JPEG	stratna	tylko dla 8-bit RGB; nie używać dla DEM lub klasyfikacji

NetCDF i HDF5 — wielowymiarowość

Dla danych z osią czasu, głębokości lub innym dodatkowym wymiarem GeoTIFF nie wystarcza (jest 2D + pasma). Wówczas używa się NetCDF lub HDF5 — formatów hierarchicznych, które zapisują tablice n-wymiarowe wraz z opisem osi.

Przykład: pole opadu z modelu klimatycznego o wymiarach (czas: 365, lat: 100, lon: 100). Wszystkie te dane mieszczą się w jednym pliku NetCDF z osiami time, lat, lon i metadanymi opisującymi jednostki, długość kroku czasowego itd. W Pythonie czyta się je biblioteką xarray.

ASCII Grid — czytelny, ale ciężki

.asc to format tekstowy — można go otworzyć w edytorze i przeczytać. Nadaje się do prostych transferów lub debugowania, ale wartości zapisane jako tekst zajmują znacznie więcej miejsca niż w GeoTIFF (8-cyfrowa liczba zmiennoprzecinkowa to 8+ znaków/bajtów zamiast 4 dla float32). Dla produkcyjnej pracy nie zalecany.

Typy danych

Każda komórka rastra przechowuje wartość konkretnego typu — i ten typ determinuje, jakie wartości można w niej zapisać i ile bajtów zajmuje. To kluczowy parametr: pomyłka w typie danych psuje raster (overflow, utrata precyzji) lub niepotrzebnie zwiększa rozmiar pliku.

Tabela typów

Typ	Bajty/komórkę	Zakres	Typowe zastosowanie
`bool`	1*	`False` / `True`	maski binarne
`uint8`	1	0-255	klasyfikacje (do 255 klas), pasma RGB
`int8`	1	$- 128$ do $127$	rzadko (preferuje się `uint8` lub większe `int`)
`uint16`	2	0-65 535	dane satelitarne (Sentinel-2 BOA), DEM w cm
`int16`	2	$- 32768$ do $32767$	DEM globalny (SRTM), dane temperaturowe
`uint32`	4	$0$ - $\sim 4.29 \times 1 0^{9}$	identyfikatory, duże klasyfikacje
`int32`	4	$\pm 2.15 \times 1 0^{9}$	rzadziej (zwykle wystarczy `int16` lub `float32`)
`float32`	4	$\pm 3.4 \times 1 0^{38}$ , ~7 cyfr	DEM precyzyjny, NDVI, wyniki obliczeń
`float64`	8	$\pm 1.8 \times 1 0^{308}$ , ~15 cyfr	obliczenia wymagające dużej precyzji
`complex64`	8	dwa float32 (Re + Im)	dane SAR (radar) — surowe sygnały
`complex128`	16	dwa float64	jw., wyższa precyzja

* W NumPy bool zajmuje 1 bajt na element. W formatach binarnych (np. niektóre kompresje) możliwe upakowanie 8 wartości na bajt.

Co znaczy “uint” vs “int”

uint (unsigned integer) — bez znaku, tylko wartości nieujemne. uint8 mieści 0-255 (256 wartości w 1 bajcie)
int (signed integer) — ze znakiem, dodatnie i ujemne. int8 mieści -128 do 127 (też 256 wartości w 1 bajcie, ale przesunięte)

Co znaczy “float”

float32 — pojedyncza precyzja IEEE 754, 4 bajty, mantysa daje ~7 cyfr znaczących. Wystarcza do DEM, indeksów, większości danych ciągłych
float64 — podwójna precyzja, 8 bajtów, ~15 cyfr znaczących. Standard w obliczeniach naukowych w Pythonie/NumPy, ale dla rastra zwykle przesada — podwaja rozmiar pliku bez praktycznego zysku

Tylko float mieści NaN

Specjalna wartość NaN (Not a Number) istnieje wyłącznie w typach zmiennoprzecinkowych. W typach całkowitych nie ma jak zaznaczyć “brak danych” — trzeba użyć konkretnej wartości NoData (np. 255 dla uint8, -9999 dla int16) lub zapisać osobną maskę.

Pułapki typów danych

Overflow w typach całkowitych. W NumPy operacje na typach całkowitych nie rzucają wyjątku przy przepełnieniu — wartość “zawija się”:

import numpy as np
 
a = np.array([200, 200], dtype=np.uint8)
print(a + 100)
# [44 44]   <- 300 mod 256 = 44, bez ostrzeżenia (zwykle z RuntimeWarning)

Stąd przy operacjach na danych całkowitych warto najpierw rzutować na większy typ:

b = a.astype(np.int32) + 100
print(b)
# [300 300]

NoData a zakres typu. Wartość NoData musi mieścić się w typie i nie kolidować z prawdziwymi danymi. Dla uint8 (0-255) typowo używa się 255 jako NoData — pod warunkiem, że klasa “255” w danych niczego nie reprezentuje. Dla int16 typowo -9999. Dla float32 można użyć np.nan lub -9999.0 (zależy od konwencji).

Konwersja typu psuje wartości. Konwersja float32 -> uint8 obcina część ułamkową i zaokrągla do zakresu 0-255. Konwersja uint16 -> uint8 obcina górną połowę:

arr = np.array([0.7, 1.5, 300.0], dtype=np.float32)
print(arr.astype(np.uint8))
# [0 1 44]   ← 300 wyszło poza uint8, "zawinięte" do 44

Wpływ typu na rozmiar pliku

Dla rastra 10000 × 10000 (sto milionów komórek):

Typ	Rozmiar surowy (bez kompresji)
`uint8`	100 MB
`int16`	200 MB
`float32`	400 MB
`float64`	800 MB

Po zastosowaniu kompresji (LZW/Deflate) pliki bywają 2-10× mniejsze, zależnie od redundancji w danych. Mimo to dobór odpowiedniego typu jest pierwszym krokiem do utrzymania rozsądnych rozmiarów.

Reguły wyboru typu

Klasyfikacja do 255 klas → uint8
Klasyfikacja do 65 535 klas → uint16 (rzadko potrzebne)
Maska binarna → uint8 (0/1) lub bool; uint8 lepiej współpracuje z większością narzędzi
DEM lokalny w wysokiej rozdzielczości → float32 (typowy wybór; centymetrowa precyzja jest zachowana)
DEM globalny niskiej rozdzielczości → int16 w metrach (SRTM, GLO-30)
NDVI / inne indeksy $- 1 \dots 1$ → float32
Wynik modelu hydraulicznego (głębokość, prędkość) → float32
RGB/RGBA ortofotomapa → uint8 na pasmo

Reguła ogólna: wybierz najmniejszy typ, który mieści dane oraz wartość NoData, bez utraty istotnej precyzji.

Biblioteki Pythona do rastrów

Biblioteka	Rola	Kiedy używać
NumPy	tablice n-wymiarowe, operacje wektorowe	każda praca z rastrem (rastry są ndarray-ami)
rasterio	I/O rastrów + GDAL pod spodem	odczyt/zapis GeoTIFF, transform, CRS, maskowanie
xarray	tablice z etykietowanymi osiami	wielowymiarowe dane (NetCDF, czas, głębokość)
rioxarray	most xarray ←> rasterio	xarray z georeferencją, reprojekcie, geoprocessing
GDAL	niskopoziomowe API	rzadko bezpośrednio; preferuj rasterio
scipy.ndimage	filtry, morfologia, etykietowanie	wygładzanie, gradienty, regionalizacja
pyproj	transformacje współrzędnych	konwersje między CRS
whitebox / pysheds	hydrologia rastrowa	flow direction, akumulacja, delineacja zlewni
matplotlib	wizualizacja	szybkie podglądy w skryptach i notatnikach

W tym wpisie skupiamy się na trzech: NumPy jako fundamencie obliczeń, rasterio do I/O i metadanych przestrzennych, xarray + rioxarray dla wielowymiarowych danych i podejścia z etykietami osi.

Instalacja

Pełna instrukcja w Środowisko Python dla przykładów — wspólne requirements.txt z zafiksowanymi wersjami pokrywa też ten wpis (numpy, rasterio, xarray, rioxarray, matplotlib).

Skrót — pobranie i instalacja w dowolnym katalogu roboczym:

curl -O https://debever.pl/tools/requirements.txt
curl -o .python-version https://debever.pl/tools/python-version
python -m venv .venv
source .venv/bin/activate
pip install -r requirements.txt

Na macOS i Linuksie pip install rasterio zwykle działa bez dodatkowych kroków — wheel zawiera już skompilowaną bibliotekę GDAL. Na Windows czasem potrzeba użyć paczek z Christoph Gohlke’s repo lub conda-forge.

NumPy — fundament

Każdy raster w pamięci to NumPy ndarray. Wszystkie operacje (algebra, statystyki, filtry) wykonujemy na tym ndarray-u. Biblioteki rasterowe zwracają ndarray przy odczycie i przyjmują ndarray przy zapisie.

import numpy as np
 
# tworzenie rastra wprost
arr = np.zeros((100, 100), dtype=np.float32)
arr[40:60, 40:60] = 100.0  # kwadrat 20x20 z wartością 100
 
print(arr.shape)   # (100, 100)
print(arr.dtype)   # float32
print(arr.mean())  # 4.0  - 400 pikseli * 100 / 10000

rasterio — I/O i georeferencja

rasterio to “pythonowy” wrapper na GDAL. Dostarcza wzorzec context manager (z with), który otwiera plik, daje dostęp do danych i metadanych, a po wyjściu z bloku zamyka go automatycznie.

import rasterio
 
with rasterio.open("dem.tif") as src:
    print(src.width, src.height)   # liczba kolumn i rzędów
    print(src.count)               # liczba pasm
    print(src.dtypes)              # ('float32',)
    print(src.crs)                 # CRS, np. EPSG:2180
    print(src.transform)           # transformacja afiniczna
    print(src.bounds)              # (left, bottom, right, top)
    print(src.nodata)              # wartość NoData
 
    arr = src.read(1)              # odczyt 1. pasma jako ndarray

src.read() bez argumentu zwraca wszystkie pasma jako tablicę 3D (bands, height, width). src.read(1) zwraca tylko 1. pasmo jako 2D (height, width).

xarray i rioxarray — rastry z etykietowanymi osiami

xarray modeluje dane jako DataArray — tablicę z nazwami osi (x, y, time, band) i współrzędnymi (rzeczywiste wartości na każdej osi, nie indeksy 0..N). To bardzo wygodne przy większej liczbie wymiarów (np. dane czasowe) i przy operacjach typu “weź wartość dla time=2024-05-10”.

rioxarray dodaje do xarray akcesor .rio z metodami GIS (CRS, reprojection, clipping, zapis przez rasterio).

import rioxarray as rxr
 
# odczyt z rastra przez rioxarray (1. pasmo jako 2D)
da = rxr.open_rasterio("dem.tif", masked=True).squeeze()
 
print(da)               # DataArray z osiami x, y i CRS
print(da.rio.crs)       # CRS
print(da.rio.bounds())  # (left, bottom, right, top)
print(float(da.mean())) # statystyki - z automatycznym pominięciem NoData (masked=True)

masked=True zamienia NoData na NaN, dzięki czemu funkcje statystyczne (.mean(), .max()) automatycznie pomijają te komórki.

Tworzenie przykładowego rastra

Zbudujemy schematyczny DEM zlewni 10 × 10. Reguła:

dem(i, j) = 195 − 10·i + 8·|j − 5|

gdzie i to numer rzędu (0 = góra, 9 = dół), j to numer kolumny (5 to oś doliny, w której biegnie ciek). Wzór daje:

spadek N→S (każdy rząd niżej o 10 m, od 195 m w górnym rzędzie do 105 m w dolnym wzdłuż osi doliny)
V-dolinę symetrycznie wokół kolumny 5 (cele oddalone o |j-5| kolumn są wyżej o 8·|j-5| m niż dno doliny)

Wartości skrajne: minimum 105 m (komórka (9, 5) — ujście zlewni), maksimum 235 m (komórka (0, 0) — lewy-górny narożnik). Oś doliny w kolumnie 5 dzieli planszę asymetrycznie: po lewej jest 5 kolumn (j = 0..4), po prawej tylko 4 (j = 6..9), więc komórka (0, 9) ma “tylko” 227 m. Asymetrię widać też w pełnej macierzy poniżej.

import numpy as np
 
i, j = np.meshgrid(np.arange(10), np.arange(10), indexing="ij")
dem = (195.0 - 10.0 * i + 8.0 * np.abs(j - 5)).astype(np.float32)
 
print(dem.shape, dem.dtype)
print(f"min={float(dem.min())} max={float(dem.max())}")

Wynik:

(10, 10) float32
min=105.0 max=235.0

Pełna macierz (każdy wiersz to jeden rząd DEM, kolumny od zachodu do wschodu):

235 227 219 211 203 195 203 211 219 227
225 217 209 201 193 185 193 201 209 217
215 207 199 191 183 175 183 191 199 207
205 197 189 181 173 165 173 181 189 197
195 187 179 171 163 155 163 171 179 187
185 177 169 161 153 145 153 161 169 177
175 167 159 151 143 135 143 151 159 167
165 157 149 141 133 125 133 141 149 157
155 147 139 131 123 115 123 131 139 147
145 137 129 121 113 105 113 121 129 137

Dolinę widać “na pierwszy rzut oka”: w każdym wierszu minimum jest w kolumnie 5, wartości rosną symetrycznie ku krawędziom. Ciek (najniższe komórki w każdym rzędzie) tworzy pionową linię w środku.

Transformacja afiniczna i CRS

Macierz NumPy to dopiero połowa rastra. Druga połowa to transformacja afiniczna — informacja, jak (rząd, kolumna) przelicza się na (x, y) we współrzędnych geograficznych.

W rasterio transformację najprościej tworzy się funkcją from_origin:

from rasterio.transform import from_origin
 
pixel_size = 100.0  # 100 m × 100 m
transform = from_origin(west=500_000, north=400_000, xsize=pixel_size, ysize=pixel_size)
crs = "EPSG:2180"   # PUWG 1992
 
print(transform)

Transformacja jest reprezentowana jako:

| 100.00,   0.00, 500000.00 |
|   0.00,-100.00, 400000.00 |
|   0.00,   0.00,      1.00 |

Trzy kluczowe parametry: rozmiar piksela X (100 m), rozmiar piksela Y (-100 m, znak ujemny bo Y rośnie ku północy a wiersze rosną ku południowi), i współrzędne lewego-górnego rogu (500 000, 400 000) w EPSG:2180.

Lewy-górny róg, nie środek piksela

Konwencja w rasterio (i GDAL): transform odnosi się do lewego-górnego rogu lewego-górnego piksela, nie do jego środka. Inne konwencje (np. ESRI) bywają różne — sprawdzaj dokumentację.

Zapis DEM jako GeoTIFF

import rasterio
from rasterio.transform import from_origin
 
pixel_size = 100.0
transform = from_origin(west=500_000, north=400_000, xsize=pixel_size, ysize=pixel_size)
crs = "EPSG:2180"
 
with rasterio.open(
    "/tmp/dem.tif", "w",
    driver="GTiff",
    height=dem.shape[0],
    width=dem.shape[1],
    count=1,
    dtype=dem.dtype,
    crs=crs,
    transform=transform,
    nodata=-9999.0,
    compress="lzw",
) as dst:
    dst.write(dem, 1)
 
print("Zapisano /tmp/dem.tif")

Argumenty rasterio.open(..., "w", ...):

driver="GTiff" — GeoTIFF (najczęstszy wybór)
height, width — wymiary z macierzy NumPy
count — liczba pasm (tu 1)
dtype — typ danych (tu float32 z dem.dtype)
crs, transform — opisane wyżej
nodata — wartość NoData; -9999.0 to konwencja dla DEM (poza rzeczywistym zakresem wysokości)
compress="lzw" — bezstratna kompresja

Plik /tmp/dem.tif można teraz otworzyć w QGIS, ArcGIS, lub odczytać z powrotem w Pythonie.

Operacje na rastrach

Najczęstsze operacje na rastrach dzielą się na kilka kategorii:

Statystyki globalne i lokalne — min, max, średnia, percentyle, histogram
Algebra rastrowa — operacje arytmetyczne na pikselach
Reklasyfikacja — przekształcanie wartości w klasy lub inne wartości wg progów
Maskowanie — wybieranie komórek spełniających warunek; obsługa NoData
Filtrowanie / pochodne — gradient, średnia ruchoma
Resampling i reprojection — zmiana rozdzielczości lub układu współrzędnych
Wycinanie i mozaikowanie — przycinanie do poligonu, łączenie kilku rastrów

Statystyki

print(f"min: {float(dem.min())}")
print(f"max: {float(dem.max())}")
print(f"mean: {float(dem.mean()):.2f}")
print(f"std: {float(dem.std()):.2f}")
print(f"median: {float(np.median(dem)):.2f}")
print(f"p5, p95: {np.percentile(dem, [5, 95])}")

Dla rastrów z wartościami NoData używaj np.nanmin, np.nanmean (po wcześniejszej zamianie NoData na NaN) albo np.ma.masked_array:

# wariant 1: NaN
data = dem.copy()
data[dem == -9999] = np.nan
print(np.nanmean(data))
 
# wariant 2: masked array
masked = np.ma.masked_equal(dem, -9999)
print(masked.mean())

Algebra rastrowa

Operacje arytmetyczne działają na całych macierzach (broadcasting NumPy):

# konwersja jednostek (m -> km)
dem_km = dem / 1000.0
 
# odchylenie od średniej
deviation = dem - dem.mean()
 
# suma dwóch rastrów (wymagają tej samej geometrii!)
combined = dem + 0.0  # placeholder, zamiast tego np. rzeczywista mapa korekty

Operacje między dwoma rastrami wymagają, aby miały identyczne wymiary, transformację i CRS. Jeżeli się różnią — najpierw trzeba je wyrównać przez resampling/reprojection (sekcja niżej).

Reklasyfikacja — pokrycie terenu z DEM

Reklasyfikacja przekształca ciągłe wartości (wysokość) w klasy dyskretne (kategorie pokrycia terenu). Dla naszej zlewni przyjmijmy umowne progi:

Klasa	Kod	Zakres wysokości	Opis
woda	1	$h < 130$ m	dno doliny przy ujściu
łąka	2	$130 \leq h < 170$ m	dolina
las	3	$170 \leq h < 210$ m	zbocza
grunty wyżynne	4	$h \geq 210$ m	szczytowe partie

landcover = np.where(dem < 130, 1,
              np.where(dem < 170, 2,
                np.where(dem < 210, 3, 4))).astype(np.uint8)
 
print(landcover)
print("rozkład klas:", dict(zip(*np.unique(landcover, return_counts=True))))

Wynik:

[[4 4 4 4 3 3 3 4 4 4]
 [4 4 3 3 3 3 3 3 3 4]
 [4 3 3 3 3 3 3 3 3 3]
 [3 3 3 3 3 2 3 3 3 3]
 [3 3 3 3 2 2 2 3 3 3]
 [3 3 2 2 2 2 2 2 2 3]
 [3 2 2 2 2 2 2 2 2 2]
 [2 2 2 2 2 1 2 2 2 2]
 [2 2 2 2 1 1 1 2 2 2]
 [2 2 1 1 1 1 1 1 1 2]]
rozkład klas: {1: 11, 2: 39, 3: 39, 4: 11}

np.digitize — alternatywa

Dla wielu progów czytelniejsze bywa np.digitize(dem, bins=[130, 170, 210]) zamiast zagnieżdżonych np.where. Zwraca indeks przedziału (0, 1, 2 lub 3), z którego po przesunięciu +1 dostajemy te same klasy.

Maskowanie — wybór komórek

Maska to tablica bool o tych samych wymiarach co raster, z True tam, gdzie warunek jest spełniony.

# maska komórek dolinnych (poniżej średniej zlewni)
mean_h = dem.mean()
maska_dolina = dem < mean_h
print(f"komórek poniżej średniej ({mean_h:.0f} m): {maska_dolina.sum()}")
 
# wartości tylko dla doliny (reszta NaN)
dem_dolina = np.where(maska_dolina, dem, np.nan)

Maski mnoży się logicznie (&, |, ~):

# zbocza wschodnie poniżej 200 m
maska = (dem < 200) & (np.arange(10) > 5)

Pochodne DEM — gradient i nachylenie

np.gradient(arr, dy, dx) oblicza pochodną cząstkową w każdym kierunku. Dla DEM-u o pikselu 100 m:

# zwraca (dz/dy, dz/dx) - kolejność osi: rzędy potem kolumny
gy, gx = np.gradient(dem, 100.0)
 
# nachylenie jako tan(α)
slope_tan = np.sqrt(gy**2 + gx**2)
 
# nachylenie w stopniach
slope_deg = np.degrees(np.arctan(slope_tan))
 
print(f"slope min/max [°]: {float(slope_deg.min()):.2f} / {float(slope_deg.max()):.2f}")
print(f"slope mean [°]: {float(slope_deg.mean()):.2f}")

Dla naszej zlewni:

slope min/max [°]: 5.71 / 7.30
slope mean [°]: 7.14

Wartości są ledwie zróżnicowane: 5.71° wzdłuż osi doliny (gdzie tylko składowa y wnosi nachylenie, bo gradient w X w środku V wynosi 0) i 7.30° poza nią (gdzie do składowej y dochodzi ±0.08 w X). Średnia 7.14° odzwierciedla proporcję — tylko 10 z 100 komórek leży na osi doliny. W rzeczywistym DEM-ie nachylenia byłyby znacznie bardziej zmienne — od prawie zera w dolinach po dziesiątki stopni na zboczach.

To nie jest "prawdziwa" hydrologia

np.gradient daje matematyczny gradient, dobry do nachylenia i orientacji stoków. Do hydrologii rastrowej (kierunek spływu D8/D∞, akumulacja, delineacja zlewni) używaj wyspecjalizowanych bibliotek: pysheds, whitebox, GRASS GIS lub Whitebox Tools. Te narzędzia obsługują tzw. flow conditioning (wypełnianie zagłębień, breaching) i wieloskalowe algorytmy — czego prosty gradient NumPy nie zrobi.

Maska cieku — uproszczone wykrywanie

W naszej syntetycznej zlewni cieku nie wskazujemy ręcznie — możemy wykryć go jako “najniższą komórkę w każdym wierszu”. Ten trik działa, bo dolina została zaprojektowana symetrycznie z minimum w kolumnie 5:

stream_col = dem.argmin(axis=1)        # indeks kolumny minimum w każdym rzędzie
print("kolumny cieku:", stream_col)    # [5 5 5 5 5 5 5 5 5 5]
 
stream_mask = np.zeros_like(dem, dtype=bool)
stream_mask[np.arange(10), stream_col] = True
print(stream_mask.astype(np.uint8))

Wynik (1 = ciek, 0 = poza ciekiem):

[[0 0 0 0 0 1 0 0 0 0]
 [0 0 0 0 0 1 0 0 0 0]
 [0 0 0 0 0 1 0 0 0 0]
 [0 0 0 0 0 1 0 0 0 0]
 [0 0 0 0 0 1 0 0 0 0]
 [0 0 0 0 0 1 0 0 0 0]
 [0 0 0 0 0 1 0 0 0 0]
 [0 0 0 0 0 1 0 0 0 0]
 [0 0 0 0 0 1 0 0 0 0]
 [0 0 0 0 0 1 0 0 0 0]]

Dla rzeczywistego DEM-u ten trik nie zadziała — tam ciek wymaga prawdziwej analizy hydrologicznej (kierunki spływu + akumulacja lub maska lub wektor).

Resampling — zmiana rozdzielczości

Resampling to przeliczenie rastra na inną rozdzielczość. Najpopularniejsze metody:

Nearest neighbor — najbliższy sąsiad. Najszybszy, zachowuje wartości oryginalne. Stosowany do danych klasyfikacyjnych (kategorie nie powinny być uśredniane)
Bilinear — interpolacja liniowa po 4 sąsiadach. Łagodzi obraz, dobry dla danych ciągłych (DEM)
Bicubic — interpolacja kubiczna po 16 sąsiadach. Gładsza niż bilinear, ale może wprowadzać drobne artefakty
Average — średnia z grupy pikseli przy zmniejszaniu rozdzielczości

W rasterio resampling realizujemy przez rasterio.warp.reproject:

from rasterio.warp import reproject, Resampling
 
# zmniejszamy 10x10 -> 5x5 metodą "average"
dst = np.zeros((5, 5), dtype=np.float32)
reproject(
    source=dem,
    destination=dst,
    src_transform=transform,
    src_crs=crs,
    dst_transform=from_origin(500_000, 400_000, 200, 200),  # piksel 200 m
    dst_crs=crs,
    resampling=Resampling.average,
)
print(dst)

Nasze 10×10 z pikselem 100 m staje się 5×5 z pikselem 200 m. Dla danych klasyfikacyjnych (np. landcover) użylibyśmy Resampling.nearest.

Wycinanie do poligonu (clip)

W praktyce często chcemy ograniczyć raster do granic zlewni (poligonu wektorowego). W rasterio:

# pseudokod - wymaga poligonu jako GeoJSON-podobny dict lub geopandas.GeoSeries
# from rasterio.mask import mask
# with rasterio.open("dem.tif") as src:
#     out_image, out_transform = mask(src, [granica_zlewni], crop=True)

W rioxarray:

# da_clipped = da.rio.clip([granica_zlewni], crs="EPSG:2180")

Tu nie wykonujemy tego — wymagałoby wektorowego poligonu i biblioteki shapely/geopandas, co wykracza poza zakres wpisu.

Pełny przykład — DEM jako baza pochodnych

Łączymy elementy: budujemy DEM, derywujemy z niego trzy mapy (pokrycie terenu, nachylenie, maska cieku) i zapisujemy wszystkie cztery jako GeoTIFF z tym samym CRS i transformacją:

"""Schematyczna zlewnia 10x10: DEM + pochodne.
 
Skrypt:
1. Buduje syntetyczny DEM (V-dolina + spadek N-S).
2. Reklasyfikuje DEM na pokrycie terenu (4 klasy).
3. Wylicza nachylenie ze gradientu.
4. Wykrywa ciek jako najniższą komórkę w wierszu.
5. Zapisuje wszystkie cztery rastry jako GeoTIFF.
"""
 
import numpy as np
import rasterio
from rasterio.transform import from_origin
 
 
def build_dem() -> np.ndarray:
    """Zwraca syntetyczny DEM 10x10 jako float32.
 
    Reguła: dem(i, j) = 195 - 10*i + 8*|j-5|.
    """
    i, j = np.meshgrid(np.arange(10), np.arange(10), indexing="ij")
    return (195.0 - 10.0 * i + 8.0 * np.abs(j - 5)).astype(np.float32)
 
 
def classify_landcover(dem: np.ndarray) -> np.ndarray:
    """Reklasyfikuje DEM na klasy pokrycia terenu (uint8).
 
    1 = woda (h < 130 m), 2 = laka (130-170 m),
    3 = las (170-210 m), 4 = grunty wyzynne (h >= 210 m).
    """
    return np.where(dem < 130, 1,
             np.where(dem < 170, 2,
               np.where(dem < 210, 3, 4))).astype(np.uint8)
 
 
def compute_slope_deg(dem: np.ndarray, pixel_size: float) -> np.ndarray:
    """Oblicza nachylenie w stopniach z gradientu DEM."""
    gy, gx = np.gradient(dem, pixel_size)
    slope_tan = np.sqrt(gy**2 + gx**2)
    return np.degrees(np.arctan(slope_tan)).astype(np.float32)
 
 
def detect_stream(dem: np.ndarray) -> np.ndarray:
    """Maska cieku: najnizsza komorka w kazdym wierszu (uint8 0/1)."""
    stream_col = dem.argmin(axis=1)
    mask = np.zeros_like(dem, dtype=np.uint8)
    mask[np.arange(dem.shape[0]), stream_col] = 1
    return mask
 
 
def save_raster(
    data: np.ndarray,
    path: str,
    transform,
    crs: str,
    nodata,
) -> None:
    """Zapisuje 1-pasmowy raster jako GeoTIFF z kompresja LZW."""
    with rasterio.open(
        path, "w",
        driver="GTiff",
        height=data.shape[0],
        width=data.shape[1],
        count=1,
        dtype=data.dtype,
        crs=crs,
        transform=transform,
        nodata=nodata,
        compress="lzw",
    ) as dst:
        dst.write(data, 1)
 
 
def main() -> None:
    pixel_size = 100.0
    crs = "EPSG:2180"
    transform = from_origin(west=500_000, north=400_000,
                            xsize=pixel_size, ysize=pixel_size)
 
    dem = build_dem()
    landcover = classify_landcover(dem)
    slope = compute_slope_deg(dem, pixel_size)
    stream = detect_stream(dem)
 
    save_raster(dem,       "/tmp/dem.tif",       transform, crs, nodata=-9999.0)
    save_raster(landcover, "/tmp/landcover.tif", transform, crs, nodata=0)
    save_raster(slope,     "/tmp/slope.tif",     transform, crs, nodata=-9999.0)
    save_raster(stream,    "/tmp/stream.tif",    transform, crs, nodata=255)
 
    print(f"DEM:        min={float(dem.min())} max={float(dem.max())} mean={float(dem.mean()):.2f}")
    classes, counts = np.unique(landcover, return_counts=True)
    print(f"Pokrycie:   klasy={dict(zip(classes.tolist(), counts.tolist()))}")
    print(f"Nachylenie: min={float(slope.min()):.2f} max={float(slope.max()):.2f} sredia={float(slope.mean()):.2f} [stopnie]")
    print(f"Ciek:       komorek={int(stream.sum())}")
    print("Zapisano: /tmp/dem.tif, /tmp/landcover.tif, /tmp/slope.tif, /tmp/stream.tif")
 
 
if __name__ == "__main__":
    main()

Uruchomienie (python3 zlewnia.py) wypisuje:

DEM:        min=105.0 max=235.0 mean=170.00
Pokrycie:   klasy={1: 11, 2: 39, 3: 39, 4: 11}
Nachylenie: min=5.71 max=7.30 sredia=7.14 [stopnie]
Ciek:       komorek=10
Zapisano: /tmp/dem.tif, /tmp/landcover.tif, /tmp/slope.tif, /tmp/stream.tif

Odczyt z xarray + rioxarray

Plik zapisany przez rasterio możemy odczytać dowolną biblioteką. Z xarray dostajemy DataArray z osiami i CRS:

import rioxarray as rxr
 
dem_da = rxr.open_rasterio("/tmp/dem.tif", masked=True).squeeze()
 
print(dem_da)
print(f"CRS: {dem_da.rio.crs}")
print(f"Bounds: {dem_da.rio.bounds()}")
print(f"Pixel size: {dem_da.rio.resolution()}")
print(f"Mean (NoData ignored): {float(dem_da.mean()):.2f}")

masked=True automatycznie zamienia wartości NoData na NaN, dlatego .mean() pomija je przy liczeniu. Dla naszego DEM-u i tak NoData nie ma — wszystkie komórki są wypełnione.

Sprawdzenie spójności metadanych

Wszystkie cztery rastry mają tę samą geometrię (10×10, pikselami 100 m, ten sam CRS i transform). To krytyczne, jeśli chcemy je nakładać w QGIS lub łączyć w obliczeniach:

import rasterio
 
paths = ["/tmp/dem.tif", "/tmp/landcover.tif", "/tmp/slope.tif", "/tmp/stream.tif"]
for p in paths:
    with rasterio.open(p) as src:
        print(f"{p}: shape={src.shape} dtype={src.dtypes[0]} "
              f"crs={src.crs} bounds={tuple(round(b, 1) for b in src.bounds)}")

Wynik:

/tmp/dem.tif: shape=(10, 10) dtype=float32 crs=EPSG:2180 bounds=(500000.0, 399000.0, 501000.0, 400000.0)
/tmp/landcover.tif: shape=(10, 10) dtype=uint8 crs=EPSG:2180 bounds=(500000.0, 399000.0, 501000.0, 400000.0)
/tmp/slope.tif: shape=(10, 10) dtype=float32 crs=EPSG:2180 bounds=(500000.0, 399000.0, 501000.0, 400000.0)
/tmp/stream.tif: shape=(10, 10) dtype=uint8 crs=EPSG:2180 bounds=(500000.0, 399000.0, 501000.0, 400000.0)

Wszystkie cztery rastry pokrywają identyczny prostokąt (1000×1000 m) i można je nakładać piksel-na-piksel.

Dobre praktyki

Nie zmieniaj typu w pamięci bez potrzeby. arr.astype(np.float64) na całym rastrze niepotrzebnie zwiększa zużycie RAM. Pracuj w typie z pliku, konwertuj tylko gdy konieczne (np. operacje wymagające precyzji).
Zawsze ustawiaj NoData. Bez tej informacji narzędzia liczą statystyki łącznie z “pustymi” komórkami i mieszają je z prawdziwymi danymi.
Pasuj typ do zakresu danych. Pokrycie terenu (do 255 klas) zapisuj w uint8, nie int32 — mniejszy plik.
Używaj kompresji. compress="lzw" lub "deflate" daje kilkukrotnie mniejsze pliki bez utraty informacji.
Sprawdzaj geometrię (CRS, transform, shape) przy operacjach na wielu rastrach. Niezgodność prowadzi do trudnych do wykrycia błędów: dane “wyglądają OK”, ale piksele nie odpowiadają sobie.
Zamykaj pliki po zakończonej operacji. Używaj with rasterio.open(...) as src: zamiast ręcznego src.close() — context manager pilnuje zamknięcia nawet przy wyjątku.
Wizualizuj raster przed dalszym przetwarzaniem. matplotlib.pyplot.imshow(arr, cmap="terrain") to najszybszy sposób, żeby wykryć “raster do góry nogami”, brakujące dane lub nieoczekiwany zakres wartości.

Zobacz także

GIS — indeks GIS
Python - podstawy — fundamenty Pythona
Testowanie kodu w Pythonie — testy dla zlewnia.py z tego wpisu (refactor pod testowalność, fixtures, round-trip I/O, pełny zestaw testów)
Kartograf — narzędzie do pobierania NMT z GUGiK (rastry wysokościowe w EPSG:2180)
IMGWTools — pobieranie danych meteo z IMGW
Hydrolog — obliczenia hydrologiczne
rasterio docs — dokumentacja rasterio
xarray docs — dokumentacja xarray
rioxarray docs — most xarray ←> rasterio
GDAL — biblioteka pod spodem rasterio i wszystkich narzędzi rastrowych
Cloud Optimized GeoTIFF — specyfikacja COG
PROJ.org — biblioteka transformacji układów współrzędnych
pysheds — hydrologia rastrowa w Pythonie (flow direction, akumulacja, zlewnie)
Whitebox Tools — szerokie narzędzie do analiz GIS, dostępne też z poziomu Pythona

Baza wiedzy

Przeglądaj

Praca z rastrami w Pythonie

Praca z rastrami w Pythonie

Czym jest raster

Anatomia rastra

Raster vs wektor

Popularne formaty rastrów

GeoTIFF — domyślny wybór

NetCDF i HDF5 — wielowymiarowość

ASCII Grid — czytelny, ale ciężki

Typy danych

Tabela typów

Co znaczy “uint” vs “int”

Co znaczy “float”

Pułapki typów danych

Wpływ typu na rozmiar pliku

Reguły wyboru typu

Biblioteki Pythona do rastrów

Instalacja

NumPy — fundament

rasterio — I/O i georeferencja

xarray i rioxarray — rastry z etykietowanymi osiami

Tworzenie przykładowego rastra

Transformacja afiniczna i CRS

Zapis DEM jako GeoTIFF

Operacje na rastrach

Statystyki

Algebra rastrowa

Reklasyfikacja — pokrycie terenu z DEM

Maskowanie — wybór komórek

Pochodne DEM — gradient i nachylenie

Maska cieku — uproszczone wykrywanie

Resampling — zmiana rozdzielczości

Wycinanie do poligonu (clip)

Pełny przykład — DEM jako baza pochodnych

Odczyt z xarray + rioxarray

Sprawdzenie spójności metadanych

Dobre praktyki

Zobacz także

Graf

Spis treści

Odnośniki zwrotne