1import xarray as xr
2import numpy as np
3import matplotlib.pyplot as plt
4import matplotlib.colors as mcolors
5import matplotlib.ticker as mticker
6import cartopy.crs as ccrs
7import cartopy.feature as cfeature
8
9# ── Store location ──────────────────────────────────────────────
10ZARR_PATH = 's3://us-west-2.opendata.source.coop/earthblox/cciwr/data/cciwr.zarr'
11
12# ── Open all three groups (anonymous S3 access, lazy — no data loaded yet) ──
13storage_options = {'anon': True}
14ds_rc    = xr.open_zarr(ZARR_PATH, group='RC',    consolidated=True, storage_options=storage_options)
15ds_futur = xr.open_zarr(ZARR_PATH, group='FUTUR', consolidated=True, storage_options=storage_options)
16ds_histo = xr.open_zarr(ZARR_PATH, group='HISTO', consolidated=True, storage_options=storage_options)
17
18ds_rc

1import xarray as xr
2import numpy as np
3import matplotlib.pyplot as plt
4import matplotlib.colors as mcolors
5import matplotlib.ticker as mticker
6import cartopy.crs as ccrs
7import cartopy.feature as cfeature
8
9# ── Store location ──────────────────────────────────────────────
10ZARR_PATH = 's3://us-west-2.opendata.source.coop/earthblox/cciwr/data/cciwr.zarr'
11
12# ── Open all three groups (anonymous S3 access, lazy — no data loaded yet) ──
13storage_options = {'anon': True}
14ds_rc    = xr.open_zarr(ZARR_PATH, group='RC',    consolidated=True, storage_options=storage_options)
15ds_futur = xr.open_zarr(ZARR_PATH, group='FUTUR', consolidated=True, storage_options=storage_options)
16ds_histo = xr.open_zarr(ZARR_PATH, group='HISTO', consolidated=True, storage_options=storage_options)
17
18ds_rc

1# ── User parameters ─────────────────────────────────────────────────────────
2POI_LAT  = 45.75   # Lyon city centre (snapped to nearest 0.5° cell)
3POI_LON  = 4.75
4SSP      = 'ssp585'
5HORIZON  = 2085
6SEASON   = 'yearly'
7
8# ── Select the nearest grid cell for the POI ────────────────────────────────
9# lat is descending, so method='nearest' works correctly with sel()
10poi_rc = ds_rc['qr'].sel(
11    ssp=SSP,
12    horizon=HORIZON,
13    season=SEASON,
14    lat=POI_LAT,
15    lon=POI_LON,
16)
17# poi_rc shape: (ghm, gcm)  — all ensemble members at this cell
18
19# Stack ghm × gcm into a flat ensemble dimension for quantile computation
20ensemble = poi_rc.stack(member=('ghm', 'gcm'))
21
22q10  = float(ensemble.quantile(0.10))
23q25  = float(ensemble.quantile(0.25))
24med  = float(ensemble.quantile(0.50))
25q75  = float(ensemble.quantile(0.75))
26q90  = float(ensemble.quantile(0.90))
27
28snapped_lat, snapped_lon = float(poi_rc.lat), float(poi_rc.lon)
29
30print(f'Scenario     : {SSP}, horizon {HORIZON}, season {SEASON}')
31print(f'Variable     : qr (groundwater recharge, relative change %)')
32print(f'  25th percentile : {100 * q25:+.1f} %')
33print(f'  Median          : {100 * med:+.1f} %')
34print(f'  75th percentile : {100 * q75:+.1f} %')

1# ── User parameters ─────────────────────────────────────────────────────────
2POI_LAT  = 45.75   # Lyon city centre (snapped to nearest 0.5° cell)
3POI_LON  = 4.75
4SSP      = 'ssp585'
5HORIZON  = 2085
6SEASON   = 'yearly'
7
8# ── Select the nearest grid cell for the POI ────────────────────────────────
9# lat is descending, so method='nearest' works correctly with sel()
10poi_rc = ds_rc['qr'].sel(
11    ssp=SSP,
12    horizon=HORIZON,
13    season=SEASON,
14    lat=POI_LAT,
15    lon=POI_LON,
16)
17# poi_rc shape: (ghm, gcm)  — all ensemble members at this cell
18
19# Stack ghm × gcm into a flat ensemble dimension for quantile computation
20ensemble = poi_rc.stack(member=('ghm', 'gcm'))
21
22q10  = float(ensemble.quantile(0.10))
23q25  = float(ensemble.quantile(0.25))
24med  = float(ensemble.quantile(0.50))
25q75  = float(ensemble.quantile(0.75))
26q90  = float(ensemble.quantile(0.90))
27
28snapped_lat, snapped_lon = float(poi_rc.lat), float(poi_rc.lon)
29
30print(f'Scenario     : {SSP}, horizon {HORIZON}, season {SEASON}')
31print(f'Variable     : qr (groundwater recharge, relative change %)')
32print(f'  25th percentile : {100 * q25:+.1f} %')
33print(f'  Median          : {100 * med:+.1f} %')
34print(f'  75th percentile : {100 * q75:+.1f} %')

1# ── Percentile-box chart across all three horizons ───────────────────────────
2# Reuses POI_LAT, POI_LON, SSP, SEASON from Case study 1 parameters.
3
4VARIABLE = 'qr'          # variable to chart
5HORIZONS       = [2030, 2055, 2085]
6HORIZON_LABELS = ['2016–2045', '2041–2070', '2071–2100']
7
8VAR_LABELS = {
9    'qtot':       'Total Runoff',
10    'qr':         'Groundwater Recharge',
11    'evap_total': 'Actual Evapotranspiration',
12}
13SSP_LABELS = {
14    'ssp126': 'SSP1 (RCP 2.6)',
15    'ssp370': 'SSP3 (RCP 7.0)',
16    'ssp585': 'SSP5 (RCP 8.5)',
17}
18
19# ── Collect percentiles for every horizon ────────────────────────────────────
20horizon_stats = {}
21snapped_lat_chart = snapped_lon_chart = None
22
23for h in HORIZONS:
24    # Exact match on string/categorical dims first; nearest-neighbour only for lat/lon
25    poi_h = (
26        ds_rc[VARIABLE]
27        .sel(ssp=SSP, horizon=h, season=SEASON)
28        .sel(lat=POI_LAT, lon=POI_LON, method='nearest')
29    )
30    ensemble_h = 100.0 * poi_h.stack(member=('ghm', 'gcm')).values
31
32    if snapped_lat_chart is None:
33        snapped_lat_chart = float(poi_h.lat)
34        snapped_lon_chart = float(poi_h.lon)
35
36    horizon_stats[h] = {
37        'q10': float(np.nanpercentile(ensemble_h, 10)),
38        'q25': float(np.nanpercentile(ensemble_h, 25)),
39        'med': float(np.nanpercentile(ensemble_h, 50)),
40        'q75': float(np.nanpercentile(ensemble_h, 75)),
41        'q90': float(np.nanpercentile(ensemble_h, 90)),
42        'n':   int(np.sum(~np.isnan(ensemble_h))),
43    }
44
45# ── Chart ────────────────────────────────────────────────────────────────────
46fig, ax = plt.subplots(figsize=(8, 4))
47
48x     = np.arange(len(HORIZONS))
49width = 0.55
50
51COLOR_OUTER  = '#c8c8c8'   # 10th – 90th percentile band
52COLOR_IQR    = '#787878'   # 25th – 75th IQR
53COLOR_MEDIAN = '#303030'   # median line
54
55for i, h in enumerate(HORIZONS):
56    s = horizon_stats[h]
57
58    # 10th – 90th band (lightest)
59    ax.bar(i, s['q90'] - s['q10'], bottom=s['q10'],
60           color=COLOR_OUTER, width=width, zorder=2, linewidth=0)
61
62    # 25th – 75th IQR (medium)
63    ax.bar(i, s['q75'] - s['q25'], bottom=s['q25'],
64           color=COLOR_IQR, width=width, zorder=3, linewidth=0)
65
66    # Median line
67    ax.hlines(s['med'], i - width / 2, i + width / 2,
68              colors=COLOR_MEDIAN, linewidths=2.5, zorder=4)
69
70# Reference line at 0 %
71ax.axhline(0, color='black', linewidth=0.9, zorder=1)
72
73# Axes formatting
74ax.set_xticks(x)
75ax.set_xticklabels(HORIZON_LABELS, fontsize=11)
76ax.yaxis.set_major_formatter(plt.FuncFormatter(lambda v, _: f'{v:.0f}%'))
77ax.set_ylabel(
78    f'Relative change in\n{VAR_LABELS.get(VARIABLE, VARIABLE)} since 1985–2014 [%]',
79    fontsize=10,
80)
81ax.tick_params(axis='y', labelsize=10)
82ax.set_xlim(-0.55, len(HORIZONS) - 0.45)
83ax.set_ylim(-100, 100)
84
85# Major + minor horizontal grid
86ax.yaxis.set_minor_locator(plt.MultipleLocator(5))
87ax.yaxis.grid(True, which='major', linewidth=0.5, color='lightgrey', zorder=0)
88ax.yaxis.grid(True, which='minor', linewidth=0.25, color='#e8e8e8', zorder=0)
89ax.set_axisbelow(True)
90ax.spines[['top', 'right']].set_visible(False)
91
92# Title
93ssp_label = SSP_LABELS.get(SSP, SSP)
94n_members = horizon_stats[HORIZONS[0]]['n']
95ax.set_title(
96    f'Lat: {snapped_lat_chart}°, Lon: {snapped_lon_chart}° – '
97    f'Seasonality: {SEASON.capitalize()} – {ssp_label}',
98    fontsize=11, pad=10,
99)
100
101# Footer — ensemble size + legend explanation
102ax.text(
103    0.5, -0.10,
104    f'Based on {n_members} ensemble members',
105    ha='center', va='top', transform=ax.transAxes,
106    fontstyle='italic', fontsize=9,
107)
108ax.text(
109    0.5, -0.18,
110    'Outer box: 10th–90th percentile  ·  Inner box: 25th–75th percentile  ·  Line: median',
111    ha='center', va='top', transform=ax.transAxes,
112    fontstyle='italic', fontsize=8, color='#555555',
113)
114
115plt.tight_layout()
116plt.show()

1# ── Percentile-box chart across all three horizons ───────────────────────────
2# Reuses POI_LAT, POI_LON, SSP, SEASON from Case study 1 parameters.
3
4VARIABLE = 'qr'          # variable to chart
5HORIZONS       = [2030, 2055, 2085]
6HORIZON_LABELS = ['2016–2045', '2041–2070', '2071–2100']
7
8VAR_LABELS = {
9    'qtot':       'Total Runoff',
10    'qr':         'Groundwater Recharge',
11    'evap_total': 'Actual Evapotranspiration',
12}
13SSP_LABELS = {
14    'ssp126': 'SSP1 (RCP 2.6)',
15    'ssp370': 'SSP3 (RCP 7.0)',
16    'ssp585': 'SSP5 (RCP 8.5)',
17}
18
19# ── Collect percentiles for every horizon ────────────────────────────────────
20horizon_stats = {}
21snapped_lat_chart = snapped_lon_chart = None
22
23for h in HORIZONS:
24    # Exact match on string/categorical dims first; nearest-neighbour only for lat/lon
25    poi_h = (
26        ds_rc[VARIABLE]
27        .sel(ssp=SSP, horizon=h, season=SEASON)
28        .sel(lat=POI_LAT, lon=POI_LON, method='nearest')
29    )
30    ensemble_h = 100.0 * poi_h.stack(member=('ghm', 'gcm')).values
31
32    if snapped_lat_chart is None:
33        snapped_lat_chart = float(poi_h.lat)
34        snapped_lon_chart = float(poi_h.lon)
35
36    horizon_stats[h] = {
37        'q10': float(np.nanpercentile(ensemble_h, 10)),
38        'q25': float(np.nanpercentile(ensemble_h, 25)),
39        'med': float(np.nanpercentile(ensemble_h, 50)),
40        'q75': float(np.nanpercentile(ensemble_h, 75)),
41        'q90': float(np.nanpercentile(ensemble_h, 90)),
42        'n':   int(np.sum(~np.isnan(ensemble_h))),
43    }
44
45# ── Chart ────────────────────────────────────────────────────────────────────
46fig, ax = plt.subplots(figsize=(8, 4))
47
48x     = np.arange(len(HORIZONS))
49width = 0.55
50
51COLOR_OUTER  = '#c8c8c8'   # 10th – 90th percentile band
52COLOR_IQR    = '#787878'   # 25th – 75th IQR
53COLOR_MEDIAN = '#303030'   # median line
54
55for i, h in enumerate(HORIZONS):
56    s = horizon_stats[h]
57
58    # 10th – 90th band (lightest)
59    ax.bar(i, s['q90'] - s['q10'], bottom=s['q10'],
60           color=COLOR_OUTER, width=width, zorder=2, linewidth=0)
61
62    # 25th – 75th IQR (medium)
63    ax.bar(i, s['q75'] - s['q25'], bottom=s['q25'],
64           color=COLOR_IQR, width=width, zorder=3, linewidth=0)
65
66    # Median line
67    ax.hlines(s['med'], i - width / 2, i + width / 2,
68              colors=COLOR_MEDIAN, linewidths=2.5, zorder=4)
69
70# Reference line at 0 %
71ax.axhline(0, color='black', linewidth=0.9, zorder=1)
72
73# Axes formatting
74ax.set_xticks(x)
75ax.set_xticklabels(HORIZON_LABELS, fontsize=11)
76ax.yaxis.set_major_formatter(plt.FuncFormatter(lambda v, _: f'{v:.0f}%'))
77ax.set_ylabel(
78    f'Relative change in\n{VAR_LABELS.get(VARIABLE, VARIABLE)} since 1985–2014 [%]',
79    fontsize=10,
80)
81ax.tick_params(axis='y', labelsize=10)
82ax.set_xlim(-0.55, len(HORIZONS) - 0.45)
83ax.set_ylim(-100, 100)
84
85# Major + minor horizontal grid
86ax.yaxis.set_minor_locator(plt.MultipleLocator(5))
87ax.yaxis.grid(True, which='major', linewidth=0.5, color='lightgrey', zorder=0)
88ax.yaxis.grid(True, which='minor', linewidth=0.25, color='#e8e8e8', zorder=0)
89ax.set_axisbelow(True)
90ax.spines[['top', 'right']].set_visible(False)
91
92# Title
93ssp_label = SSP_LABELS.get(SSP, SSP)
94n_members = horizon_stats[HORIZONS[0]]['n']
95ax.set_title(
96    f'Lat: {snapped_lat_chart}°, Lon: {snapped_lon_chart}° – '
97    f'Seasonality: {SEASON.capitalize()} – {ssp_label}',
98    fontsize=11, pad=10,
99)
100
101# Footer — ensemble size + legend explanation
102ax.text(
103    0.5, -0.10,
104    f'Based on {n_members} ensemble members',
105    ha='center', va='top', transform=ax.transAxes,
106    fontstyle='italic', fontsize=9,
107)
108ax.text(
109    0.5, -0.18,
110    'Outer box: 10th–90th percentile  ·  Inner box: 25th–75th percentile  ·  Line: median',
111    ha='center', va='top', transform=ax.transAxes,
112    fontstyle='italic', fontsize=8, color='#555555',
113)
114
115plt.tight_layout()
116plt.show()

1# ── User parameters ─────────────────────────────────────────────────────────
2SSP     = 'ssp585'
3HORIZON = 2085
4SEASON  = 'yearly'
5
6# Europe bounding box / Change by your area of interest.
7LAT_MIN, LAT_MAX = 34.0, 72.0
8LON_MIN, LON_MAX = -25.0, 45.0
9
10# ── Spatial slice — lat is descending so min/max order is reversed ───────────
11europe_rc = ds_rc['qtot'].sel(
12    ssp=SSP,
13    horizon=HORIZON,
14    season=SEASON,
15    lat=slice(LAT_MAX, LAT_MIN),   # descending axis
16    lon=slice(LON_MIN, LON_MAX),
17)
18# europe_rc shape: (ghm, gcm, lat, lon)
19
20# Ensemble median across all GHM × GCM pairs
21median_rc = europe_rc.median(dim=['ghm', 'gcm']).compute()
22# median_rc shape: (lat, lon)
23
24# ── Diverging colormap centred on 0 ─────────────────────────────────────────
25vmax = float(np.nanpercentile(np.abs(median_rc.values), 95))
26norm = mcolors.TwoSlopeNorm(vmin=-vmax, vcenter=0, vmax=vmax)
27
28# ── Map ──────────────────────────────────────────────────────────────────────
29fig, ax = plt.subplots(
30    figsize=(8, 7),
31    subplot_kw={'projection': ccrs.PlateCarree()},
32)
33ax.set_extent([LON_MIN, LON_MAX, LAT_MIN, LAT_MAX], crs=ccrs.PlateCarree())
34
35img = ax.pcolormesh(
36    median_rc.lon.values,
37    median_rc.lat.values,
38    median_rc.values,
39    norm=norm,
40    cmap='RdBu',
41    transform=ccrs.PlateCarree(),
42)
43
44# Overlay borders and coastlines
45ax.add_feature(cfeature.BORDERS, linewidth=0.7, edgecolor='black')
46ax.add_feature(cfeature.COASTLINE, linewidth=0.8)
47ax.add_feature(cfeature.RIVERS, linewidth=0.4, edgecolor='steelblue', alpha=0.6)
48ax.gridlines(draw_labels=False, linewidth=0.3, color='grey', alpha=0.5)
49
50cbar = plt.colorbar(img, ax=ax, orientation='vertical', pad=0.02, fraction=0.03)
51cbar.set_label('Relative change in total runoff', fontsize=10)
52cbar.formatter = mticker.PercentFormatter(xmax=1.0, decimals=0)
53cbar.update_ticks()
54
55ax.set_title(
56    f'Ensemble-median RC — qtot | {SSP}, horizon {HORIZON}, {SEASON}\n'
57    f'(n = {len(ds_rc.ghm) * len(ds_rc.gcm)} ensemble members)',
58    fontsize=11,
59)
60
61plt.tight_layout()
62plt.show()

1# ── User parameters ─────────────────────────────────────────────────────────
2SSP     = 'ssp585'
3HORIZON = 2085
4SEASON  = 'yearly'
5
6# Europe bounding box / Change by your area of interest.
7LAT_MIN, LAT_MAX = 34.0, 72.0
8LON_MIN, LON_MAX = -25.0, 45.0
9
10# ── Spatial slice — lat is descending so min/max order is reversed ───────────
11europe_rc = ds_rc['qtot'].sel(
12    ssp=SSP,
13    horizon=HORIZON,
14    season=SEASON,
15    lat=slice(LAT_MAX, LAT_MIN),   # descending axis
16    lon=slice(LON_MIN, LON_MAX),
17)
18# europe_rc shape: (ghm, gcm, lat, lon)
19
20# Ensemble median across all GHM × GCM pairs
21median_rc = europe_rc.median(dim=['ghm', 'gcm']).compute()
22# median_rc shape: (lat, lon)
23
24# ── Diverging colormap centred on 0 ─────────────────────────────────────────
25vmax = float(np.nanpercentile(np.abs(median_rc.values), 95))
26norm = mcolors.TwoSlopeNorm(vmin=-vmax, vcenter=0, vmax=vmax)
27
28# ── Map ──────────────────────────────────────────────────────────────────────
29fig, ax = plt.subplots(
30    figsize=(8, 7),
31    subplot_kw={'projection': ccrs.PlateCarree()},
32)
33ax.set_extent([LON_MIN, LON_MAX, LAT_MIN, LAT_MAX], crs=ccrs.PlateCarree())
34
35img = ax.pcolormesh(
36    median_rc.lon.values,
37    median_rc.lat.values,
38    median_rc.values,
39    norm=norm,
40    cmap='RdBu',
41    transform=ccrs.PlateCarree(),
42)
43
44# Overlay borders and coastlines
45ax.add_feature(cfeature.BORDERS, linewidth=0.7, edgecolor='black')
46ax.add_feature(cfeature.COASTLINE, linewidth=0.8)
47ax.add_feature(cfeature.RIVERS, linewidth=0.4, edgecolor='steelblue', alpha=0.6)
48ax.gridlines(draw_labels=False, linewidth=0.3, color='grey', alpha=0.5)
49
50cbar = plt.colorbar(img, ax=ax, orientation='vertical', pad=0.02, fraction=0.03)
51cbar.set_label('Relative change in total runoff', fontsize=10)
52cbar.formatter = mticker.PercentFormatter(xmax=1.0, decimals=0)
53cbar.update_ticks()
54
55ax.set_title(
56    f'Ensemble-median RC — qtot | {SSP}, horizon {HORIZON}, {SEASON}\n'
57    f'(n = {len(ds_rc.ghm) * len(ds_rc.gcm)} ensemble members)',
58    fontsize=11,
59)
60
61plt.tight_layout()
62plt.show()

1# ── User parameters ─────────────────────────────────────────────────────────
2SSP       = 'ssp585'
3HORIZON   = 2085
4SEASON    = 'yearly'    # could also do seasonal slices (e.g. 'DJF', 'MAM', 'JJA', 'SON') for a more seasonal perspective
5VARIABLE  = 'qr'        # groundwater recharge
6THRESHOLD = -0.20       # drying stress: count members with RC below -20%
7
8# Europe bounding box (same extent as Case study 2).
9LAT_MIN, LAT_MAX = 34.0, 72.0
10LON_MIN, LON_MAX = -25.0, 45.0
11
12# ── Spatial slice — lat is descending so min/max order is reversed ───────────
13europe_rc = ds_rc[VARIABLE].sel(
14    ssp=SSP,
15    horizon=HORIZON,
16    season=SEASON,
17    lat=slice(LAT_MAX, LAT_MIN),   # descending axis
18    lon=slice(LON_MIN, LON_MAX),
19)
20# europe_rc shape: (ghm, gcm, lat, lon)
21
22# ── Model consensus: fraction of members with RC below the drying threshold ──
23# Count members satisfying the stress, ignoring NaNs, then divide by valid count.
24# The number of valid members varies from pixel to pixel (it is NOT always 20),
25# so we normalise by the per-pixel valid count rather than a fixed denominator.
26below     = (europe_rc < THRESHOLD)
27n_below   = below.sum(dim=['ghm', 'gcm'])
28n_valid   = europe_rc.notnull().sum(dim=['ghm', 'gcm'])
29consensus = (n_below / n_valid).where(n_valid > 0).compute()   # 0–1, NaN where no data
30# consensus shape: (lat, lon)
31
32# ── Sequential colormap (0–100% agreement) ──────────────────────────────────
33norm = mcolors.Normalize(vmin=0.0, vmax=1.0)
34
35# ── Map ──────────────────────────────────────────────────────────────────────
36fig, ax = plt.subplots(
37    figsize=(8, 7),
38    subplot_kw={'projection': ccrs.PlateCarree()},
39)
40ax.set_extent([LON_MIN, LON_MAX, LAT_MIN, LAT_MAX], crs=ccrs.PlateCarree())
41
42img = ax.pcolormesh(
43    consensus.lon.values,
44    consensus.lat.values,
45    consensus.values,
46    norm=norm,
47    cmap='YlOrRd',
48    transform=ccrs.PlateCarree(),
49)
50
51# Overlay borders and coastlines
52ax.add_feature(cfeature.BORDERS, linewidth=0.7, edgecolor='black')
53ax.add_feature(cfeature.COASTLINE, linewidth=0.8)
54ax.add_feature(cfeature.RIVERS, linewidth=0.4, edgecolor='steelblue', alpha=0.6)
55ax.gridlines(draw_labels=False, linewidth=0.3, color='grey', alpha=0.5)
56
57cbar = plt.colorbar(img, ax=ax, orientation='vertical', pad=0.02, fraction=0.03)
58cbar.set_label(f'Model consensus on RC < {THRESHOLD:+.0%} (share of members)', fontsize=10)
59cbar.formatter = mticker.PercentFormatter(xmax=1.0, decimals=0)
60cbar.update_ticks()
61
62ax.set_title(
63    f'Model consensus on drying — {VAR_LABELS.get(VARIABLE, VARIABLE)} '
64    f'RC < {THRESHOLD:+.0%}\n'
65    f'{SSP}, horizon {HORIZON}, {SEASON}',
66    fontsize=11,
67)
68
69plt.tight_layout()
70plt.show()

1# ── User parameters ─────────────────────────────────────────────────────────
2SSP       = 'ssp585'
3HORIZON   = 2085
4SEASON    = 'yearly'    # could also do seasonal slices (e.g. 'DJF', 'MAM', 'JJA', 'SON') for a more seasonal perspective
5VARIABLE  = 'qr'        # groundwater recharge
6THRESHOLD = -0.20       # drying stress: count members with RC below -20%
7
8# Europe bounding box (same extent as Case study 2).
9LAT_MIN, LAT_MAX = 34.0, 72.0
10LON_MIN, LON_MAX = -25.0, 45.0
11
12# ── Spatial slice — lat is descending so min/max order is reversed ───────────
13europe_rc = ds_rc[VARIABLE].sel(
14    ssp=SSP,
15    horizon=HORIZON,
16    season=SEASON,
17    lat=slice(LAT_MAX, LAT_MIN),   # descending axis
18    lon=slice(LON_MIN, LON_MAX),
19)
20# europe_rc shape: (ghm, gcm, lat, lon)
21
22# ── Model consensus: fraction of members with RC below the drying threshold ──
23# Count members satisfying the stress, ignoring NaNs, then divide by valid count.
24# The number of valid members varies from pixel to pixel (it is NOT always 20),
25# so we normalise by the per-pixel valid count rather than a fixed denominator.
26below     = (europe_rc < THRESHOLD)
27n_below   = below.sum(dim=['ghm', 'gcm'])
28n_valid   = europe_rc.notnull().sum(dim=['ghm', 'gcm'])
29consensus = (n_below / n_valid).where(n_valid > 0).compute()   # 0–1, NaN where no data
30# consensus shape: (lat, lon)
31
32# ── Sequential colormap (0–100% agreement) ──────────────────────────────────
33norm = mcolors.Normalize(vmin=0.0, vmax=1.0)
34
35# ── Map ──────────────────────────────────────────────────────────────────────
36fig, ax = plt.subplots(
37    figsize=(8, 7),
38    subplot_kw={'projection': ccrs.PlateCarree()},
39)
40ax.set_extent([LON_MIN, LON_MAX, LAT_MIN, LAT_MAX], crs=ccrs.PlateCarree())
41
42img = ax.pcolormesh(
43    consensus.lon.values,
44    consensus.lat.values,
45    consensus.values,
46    norm=norm,
47    cmap='YlOrRd',
48    transform=ccrs.PlateCarree(),
49)
50
51# Overlay borders and coastlines
52ax.add_feature(cfeature.BORDERS, linewidth=0.7, edgecolor='black')
53ax.add_feature(cfeature.COASTLINE, linewidth=0.8)
54ax.add_feature(cfeature.RIVERS, linewidth=0.4, edgecolor='steelblue', alpha=0.6)
55ax.gridlines(draw_labels=False, linewidth=0.3, color='grey', alpha=0.5)
56
57cbar = plt.colorbar(img, ax=ax, orientation='vertical', pad=0.02, fraction=0.03)
58cbar.set_label(f'Model consensus on RC < {THRESHOLD:+.0%} (share of members)', fontsize=10)
59cbar.formatter = mticker.PercentFormatter(xmax=1.0, decimals=0)
60cbar.update_ticks()
61
62ax.set_title(
63    f'Model consensus on drying — {VAR_LABELS.get(VARIABLE, VARIABLE)} '
64    f'RC < {THRESHOLD:+.0%}\n'
65    f'{SSP}, horizon {HORIZON}, {SEASON}',
66    fontsize=11,
67)
68
69plt.tight_layout()
70plt.show()