Note
Go to the end to download the full example code.
Köppen-Geiger climate classification (Beck et al.)#
The Köppen-Geiger scheme sorts the world into 30 climate classes from monthly temperature and precipitation: tropical (A), arid (B), temperate (C), cold (D) and polar (E), each subdivided by seasonality and summer warmth. Beck et al. (2023) rebuilt the maps at 1 km from station-constrained climatologies for five historical 30-year periods, and projected them to 2041-2070 and 2071-2099 under seven CMIP6 scenarios.
One route, the paper’s figshare archive, needs no account: a single 1.3 GB zip
that is read in place, one GeoTIFF at a time. load takes period=,
scenario= and resolution= (1°, 0.5°, 0.1° or 1 km); search lists
the archive’s layout without downloading it.
The classes come back as uint8 values with CF flag_values,
flag_meanings and flag_colors attributes, so esd.plotting.categorical
draws the paper’s own legend. The first figure is the Pacific Northwest at the
default 0.1°; the second is the Mount Rainier box at 1 km, present day beside
the SSP5-8.5 end-of-century projection, with a table of what changed.
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import easysnowdata as esd
print(f"easysnowdata {esd.__version__}")
easysnowdata 0.3.9.dev1+gd97689c2f
Where the product comes from, and what the route asks for.
for src in esd.catalog.get("koppen-geiger").sources:
print(f"{src.id:22} {src.title:45} {', '.join(src.requires) or 'no account'}")
figshare figshare zip no account
What the archive holds, without downloading it: one row per period, scenario and resolution.
inventory_df = esd.climate.koppen_geiger.search()
print(f"{len(inventory_df)} rasters in the archive")
print(inventory_df[inventory_df["resolution"] == "1 km"].head(8).to_string(index=False))
72 rasters in the archive
period scenario resolution resolution_m member
1901_1930 NaN 1 km 1000.0 1901_1930/koppen_geiger_0p00833333.tif
1931_1960 NaN 1 km 1000.0 1931_1960/koppen_geiger_0p00833333.tif
1961_1990 NaN 1 km 1000.0 1961_1990/koppen_geiger_0p00833333.tif
1991_2020 NaN 1 km 1000.0 1991_2020/koppen_geiger_0p00833333.tif
2041_2070 ssp119 1 km 1000.0 2041_2070/ssp119/koppen_geiger_0p00833333.tif
2041_2070 ssp126 1 km 1000.0 2041_2070/ssp126/koppen_geiger_0p00833333.tif
2041_2070 ssp245 1 km 1000.0 2041_2070/ssp245/koppen_geiger_0p00833333.tif
2041_2070 ssp370 1 km 1000.0 2041_2070/ssp370/koppen_geiger_0p00833333.tif
The Pacific Northwest at the default 0.1° for the default period, 1991-2020. The legend lists only the classes that occur in the array; the value 0 is ocean and is kept as the nodata sentinel, so it draws as nothing.
pnw_da = esd.climate.koppen_geiger.load((-125.0, 42.0, -116.0, 49.5))
print(pnw_da)
ax = esd.plotting.categorical(pnw_da, title="Köppen-Geiger classes, 1991-2020, 0.1°")
ax.figure.tight_layout()
<xarray.DataArray 'koppen_geiger_class' (latitude: 75, longitude: 90)> Size: 7kB
dask.array<getitem, shape=(75, 90), dtype=uint8, chunksize=(75, 90), chunktype=numpy.ndarray>
Coordinates:
* latitude (latitude) float64 600B 49.45 49.35 49.25 ... 42.25 42.15 42.05
* longitude (longitude) float64 720B -124.9 -124.8 -124.7 ... -116.1 -116.0
spatial_ref int64 8B 0
Attributes: (12/22)
AREA_OR_POINT: Area
scale_factor: 1.0
add_offset: 0.0
_FillValue: 0
nodata: 0
long_name: Köppen-Geiger climate class
... ...
license: CC BY 4.0
easysnowdata_version: 0.3.9.dev1+gd97689c2f
doi: 10.1038/s41597-023-02549-6
period: 1991_2020
resolution: 0.1 degree
archive_member: 1991_2020/koppen_geiger_0p1.tif
Mount Rainier at 1 km: the present period beside the end of the century under
SSP5-8.5. all_classes puts every class in the legend, so the same legend
serves both panels, in two columns to keep it shorter than the map.
The 1 km grid is geographic; the box is loaded with a 2 km margin and both
periods are drawn on its UTM zone at 1 km, nearest so the classes survive.
box = esd.parse_aoi((-122.6, 46.4, -120.9, 47.3))
grid = box.to_geobox(crs="utm", resolution=1000)
present_da = esd.climate.koppen_geiger.load(box.buffer(2000), resolution="1 km")
future_da = esd.climate.koppen_geiger.load(
box.buffer(2000), period="2071_2099", scenario="ssp585", resolution="1 km"
)
print(present_da.attrs["archive_member"], "->", future_da.attrs["archive_member"])
fig, axes = plt.subplots(1, 2, figsize=(12.5, 5.2), sharey=True)
esd.plotting.categorical(
present_da.odc.reproject(grid, resampling="nearest"),
ax=axes[0],
legend=False,
title="1991-2020",
)
esd.plotting.categorical(
future_da.odc.reproject(grid, resampling="nearest"),
ax=axes[1],
title="2071-2099, SSP5-8.5",
legend_kwargs={"all_classes": True, "ncol": 2, "fontsize": 7},
)
fig.tight_layout()
1991_2020/koppen_geiger_0p00833333.tif -> 2071_2099/ssp585/koppen_geiger_0p00833333.tif
How much of the box changes class, and to what. esd.processing.flags
reads the CF flag table back out of the attrs. Land pixels per class: the
cold dry-summer belts (Dsc, Dsb) and the summit tundra (ET) shrink, and the
hot-summer Csa, absent today, takes over most of what was warm-summer Csb.
This is a 1 km grid in degrees, so the counts are pixels, not km².
symbols = esd.processing.flags(present_da).set_index("value")["meaning"]
land = (present_da.values != 0) & (future_da.values != 0)
table_df = (
pd.DataFrame(
{
"1991-2020": pd.Series(present_da.values[land].ravel()).value_counts(),
"2071-2099 SSP5-8.5": pd.Series(
future_da.values[land].ravel()
).value_counts(),
}
)
.fillna(0)
.astype(int)
)
table_df.index = [symbols[v] for v in table_df.index]
table_df["change"] = table_df["2071-2099 SSP5-8.5"] - table_df["1991-2020"]
print(table_df.sort_values("change").to_string())
changed = np.mean(present_da.values[land] != future_da.values[land])
print(f"{changed:.0%} of the land pixels change class")
1991-2020 2071-2099 SSP5-8.5 change
Dsc 4368 275 -4093
Dsb 8477 4760 -3717
Csb 10430 8642 -1788
Cfb 318 16 -302
ET 298 70 -228
Dfb 58 0 -58
EF 17 0 -17
Dfc 10 0 -10
BSk 8 0 -8
Dsa 7 0 -7
Csc 1 0 -1
Cfa 0 2 2
Csa 10 10237 10227
81% of the land pixels change class