Wpis zakłada znajomość podstaw Pythona i NumPy w stopniu podstawowym — patrz Python - podstawy. Dla części przykładów (rastrowe I/O przez rioxarray) przyda się też tło z Praca z rastrami w Pythonie: czym jest CRS, transformacja afiniczna, NoData. Nie jest to wymóg twardy — sekcje, które tego dotyczą, są wyraźnie oznaczone.
Środowisko do uruchomienia
Bloki kodu w tym artykule zostały zweryfikowane w Pythonie 3.12.12 z xarray 2026.4.0, numpy 2.4.4, pandas 3.0.2, matplotlib 3.10.9, rioxarray 0.22.0, rasterio 1.5.0. Pełna instrukcja odtworzenia środowiska: Środowisko Python dla przykładów.
NumPy daje nam tablice n-wymiarowe — ndarray — ale tablica numpy nie wie co znaczą jej osie. Wymiar pierwszy to czas? Szerokość geograficzna? Numer pasma? Trzeba pamiętać konwencję i nie pomylić ich w kodzie. Pandas dokłada etykiety dla 1D (Series) i 2D (DataFrame), ale powyżej dwóch wymiarów się kończy. A dane meteorologiczne, klimatyczne i rastrowe naturalnie mieszkają w trzech, czterech lub pięciu wymiarach: czas × szerokość × długość × wysokość × pasmo.
xarray to biblioteka, która zachowuje moc obliczeniową NumPy, dokłada etykiety wymiarów i współrzędne na każdej osi i robi to w sposób skalowalny do wielu wymiarów. To naturalne narzędzie dla pól meteorologicznych, danych klimatycznych, wyników modeli oceanograficznych — ale też dla zwykłej, czasowej serii pomiarów z wielu stacji.
Wpis omawia: czym xarray jest i czym nie jest, dwie centralne struktury danych (DataArray, Dataset), wczytywanie i zapis (NetCDF, GeoTIFF), indeksowanie po etykietach, arytmetykę z broadcastingiem, agregacje, operacje na osi czasu, leniwe obliczenia z dask — a na końcu pełne studium przypadku syntetycznego pola opadów nad Polską z agregacją miesięczną, oknem kroczącym i wizualizacją.
Czym jest xarray
W jednym zdaniu: xarray rozszerza NumPy o nazwy wymiarów i etykiety na każdej osi. Tablica numpy (365, 20, 20) to ślepe trzy osie. Ta sama tablica w xarray to (time, lat, lon) — z datą na osi time, szerokością geograficzną na lat, długością na lon — i samodokumentującą się reprezentacją w REPL-u.
Motywacja
Wyobraź sobie pole dziennych opadów dla roku 2025 nad obszarem 20 × 20 punktów siatki. Surowo w NumPy:
import numpy as npprecip = np.random.default_rng(42).gamma(1.0, 3.0, size=(365, 20, 20))print(precip.shape)# (365, 20, 20)# Suma opadu w lipcu, wszystkie punkty:suma_lipiec = precip[181:212].sum(axis=0) # ← osi 0 to czas, trzeba pamietacprint(suma_lipiec.shape)# (20, 20)
Problemy:
181:212 to lipiec? Trzeba przeliczyć dni od początku roku. Albo zapomnieć i pomylić sierpień z lipcem.
axis=0 — co to było? Czas, lat czy lon? Bez komentarza nie wiadomo.
Punkt (10, 10) ma jaką szerokość/długość? Bez osobnej tablicy lats i lons w zmiennych obok — nie wiadomo.
W xarray te same dane:
import xarray as xrimport pandas as pdimport numpy as nprng = np.random.default_rng(42)times = pd.date_range("2025-01-01", periods=365, freq="D")lats = np.linspace(49.0, 55.0, 20)lons = np.linspace(14.0, 24.5, 20)precip = xr.DataArray( rng.gamma(1.0, 3.0, size=(365, 20, 20)), dims=("time", "lat", "lon"), coords={"time": times, "lat": lats, "lon": lons}, name="precipitation", attrs={"units": "mm/day"},)# Suma opadu w lipcu — po etykiecie:suma_lipiec = precip.sel(time=slice("2025-07-01", "2025-07-31")).sum(dim="time")print(suma_lipiec.shape)# (20, 20)
Kod jest dłuższy na starcie, ale:
time=slice("2025-07-01", "2025-07-31") — czytelne na pierwszy rzut oka, niezależnie od kalendarza.
dim="time" — nie ma pomyłki, którą oś agregujemy.
Każdy punkt zna swoje współrzędne geograficzne (lat, lon).
attrs={"units": "mm/day"} — dane się same dokumentują.
Dla danych jednorazowych, jednowymiarowych, ten narzut się nie opłaca. Dla danych o większej liczbie wymiarów lub takich, które żyją w potokach przetwarzania przez dłuższy czas — bardzo szybko zwraca się sam.
xarray vs NumPy vs pandas
Cecha
NumPy ndarray
pandas DataFrame
xarray DataArray/Dataset
Wymiarów
dowolnie wiele
2 (indeks + kolumny)
dowolnie wiele
Etykietowane osie
nie
tak (2)
tak (wszystkie)
Heterogeniczne kolumny
nie
tak
nie (w obrębie DataArray) — Dataset ma wiele zmiennych
I/O CSV
przez pandas
natywne
przez pandas / konwersję
I/O NetCDF
nie
nie
natywne
I/O GeoTIFF
przez rasterio
nie
przez rioxarray
Etykiety czasu
nie
tak (DatetimeIndex)
tak
Broadcasting po nazwach
po pozycji
częściowo
po nazwie wymiaru
Skala (lazy/równolegle)
ograniczona
ograniczona
natywne wsparcie dask
Reguła kciuka:
2D, heterogeniczne typy w kolumnach (np. tabela: nazwa stacji jako łańcuch znaków, przepływ jako float, data jako datetime) → pandas.
Dane N-wymiarowe o jednorodnym typie (raster, pole klimatyczne, kostka czasowa) → xarray.
Surowe obliczenia bez metadanych → numpy (xarray pod spodem i tak używa numpy).
Dla kogo i kiedy
Typowi użytkownicy xarray:
klimatolodzy, meteorologowie, oceanografowie — naturalna domena (pliki NetCDF są tu standardem);
analitycy danych satelitarnych (kostki czasowe Sentinel-2/Landsat);
hydrolodzy pracujący z polami opadu, ewapotranspiracji, temperatury;
każdy, kto przetwarza wynik modelu numerycznego o wymiarach (czas, x, y, [z]).
Jeśli pracujesz wyłącznie z pojedynczymi rastrami 2D bez osi czasu — rasterio + numpy w zupełności wystarczą, a xarray to dodatkowy poziom abstrakcji bez wielkiego zwrotu. Jeśli dochodzi czas — xarray zwraca się natychmiast.
DataArray
DataArray to centralny obiekt biblioteki. Trzymaj w głowie obraz:
DataArray
├── data ← NumPy ndarray (lub dask array)
├── dims ← krotka nazw wymiarów: ("time", "lat", "lon")
├── coords ← słownik etykiet na osiach: time → ..., lat → ..., lon → ...
├── name ← nazwa zmiennej (opcjonalna)
└── attrs ← słownik metadanych (units, long_name, ...)
Najprostsza konstrukcja
import xarray as xrimport numpy as np# DataArray z czysto NumPy-owej tablicy 1D — bez wspolrzednycha = xr.DataArray([1, 2, 3])print(a)
Output:
<xarray.DataArray (dim_0: 3)> Size: 24Barray([1, 2, 3])Dimensions without coordinates: dim_0
xarray automatycznie nazwał wymiar dim_0 (bo nie podaliśmy nazwy) i nie ma żadnych współrzędnych. To minimum.
W reprezentacji obiektu widać wszystko: kształt, dtype, wymiary, współrzędne i jednostki. Nie trzeba print(arr.shape); print(arr.dtype); print(times); print(sites) jak w NumPy — jeden print(da) pokazuje pełen kontekst.
dims, coords i name — formalnie
dims — krotka nazw wymiarów; uporządkowana (kolejność ma znaczenie dla .values i operacji typu transpose).
coords — obiekt słownikopodobny zawierający etykiety na osiach. Każda współrzędna ma własny wymiar lub jest „skalarna” (jedna wartość bez wymiaru).
name — opcjonalna nazwa zmiennej; jeśli jest ustawiona, jest wyświetlana w reprezentacji obiektu i propaguje się przy konwersjach (np. do pandas.Series.name).
attrs — słownik z metadanymi; przechowuje konwencję CF (units, long_name, standard_name, …) lub własne pola.
Coords ze swoim wymiarem vs bez wymiaru
Współrzędna może mieć ten sam wymiar co dane (np. time na osi time) — wtedy jest indeksowalna (sel(time=...)). Może też być skalarna lub mieć osobny wymiar. Skalarnych coords używa się np. do oznaczenia, że cała tablica reprezentuje pomiar z jednej stacji: coords={"station_id": "WW01"}.
Praktyka — odczytywanie metadanych
W przykładach z plikami NetCDF (klimatologia) typowe atrybuty to:
import xarray as xrimport numpy as npda = xr.DataArray( np.zeros((3,)), dims="time", coords={"time": [0, 1, 2]}, name="precipitation", attrs={ "units": "kg m-2 s-1", "long_name": "Precipitation flux", "standard_name": "precipitation_flux", "_FillValue": -9999.0, },)print(da.attrs["units"]) # kg m-2 s-1print(da.attrs.get("units")) # to samo, bezpieczniej (None jak brak)print(da.attrs.get("missing")) # None — brak klucza
Atrybut units to konwencja CF (Climate and Forecast metadata conventions) — patrz oficjalna dokumentacja CF. xarray nie wymusza tej konwencji (każdy atrybut to po prostu wpis w słowniku), ale narzędzia w ekosystemie zakładają, że jednostki znajdują się w units.
Dataset
Dataset to kolekcja DataArray-jów dzielących te same wymiary i współrzędne. Analogicznie do pandas: DataArray ↔ Series, Dataset ↔ DataFrame.
Po co? W praktyce dane fizyczne często idą w komplecie: temperatura + opad + ciśnienie + wilgotność, wszystkie po tej samej siatce (time, lat, lon). Zamiast czterech osobnych obiektów DataArray trzymanych w słowniku w głowie, wszystko jest w jednym obiekcie Dataset.
xarray natywnie obsługuje kilka formatów wejścia/wyjścia. Bez dodatkowych pakietów (czyli „surowy” xarray) jest ograniczony — produkcyjnie używa się go z netCDF4, h5netcdf, zarr, rioxarray. W tym artykule pokażemy NetCDF (klasyczny standard) i GeoTIFF (przez rioxarray).
Wymagane backendy
Sam pakiet xarraynie umie zapisać NetCDF, Zarr czy GeoTIFF — potrzebuje silnika. W środowisku requirements.txt z Środowisko Python dla przykładów zainstalowane są rasterio i rioxarray (GeoTIFF). NetCDF wymaga dodatkowo netCDF4, h5netcdf lub scipy (każdy z nich potrafi obsłużyć NetCDF, z drobnymi różnicami). Sprawdzić, co jest dostępne:
import xarray as xrprint(xr.backends.list_engines())
Formaty
Format
Rozszerzenie
Charakterystyka
Pakiet
NetCDF4
.nc
binarny, hierarchiczny, samoopisujący; standard klimat/meteo
netCDF4 lub h5netcdf
NetCDF3
.nc
starsza wersja, ograniczenia, ale bardzo szeroka kompatybilność
scipy.io.netcdf
Zarr
.zarr/
format chunkowany, dobry do chmury, kompatybilny z dask
zarr
GeoTIFF
.tif
jednowymiarowy raster (+pasma), georeferencja
rioxarray + rasterio
HDF5
.h5
hierarchiczny, używany w danych satelitarnych
h5py lub netCDF4
NetCDF — zapis i odczyt
NetCDF (Network Common Data Format) to binarny, samoopisujący format hierarchiczny rozwijany przez UCAR/Unidata. Zapisuje tablice n-wymiarowe wraz z opisem osi, wymiarów, jednostek i innych metadanych. To format pierwszego wyboru dla danych klimatycznych, meteorologicznych i oceanograficznych. Pliki NetCDF respektują konwencje CF (Climate and Forecast), które standaryzują nazwy atrybutów (units, standard_name, axis, …) — dzięki temu plik z jednego ośrodka jest czytelny dla narzędzi z drugiego.
Zapis (zakładając, że w środowisku jest jeden z backendów netCDF4/h5netcdf/scipy):
import xarray as xrimport numpy as npimport pandas as pdrng = np.random.default_rng(42)times = pd.date_range("2025-01-01", periods=10)ds = xr.Dataset( data_vars={ "temp": (("time",), rng.normal(10, 5, 10)), "rain": (("time",), rng.gamma(1.0, 2.0, 10)), }, coords={"time": times}, attrs={"title": "Przykładowe dane", "Conventions": "CF-1.8"},)# Zapis — wymaga zainstalowanego backendu NetCDFds.to_netcdf("/tmp/dane.nc")# Odczytds_read = xr.open_dataset("/tmp/dane.nc")print(ds_read)ds_read.close()
Środowisko bazy nie ma backendu NetCDF
W aktualnym requirements.txt (sesja 2026-05-12) brakuje netCDF4/h5netcdf/scipy. Wywołanie ds.to_netcdf(...) rzuci ValueError: cannot write NetCDF files because none of the suitable backend libraries (netCDF4, h5netcdf, scipy) are installed. Aby uruchomić powyższy blok lokalnie, doinstaluj jeden z backendów:
pip install netCDF4# albo: pip install h5netcdf
Otwarcie pliku — open_dataset vs load_dataset
# open_dataset — leniwy (lazy); plik trzymany jest otwarty, dane czytane sa na zadanieds = xr.open_dataset("/tmp/dane.nc")# load_dataset — czyta WSZYSTKO do pamieci i zamyka plikds_full = xr.load_dataset("/tmp/dane.nc")
Reguła:
open_dataset — dla plików większych niż pamięć, gdy chcemy pracować na podzbiorze.
load_dataset — dla małych plików, gdy potem manipulujemy strukturą katalogów.
Po open_dataset(...) zamykaj przez .close() albo użyj menedżera kontekstu with xr.open_dataset(...) as ds:.
GeoTIFF przez rioxarray
GeoTIFF nie jest „rodzimym” formatem xarray (w pliku jest tylko 2D raster z opcjonalnymi pasmami), ale rioxarray dorzuca odpowiednie nakładki. Sposób — patrz też Praca z rastrami w Pythonie.
rioxarray wymaga osi nazwanych x i y (konwencja zgodna z GeoTIFF). Akcesor .rio udostępnia metody specyficzne dla rastra: .write_crs(), .to_raster(), .bounds(), .clip(), .reproject().
Konwersja do/z pandas
import xarray as xrimport pandas as pdimport numpy as nptimes = pd.date_range("2025-01-01", periods=5)da = xr.DataArray( np.arange(15).reshape(5, 3), dims=("time", "site"), coords={"time": times, "site": ["A", "B", "C"]}, name="val",)# to_pandas — dla 1D zwroci Series, dla 2D DataFramedf = da.to_pandas()print(type(df).__name__)# DataFrameprint(df.head(3))
Output:
site A B Ctime2025-01-01 0 1 22025-01-02 3 4 52025-01-03 6 7 8
to_dataframe daje „długą” reprezentację (long form):
print(da.to_dataframe().head(6))
Output:
valtime site2025-01-01 A 0 B 1 C 22025-01-02 A 3 B 4 C 5
<xarray.Dataset> Size: 112BDimensions: (time: 4, site: 2)Coordinates: * time (time) datetime64[us] 32B 2025-01-01 2025-01-02 ... * site (site) object 16B 'A' 'B'Data variables: value (time, site) float64 64B 1.0 nan 2.0 nan nan 3.0 nan 4.0
Uwaga: brakujące kombinacje (np. stacja B w dniu 2025-01-01) zostają wypełnione wartością NaN. To naturalne, bo xarray zakłada „pełną siatkę” wymiarów.
Indeksowanie i selekcja
Indeksowanie to miejsce, w którym xarray najmocniej odstaje od NumPy — i najmocniej zarabia na sobie. Dwie podstawowe metody:
.isel(...) — integer-based, jak arr[i, j] w NumPy; indeks pozycyjny 0..N.
.sel(...) — label-based, po etykietach współrzędnych: sel(time="2025-07-01").
Przykładowe dane
import xarray as xrimport numpy as npimport pandas as pdtimes = pd.date_range("2025-01-01", periods=10, freq="D")da = xr.DataArray( np.arange(30).reshape(10, 3), dims=("time", "site"), coords={"time": times, "site": ["A", "B", "C"]}, name="val",)
isel — po pozycji
# Pojedynczy element — przeciecie indeksow:print(da.isel(time=0, site=0).values)# 0# Wycinek wzdluz osi czasu — pierwsze 3 dni:print(da.isel(time=slice(0, 3)).shape)# (3, 3)# Ostatni element:print(da.isel(time=-1).values)# [27 28 29]
isel jest blisko czystego NumPy — jeśli znasz arr[0:3], arr[-1], to też wiesz, jak działa isel(x=slice(0,3)) i isel(x=-1).
slice w sel jest obustronnie domknięty — inaczej niż w czystym Pythonie. To wynika z faktu, że indeks jest etykietowy, a nie pozycyjny: użytkownik mówi „od 2-go do 5-go”, co intuicyjnie obejmuje oba krańce.
Tolerancja — method="nearest" i tolerance
# Najbliższa data — przy zapytaniu o "2025-01-03 13:00" wraca dzien 04print(da.sel(time="2025-01-03 13:00", method="nearest").values)# [ 9 10 11]# Z ograniczeniem maksymalnej odleglosci (24 h):print(da.sel(time="2025-01-03 13:00", method="nearest", tolerance=pd.Timedelta("24h")).values)# [ 9 10 11]
Opcji method="nearest" używa się głównie do dopasowania niedokładnych dat lub współrzędnych geograficznych (np. „najbliższa stacja do punktu 51.5°N, 17.2°E”).
.loc[] — jeszcze inny sposób
Notacja loc[] (znana z pandas) działa też w xarray, używając etykiet:
# Pojedynczy element po etykietach:print(da.loc["2025-01-03", "B"].values)# 7# Wycinek:print(da.loc["2025-01-02":"2025-01-04", ["A", "C"]].values)
Output:
[[ 3 5] [ 6 8] [ 9 11]]
loc[] jest zwięzła, ale w mojej praktyce mniej czytelna niż sel(time=..., site=...). To kwestia gustu.
Indeksowanie warunkowe (boolean indexing)
# Wszystkie dni od 5. do konca:mask = da.time.dt.day >= 5print(da.sel(time=mask).shape)# (6, 3)# Wartosci wieksze niz 15:print(da.where(da > 15).values[:3])
Output:
[[nan nan nan] [nan nan nan] [nan nan nan]]
.where(cond) zwraca taki sam kształt z NaN tam, gdzie warunek nie jest spełniony. Drugi argument pozwala podać wartość zastępczą:
print(da.where(da > 15, 0).values[:3])
Output:
[[0 0 0] [0 0 0] [0 0 0]]
Indeksowanie przez listy etykiet (fancy indexing)
# Dwie wybrane stacje:print(da.sel(site=["A", "C"]).shape)# (10, 2)# Trzy wybrane dni:print(da.sel(time=["2025-01-01", "2025-01-05", "2025-01-10"]).shape)# (3, 3)
Reindex — zmiana siatki bez interpolacji
import xarray as xrimport numpy as npa = xr.DataArray([1.0, 2.0, 4.0, 8.0, 16.0], dims="x", coords={"x": [0, 10, 20, 30, 40]})# Reindex do nowej siatki — w nowych pozycjach wstawia NaN:b = a.reindex(x=[0, 5, 10, 15, 20])print(b.values)# [ 1. nan 2. nan 4.]# Z fill_value:c = a.reindex(x=[0, 5, 10, 15, 20], fill_value=-1.0)print(c.values)# [ 1. -1. 2. -1. 4.]# Z method="nearest" (najblizsza wartosc z oryginalu):d = a.reindex(x=[0, 5, 10, 15, 20], method="nearest")print(d.values)# [1. 2. 2. 4. 4.]
reindex jest „przepisywaniem na nową siatkę” bez interpolacji. Dla prawdziwej interpolacji liniowej lub typu spline służy .interp() — ale ta wymaga scipy.
Operacje na danych
Arytmetyka i broadcasting
NumPy stosuje broadcasting po kształcie (ostatnie wymiary muszą się zgadzać). xarray stosuje go po nazwie wymiaru — co jest bezpieczniejsze:
import xarray as xrimport numpy as npa = xr.DataArray(np.array([[1, 2, 3], [4, 5, 6]]), dims=("x", "y"), coords={"x": [10, 20], "y": [100, 200, 300]})# Wektor po osi y:b = xr.DataArray([10, 20, 30], dims="y", coords={"y": [100, 200, 300]})print((a + b).values)
Output:
[[11 22 33] [14 25 36]]
xarray rozpoznał, że b ma wymiar y i rozszerzył go na wymiar x z a. Nie trzeba b.reshape(1, 3) ani b[None, :].
Automatyczne wyrównanie po coords
Gdy dwie tablice mają wspólny wymiar, ale różne współrzędne, xarray automatycznie znajduje przecięcie:
To zachowanie chroni przed dodaniem do siebie danych z dwóch różnych dat lub punktów siatki. W NumPy np.array([1,2,3]) + np.array([100,200,300]) „udałoby się” bez ostrzeżenia, mimo że dane mogą reprezentować różne lokalizacje.
Agregacje
import xarray as xrimport numpy as npa = xr.DataArray(np.array([[1, 2, 3], [4, 5, 6]]), dims=("x", "y"), coords={"x": [10, 20], "y": [100, 200, 300]})print(float(a.mean())) # globalna srednia# 3.5print(a.mean(dim="x").values) # srednia po x — wynik wzdluz y# [2.5 3.5 4.5]print(a.mean(dim="y").values) # srednia po y — wynik wzdluz x# [2. 5.]print(int(a.sum())) # 21print(float(a.std())) # 1.707825127659933print(int(a.min()), int(a.max())) # 1 6print(float(a.quantile(0.5))) # 3.5
Dostępne agregacje: mean, sum, std, var, min, max, median, quantile, prod, count, argmin, argmax. Każda przyjmuje dim= (lub krotkę dim).
Liczenie tylko wartości (count) i NaN
import xarray as xrimport numpy as nparr = xr.DataArray([1.0, np.nan, 3.0, 4.0], dims="x")print(int(arr.count())) # 3 — liczba wartosci nie-NaNprint(float(arr.mean())) # 2.666... — domyslnie skipna=Trueprint(float(arr.mean(skipna=False))) # nan
Większość agregacji domyślnie pomija NaN (skipna=True). Dla typów bool lub int przyjmowane jest skipna=False — wtedy nie ma czego pomijać.
Uniwersalne funkcje (ufunc) z NumPy działają bezpośrednio:
import xarray as xrimport numpy as nparr = xr.DataArray([1.0, 4.0, 9.0, 16.0], dims="x")print(np.sqrt(arr).values)# [1. 2. 3. 4.]print(np.log(arr).values[:2])# [0. 1.38629436]
Wynik zachowuje typ DataArray z osiami — niezależnie od tego, że woła go funkcja z NumPy.
Atrybuty po operacjach — uważnie
Zachowanie attrs przy operacjach to częste źródło niespójności. W xarray 2026.4.0 domyślne zachowanie jest takie: pojedyncza operacja (typu * 2, .mean(), .isel(...)) zachowujeattrs; operacja na dwóch obiektach z różnymi atrybutami daje pusty attrs:
Globalnie można to zmienić przez xr.set_options(keep_attrs=True/False):
with xr.set_options(keep_attrs=False): print((a * 2).attrs) # {}
Zachowanie zależne od wersji
Historycznie xarray miał domyślnie keep_attrs=False — atrybuty były kasowane po arytmetyce. W nowszych wersjach zachowanie się zmieniło. Jeśli kod migruje ze starszej wersji, sprawdź ten szczegół eksperymentalnie albo wymuś jawnie przez xr.set_options.
Akcesor .str to przeniesienie odpowiednika z pandas — przydatny do filtrowania po nazwach stacji, identyfikatorach itp.
Wymiar czasowy
To miejsce, w którym xarray (razem z pandas pod spodem) świeci najjaśniej. Dane w hydrologii i meteorologii to przede wszystkim szeregi czasowe — dobowe, godzinne, 5-minutowe — a operacje typu „przelicz dzienne sumy na miesięczne” są stałym elementem przepływu pracy.
DatetimeIndex w coords
Wystarczy, że współrzędna na osi time jest typu datetime64 (najprościej: pandas.date_range), żeby xarray rozpoznał oś jako czasową i udostępnił akcesor .dt:
Akcesor .dt udostępnia standardowy zestaw z pandas: year, month, day, hour, dayofyear, weekday, quarter, season (DJF/MAM/JJA/SON), is_month_start, …
resample — z jednego kroku do drugiego
import xarray as xrimport pandas as pdimport numpy as nprng = np.random.default_rng(42)times = pd.date_range("2025-01-01", periods=100, freq="D")da = xr.DataArray(rng.gamma(1.0, 1.0, 100), dims="time", coords={"time": times}, name="precip")# Z dziennych do miesiecznych — suma:miesieczne_sumy = da.resample(time="MS").sum()print(miesieczne_sumy.shape)# (4,)print(miesieczne_sumy.values)# [28.38107635 20.87407208 33.05586341 7.59437749]# Z dziennych do tygodniowych — srednia:tygodniowe_srednie = da.resample(time="W").mean()print(tygodniowe_srednie.shape)# (15,)
Częste „kody częstotliwości” (pandas/pandas 3.0):
Kod
Znaczenie
D
dziennie
h
godzinnie (uwaga: małe h w pandas 3.0)
W
tygodniowo (z domyślnym końcem niedziela)
MS
początek miesiąca (Month Start)
ME
koniec miesiąca (Month End)
YS
początek roku (Year Start)
YE
koniec roku (Year End)
QS
początek kwartału
Zmiana kodów w pandas 3.x
Starsze pandasy używały H (godziny), M (Month End), Y (Year End). W pandas 3.0 (które idzie w requirements.txt) część z nich jest „przepisana”: małe h dla godzin, ME/YE z literą E na końcu miesiąca/roku. Stare kody w xarray rzucają FutureWarning albo błąd. Jeśli kopiujesz przykład z internetu sprzed roku — sprawdź ten szczegół.
Dostępne agregacje po resample: mean, sum, max, min, std, count, median, first, last. Można też wywołać własną funkcję przez .reduce(np.func, dim="time") lub .map(func).
groupby — agregacja po atrybucie
groupby agreguje wzdłuż dowolnej współrzędnej lub jej pochodnej:
import xarray as xrimport pandas as pdimport numpy as nprng = np.random.default_rng(42)times = pd.date_range("2025-01-01", periods=365, freq="D")da = xr.DataArray(rng.normal(10, 5, 365), dims="time", coords={"time": times})# Średnia po sezonach (DJF/MAM/JJA/SON):season_means = da.groupby("time.season").mean()print(season_means.values)# np.: [DJF=..., MAM=..., JJA=..., SON=...]print(season_means.coords["season"].values)# ['DJF' 'JJA' 'MAM' 'SON']# Średnia po miesiącu kalendarzowym:month_means = da.groupby("time.month").mean()print(month_means.shape)# (12,)
"time.season", "time.month", "time.year", "time.dayofyear" to skróty do pól akcesora .dt — xarray rozpoznaje je w groupby.
groupby_bins — kategorie po wartościach
Gdy chcemy zgrupować po przedziałach wartości (np. „lat 49-51, 51-53, 53-55”):
Wynik ma nową współrzędną lat_bins z obiektami pandas.Interval jako etykietami — czytelne, ale do zapisu do NetCDF trzeba je konwertować (np. na środki przedziałów).
rolling — okno przesuwne
import xarray as xrimport numpy as nparr = xr.DataArray(np.arange(10.0, dtype=float), dims="x")# Okno 3-elementowe, srednia:roll = arr.rolling(x=3).mean()print(roll.values)# [nan nan 1. 2. 3. 4. 5. 6. 7. 8.]# pierwsze dwie pozycje to NaN — brak pelnego okna# Z center=True — okno wycentrowane:roll_c = arr.rolling(x=3, center=True).mean()print(roll_c.values)# [nan 1. 2. 3. 4. 5. 6. 7. 8. nan]# Z min_periods=1 — wynik nawet dla niepelnego okna:roll_m = arr.rolling(x=3, min_periods=1).mean()print(roll_m.values)# [0. 0.5 1. 2. 3. 4. 5. 6. 7. 8. ]
rolling udostępnia operacje: mean, sum, min, max, std, var, count, median. Można też użyć .reduce(np.func) dla własnej funkcji.
Typowe zastosowanie hydrologiczne: 7-dniowa średnia ruchoma stanu wody (do wygładzania szumu) albo 30-dniowa suma opadu (do oceny zwilżenia zlewni).
Sumy miesięczne (mm):[ 16.4 24.9 67.4 95.7 125.7 150.8 157.5 145.7 83.2 57.1 28.4 24.7]Najmokrejszy miesiac: 2025-07-01T00:00:00.000000 z 157.5 mm
Praca z dużymi danymi — dask
xarray od początku był projektowany pod leniwą ewaluację (lazy evaluation) dla danych nie mieszczących się w pamięci. Zaplecze, którego używa się do tego, to dask — biblioteka do równoległych i leniwych obliczeń na tablicach.
Idea chunków
Tablica (time=10000, lat=720, lon=1440) w float32 ma rozmiar 41 GB. Nie zmieści się w RAM przeciętnego komputera. dask dzieli ją na chunki (np. (time=100, lat=720, lon=1440) — po 414 MB) i operuje na nich kolejno, trzymając jednocześnie w pamięci tylko kilka chunków.
xarray integruje to przezroczyście — po wywołaniu .chunk(...) operacje stają się leniwe:
# Lazy chunking — wymaga zainstalowanego pakietu daskimport xarray as xrimport numpy as npds = xr.Dataset( {"v": (("time", "x"), np.random.rand(1000, 100))}, coords={"time": np.arange(1000), "x": np.arange(100)},)# bez dask:try: chunked = ds.chunk({"time": 100})except ImportError as e: print("brak dask:", e)# brak dask: chunk manager 'dask' is not available...
Środowisko bazy nie ma dask
W requirements.txt brakuje dask. Powyższy chunk rzuci ImportError. Aby uruchomić:
pip install dask
Po zainstalowaniu — typowy przebieg pracy
# Wymaga: pip install daskimport xarray as xr# Otwarcie pliku z chunkowaniem — dane czytane sa leniwieds = xr.open_dataset("dane.nc", chunks={"time": 100})print(ds["temperatura"].data)# <dask.array.Array> — nie ndarray, tylko obietnica obliczen# Operacje sa lazy:srednia_miesieczna = ds["temperatura"].resample(time="MS").mean()# tutaj NIC sie jeszcze nie policzylo# Wymusza obliczenie:wynik = srednia_miesieczna.compute()# teraz: ndarray z policzonymi wartosciami
Kiedy używać a kiedy nie
dask opłaca się, gdy:
dane nie mieszczą się w RAM (rzędu dziesiątek GB i więcej);
operacje są „kłopotliwie równoległe” (embarrassingly parallel) — niezależnie po chunkach (np. resample, większość agregacji);
mamy więcej niż jeden rdzeń CPU lub klaster.
dask nie opłaca się, gdy:
dane są małe (poniżej GB) — narzut na zarządzanie chunkami przeważa zysk z równoległości;
algorytm jest globalny (np. odwrócenie macierzy, FFT na całej tablicy) — wymaga zebrania wszystkiego w jeden chunk;
pracujemy w interaktywnym notatniku, gdzie kolejność operacji jest nieprzewidywalna.
Reguła kciuka: zacznij bez dask. Jeśli pamięć siada lub obliczenia są wolne — wtedy chunk(...) na osi czasu (najczęściej) i dask.distributed.Client() do zrównoleglenia.
Praktyka — syntetyczny grid meteorologiczny
Tu zbudujemy pełny przykład: pole dziennych opadów nad Polską dla roku 2025 (365 dni × 20 punktów szerokości × 20 punktów długości). Wszystkie kroki: synteza, zapis, odczyt, agregacja i wykres.
Ten krok wymaga pip install netCDF4 (lub h5netcdf lub scipy) — patrz wcześniej. Jeśli nie ma backendu, pomiń krok 2 i krok 3 (odczyt) — pole da jest w pamięci, dalsze kroki są niezależne.
# Wymaga zainstalowanego netCDF4 / h5netcdf / scipyda.to_dataset().to_netcdf("/tmp/precip_2025.nc")# Plik powinien miec ok. 1.2 MB (1.17 MB dla 365*20*20*8 bytes / czesc na metadane)
# Suma miesięczna w każdym punkcie siatki:miesieczne = da.resample(time="MS").sum()print(miesieczne.shape)# (12, 20, 20)print(miesieczne.time.values[:3])# ['2025-01-01T00:00:00.000000' '2025-02-01T00:00:00.000000' '2025-03-01T00:00:00.000000']# Średnia z 12 miesięcznych sum (czyli „średni miesiąc" po całej Polsce):print("Srednie miesieczne sumy (mm, srednia po obszarze):")print(miesieczne.mean(dim=("lat", "lon")).values.round(1))
Widać sezonowość: minimum w styczniu i grudniu (~19 mm), maksimum w czerwcu-lipcu (~148-153 mm). To efekt formuły 0.5 + 4.5*(1+sin(...))/2.
Krok 5 — średnia po szerokościach (groupby_bins)
# Trzy pasma szerokosci geograficznej: 49-51, 51-53, 53-55:po_szerokosci = da.groupby_bins("lat", bins=[49, 51, 53, 55]).mean()print(po_szerokosci.shape)# (365, 3, 20)print(po_szerokosci.coords["lat_bins"].values)# [(49, 51] (51, 53] (53, 55]]# Średnia roczna w każdym pasie (po czasie i długości):roczna_w_pasach = po_szerokosci.mean(dim=("time", "lon"))print(roczna_w_pasach.values.round(3))
Output:
[2.762 2.75 2.757]
W tym syntetycznym przykładzie pole opadu nie zależy od lat — średnie w pasach są niemal identyczne (różnice to fluktuacje losowe). W realnych danych zobaczylibyśmy gradient.
Krok 6 — okno kroczące 7-dniowe
# 7-dniowa suma krocząca opadów (wskaznik zwilzenia zlewni):opad_7d = da.rolling(time=7, center=True).sum()print(opad_7d.shape)# (365, 20, 20)# Maksymalna 7-dniowa suma na obszarze# (sprawdzamy, kiedy spadl najwiekszy "epizod" w pierwszym punkcie):punkt = opad_7d.isel(lat=0, lon=0)max_idx = int(punkt.argmax())print(f"Max 7-dniowa suma w punkcie (49°N, 14°E):")print(f" data: {punkt.time.values[max_idx]}")print(f" suma: {float(punkt.isel(time=max_idx)):.1f} mm")
Output (deterministyczne):
Max 7-dniowa suma w punkcie (49°N, 14°E): data: 2025-07-08T00:00:00.000000 suma: 84.3 mm
Krok 7 — wizualizacja
xarray ma wbudowane nakładki na matplotlib (DataArray.plot()), które rozpoznają typ danych i wybierają sensowny styl:
import matplotlibmatplotlib.use("Agg") # backend bez GUI — do skryptuimport matplotlib.pyplot as plt# Wykres szeregu czasowego dla jednego punktu (49.32°N, 14.55°E):fig, ax = plt.subplots(figsize=(10, 4))da.isel(lat=0, lon=0).plot(ax=ax)ax.set_title("Dzienny opad w punkcie (lat=49°N, lon=14°E)")ax.set_ylabel("Opad [mm/dzień]")fig.tight_layout()fig.savefig("/tmp/opad_punkt.png", dpi=120)plt.close(fig)print("Zapisano /tmp/opad_punkt.png")
używa współrzędnych lat/lon jako osi (nie indeksów 0..19);
dodaje pasek koloru z opisem (z attrs["units"]);
podpisuje osie nazwami wymiarów.
Wykres miesięcznych sum (uśrednionych po obszarze):
import matplotlib.pyplot as pltm = da.resample(time="MS").sum().mean(dim=("lat", "lon"))# xarray.DataArray.plot ma .line/.hist/.imshow/.pcolormesh/.scatter/.step/.contour# — nie ma natywnego .bar(). Najprosciej narysowac bar przez matplotlib:fig, ax = plt.subplots(figsize=(10, 4))miesiace = [str(t)[:7] for t in m.time.values] # "YYYY-MM"ax.bar(miesiace, m.values)ax.set_title("Srednie miesieczne sumy opadu (synteza)")ax.set_ylabel("Suma [mm]")ax.tick_params(axis="x", rotation=45)fig.tight_layout()fig.savefig("/tmp/sumy_miesieczne.png", dpi=120)plt.close(fig)print("Zapisano /tmp/sumy_miesieczne.png")
Krok 8 — zapis pojedynczego rastra (średnia roczna) do GeoTIFF
Powstały plik można otworzyć w QGIS i wyświetlić jako warstwę rastrową.
Krok 9 — w jednym łańcuchu (styl funkcyjny)
xarray pozwala na płynne łączenie operacji w łańcuch:
# Najmokrejsze 5 dni w roku w polnocnej Polsce (lat >= 53):top5 = ( da .sel(lat=slice(53, 55)) .mean(dim=("lat", "lon")) .sortby(lambda x: -x) .isel(time=slice(0, 5)))print(top5.time.values)print(top5.values.round(2))
xarray był projektowany jako rozszerzenie wzorca pandas na wiele wymiarów — wiele elementów API jest celowo symetrycznych:
pandas
xarray
DataFrame
Dataset
Series
DataArray
df.loc[idx]
da.sel(dim=idx)
df.iloc[i]
da.isel(dim=i)
df.groupby("col").mean()
da.groupby("col").mean()
df.resample("MS").sum()
da.resample(time="MS").sum()
df.rolling(7).mean()
da.rolling(time=7).mean()
df.fillna(0)
da.fillna(0)
df.dropna()
da.dropna(dim=...)
df.melt() / df.stack()
da.stack(...) / da.to_dataframe()
df.to_csv(...)
da.to_pandas().to_csv(...) (przez pd)
pd.merge(df1, df2)
xr.merge([ds1, ds2])
pd.concat([df1, df2])
xr.concat([da1, da2], dim=...)
Główne różnice:
pandas zakłada 2D z heterogenicznymi typami kolumn; xarray zakłada N-D z jednorodnym typem w obiekcie DataArray (a wiele typów obok siebie organizuje w obiekcie Dataset).
pandas indeksuje wszystkim (wiersz × kolumna); xarray indeksuje po nazwach wymiarów — bezpieczniej przy ≥3D.
pandas ma natywne CSV, JSON, SQL; xarray ma natywne NetCDF, Zarr i integrację z dask.
Dla bardziej szczegółowego materiału o pandas: pandas.
Pułapki i dobre praktyki
Kolejność wymiarów ma znaczenie
import xarray as xrimport numpy as npa = xr.DataArray(np.zeros((3, 5)), dims=("x", "y"))print(a.shape)# (3, 5)b = a.transpose("y", "x")print(b.shape)# (5, 3)
.transpose(...) zmienia kolejność wymiarów (jak np.transpose, ale po nazwie). Większość operacji w xarray jest „świadoma wymiarów” — suma po time przy układzie (lat, time, lon) daje to samo co przy (time, lat, lon). Ale przy zapisie do NetCDF, GeoTIFF lub przy wymianie z NumPy kolejność wymiarów musi być znana i zgodna z konwencją.
Konwencja CF/GeoTIFF: (time, [z], lat, lon) lub (time, [z], y, x). Trzymanie się tego ułatwia interoperacyjność.
Coords po operacjach
Po agregacji wzdłuż osi ta oś znika — ale jej współrzędna może zostać jako skalarna:
import xarray as xrimport pandas as pdimport numpy as nptimes = pd.date_range("2025-01-01", periods=5)da = xr.DataArray(np.arange(15).reshape(5, 3), dims=("time", "site"), coords={"time": times, "site": ["A", "B", "C"]})# Mean po site — wymiar site znikam = da.mean(dim="site")print(m.dims)# ('time',)print("site" in m.coords)# False
W zależności od operacji niektóre coords są zachowywane (skalarne), inne usuwane. W razie wątpliwości sprawdź przez print(result).
attrs nie zawsze propagują — wymuszaj
Patrz wcześniej. Jeśli jednostki, opis i źródło muszą zostać po obliczeniu — opakowuj w with xr.set_options(keep_attrs=True): albo ustaw globalnie na początku skryptu:
import xarray as xrxr.set_options(keep_attrs=True)
Albo wprost po operacji:
import xarray as xrimport numpy as npa = xr.DataArray([1.0, 2.0], dims="x", attrs={"units": "m"})b = a * 0.001b.attrs["units"] = "km"b.attrs["long_name"] = "Distance"
Broadcasting po coords — przy każdej arytmetyce
Powtórzmy z naciskiem: gdy dwa obiekty DataArray mają wspólny wymiar, ale różne współrzędne, xarray automatycznie znajduje przecięcie. Może to dać niespodziewanie pustą tablicę:
import xarray as xra = xr.DataArray([1, 2, 3], dims="x", coords={"x": [1, 2, 3]})b = xr.DataArray([10, 20, 30], dims="x", coords={"x": [10, 20, 30]})print((a + b).values)# [] — bo zaden x sie nie pokrywaprint((a + b).x.values)# []
Jeśli chcemy zsumować po pozycji (a nie po etykiecie), trzeba albo zignorować coords (a.values + b.values), albo wyrównać etykiety przed operacją.
NaN w typach całkowitych — niemożliwe
import xarray as xrimport numpy as npa = xr.DataArray([1, 2, 3], dims="x") # intprint(a.dtype)# int64# Probujemy ustawic NaN — wymusza konwersje na floatb = a.where(a > 1)print(b.dtype)# float64print(b.values)# [nan 2. 3.]
Tak samo jak w NumPy: NaN istnieje tylko dla typów zmiennoprzecinkowych. xarray konwertuje automatycznie. Jeśli zależy nam na typie całkowitym z „brakiem danych” — używamy konkretnej wartości NoData (np. -9999) i własnej maski.
Wydajność: wektoryzuj, unikaj pętli
Tak samo jak w NumPy i pandas, w xarray pętla po wymiarach jest bardzo wolna. Wszystko, co jest agregacją lub operacją po wszystkich elementach, zrób bezpośrednio na obiekcie DataArray:
import xarray as xrimport numpy as npimport pandas as pd# Przykladowe pole (czas x lat x lon)da = xr.DataArray( np.random.default_rng(0).normal(0, 1, size=(365, 20, 20)), dims=("time", "lat", "lon"), coords={ "time": pd.date_range("2025-01-01", periods=365), "lat": np.linspace(49, 55, 20), "lon": np.linspace(14, 24.5, 20), },)# ŹLE — petla po dniach:# wyniki = []# for t in range(365):# wyniki.append(da.isel(time=t).mean().values)# DOBRZE — wektorowo:wyniki = da.mean(dim=("lat", "lon")).values # 365 wartosci, jedno wywolanieprint(wyniki.shape)# (365,)
Pierwsza wersja może być zauważalnie wolniejsza dla siatki 20×20 i znacząco wolniejsza dla siatki 1000×1000.
Kiedy xarray, a kiedy pandas
Reguła praktyczna:
Tabela 2D z heterogenicznymi typami kolumn (np. CSV ze stacjami, datami, statusami) → pandas.
Pole numeryczne 3D i więcej z jednorodnym typem (czas × lat × lon lub czas × stacja × parametr) → xarray.
Pojedynczy raster 2D bez czasu → numpy + rasterio (xarray to przerost).
Hybryda (tabela ze stacjami + szereg pomiarowy) → pandas na zewnątrz, ewentualnie konwersja na xarray przy operacjach typu „interpoluj po czasie”.
Wersje pakietów — kompatybilność wstecz
API xarray jest stosunkowo stabilne, ale część szczegółów (m.in. domyślne wartości typu keep_attrs, niektóre nazwy parametrów, wycofywane elementy) zmienia się między wersjami. Przy migracji z xarray 2023.x lub 2024.x na 2026.x warto:
przeczytać whatsnew w oficjalnej dokumentacji;
sprawdzić FutureWarning w wynikach pytest;
jeśli kod ma być długoterminowy — zafiksować wersję (xarray==2026.4.0 w requirements.txt).
Pliki NetCDF — konwencja CF
Jeśli plik trafi do narzędzia innego niż Twoje (QGIS, CDO, Panoply, modele klimatyczne), warto trzymać się konwencji CF:
units — jednostki (np. mm, kg m-2 s-1, degC, K).
long_name — opis czytelny dla człowieka.
standard_name — jedna z zatwierdzonych nazw CF (np. precipitation_flux, air_temperature).
_FillValue — wartość oznaczająca „brak danych” (czytana automatycznie przez xarray/QGIS).
Conventions — "CF-1.8" lub nowsza wersja, w atrybutach całego datasetu.
To inwestycja, która zwraca się przy każdej wymianie pliku z inną aplikacją.