Skip to content

Commit b96ba16

Browse files
bennyistantobennyistanto
committed
Update maps
1 parent 992be6a commit b96ba16

File tree

6 files changed

+282
-7
lines changed

6 files changed

+282
-7
lines changed

‎README.md‎

Lines changed: 6 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -27,9 +27,13 @@
2727

2828
## Global Output
2929

30-
![SPI-12 computed from CHIRPS v3 global (0.05°) — December 2025](./docs/images/global-spi12-202512.png)
30+
SPI-12 (Gamma) calculated from **CHIRPS v3** at 0.05° resolution.
3131

32-
SPI-12 (Gamma) calculated from CHIRPS v3 global at 0.05° resolution.
32+
![SPI-12 computed from global CHIRPS v3 dataset (0.05°) — December 2025](./docs/images/global-spi12-202512.png)
33+
34+
SPEI-12 (Pearson III) calculated from **TerraClimate** at 0.0417° ~ 4km resolution.
35+
36+
![SPEI-12 computed from global TerraClimate dataset (0.0417° ~ 4km) — December 2024](./docs/images/global-spei12-202412.png)
3337

3438
## Credits
3539

‎docs/images/global-spei12-202412.png‎

808 KB
Loading

‎docs/images/global-spi12-202412.png‎

855 KB
Loading

‎docs/index.qmd‎

Lines changed: 12 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -168,9 +168,19 @@ WMO 11-category drought/wet classification. CF-compliant NetCDF output. Run theo
168168

169169
## Global Output
170170

171-
![SPI-12 computed from CHIRPS v3 global (0.05°) — December 2025](images/global-spi12-202512.png)
171+
::: {.panel-tabset}
172+
### SPI 12-month from CHIRPS
173+
174+
![SPI-12 computed from global CHIRPS v3 dataset (0.05°) — December 2025](images/global-spi12-202512.png)
175+
176+
SPI-12 (Gamma) calculated from **CHIRPS v3** at 0.05° resolution.
172177

173-
SPI-12 (Gamma) calculated from **CHIRPS v3 global** at 0.05° resolution.
178+
### SPEI 12-month from TerraClimate
179+
180+
![SPEI-12 computed from global TerraClimate dataset (0.0417° ~ 4km) — December 2024](images/global-spei12-202412.png)
181+
182+
SPEI-12 (Pearson III) calculated from **TerraClimate** at 0.0417° ~ 4km resolution.
183+
:::
174184

175185
---
176186

‎tests/plot_global_spei.py‎

Lines changed: 261 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,261 @@
1+
"""
2+
Plot Global SPEI-12 Map
3+
4+
Standalone script to visualize global SPEI-12 output from TerraClimate data.
5+
Produces a publication-ready map with WMO 11-category classification,
6+
equal-area projection, and country boundaries.
7+
8+
Usage:
9+
python tests/plot_global_spei.py
10+
11+
Input:
12+
global/output/netcdf/wld_cli_terraclimate_spei_pearson3_12_month.nc
13+
14+
Output:
15+
docs/images/global-spei12-202412.png
16+
"""
17+
18+
import sys
19+
from pathlib import Path
20+
21+
# Add src to path
22+
sys.path.insert(0, str(Path(__file__).parent.parent / 'src'))
23+
24+
import numpy as np
25+
import xarray as xr
26+
import matplotlib.pyplot as plt
27+
import matplotlib.colors as mcolors
28+
import matplotlib.ticker as mticker
29+
import cartopy.crs as ccrs
30+
import cartopy.feature as cfeature
31+
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
32+
33+
34+
# ============================================================================
35+
# Configuration
36+
# ============================================================================
37+
38+
# Input / output paths
39+
BASE_DIR = Path(__file__).parent.parent
40+
INPUT_FILE = BASE_DIR / 'global' / 'output' / 'netcdf' / 'wld_cli_terraclimate_spei_pearson3_12_month.nc'
41+
OUTPUT_FILE = BASE_DIR / 'docs' / 'images' / 'global-spei12-202412.png'
42+
43+
# Time slice to plot
44+
TARGET_DATE = '2024-12-01'
45+
46+
# Variable name in the NetCDF
47+
VAR_NAME = 'spei_pearson3_12_month'
48+
49+
# Index label (used in titles and colorbar)
50+
INDEX_LABEL = 'Standardized Precipitation Evapotranspiration Index (Pearson III), 12-month'
51+
DATA_SOURCE = 'Data: TerraClimate'
52+
53+
# Map settings
54+
PROJECTION = ccrs.EqualEarth() # Equal-area projection
55+
DATA_CRS = ccrs.PlateCarree() # Data is in lat/lon
56+
FIGSIZE = (12, 8)
57+
DPI = 200
58+
59+
# ============================================================================
60+
# WMO 11-Category Classification
61+
# ============================================================================
62+
63+
# Boundaries for 11 categories
64+
# -inf, -2.0, -1.5, -1.2, -0.7, -0.5, 0.5, 0.7, 1.2, 1.5, 2.0, +inf
65+
SPI_BOUNDS = [-3.5, -2.0, -1.5, -1.2, -0.7, -0.5, 0.5, 0.7, 1.2, 1.5, 2.0, 3.5]
66+
67+
SPI_COLORS = [
68+
'#760005', # Exceptionally dry (< -2.0)
69+
'#ec0013', # Extremely dry (-2.0 to -1.5)
70+
'#ffa938', # Severely dry (-1.5 to -1.2)
71+
'#fdd28a', # Moderately dry (-1.2 to -0.7)
72+
'#fefe53', # Abnormally dry (-0.7 to -0.5)
73+
'#ffffff', # Near normal (-0.5 to +0.5)
74+
'#a2fd6e', # Abnormally moist (+0.5 to +0.7)
75+
'#00b44a', # Moderately moist (+0.7 to +1.2)
76+
'#008180', # Very moist (+1.2 to +1.5)
77+
'#2a23eb', # Extremely moist (+1.5 to +2.0)
78+
'#a21fec', # Exceptionally moist (> +2.0)
79+
]
80+
81+
SPI_LABELS = [
82+
'Exceptionally\nDry',
83+
'Extremely\nDry',
84+
'Severely\nDry',
85+
'Moderately\nDry',
86+
'Abnormally\nDry',
87+
'Near\nNormal',
88+
'Abnormally\nMoist',
89+
'Moderately\nMoist',
90+
'Very\nMoist',
91+
'Extremely\nMoist',
92+
'Exceptionally\nMoist',
93+
]
94+
95+
96+
def build_colormap():
97+
"""Build a discrete colormap and norm for the WMO 11-category classification."""
98+
cmap = mcolors.ListedColormap(SPI_COLORS)
99+
norm = mcolors.BoundaryNorm(SPI_BOUNDS, cmap.N)
100+
return cmap, norm
101+
102+
103+
def load_data(filepath: Path, target_date: str, var_name: str) -> xr.DataArray:
104+
"""Load a single time slice from the global SPEI NetCDF."""
105+
print(f"Opening: {filepath.name}")
106+
print(f" (file size: {filepath.stat().st_size / 1e9:.1f} GB)")
107+
108+
# Open dataset — only decode the time slice we need
109+
ds = xr.open_dataset(filepath)
110+
111+
# Select the target date
112+
da = ds[var_name].sel(time=target_date, method='nearest')
113+
actual_time = str(da.time.values)[:10]
114+
print(f" Selected time: {actual_time}")
115+
116+
# Load into memory (single time slice)
117+
print(f" Loading data slice ({da.shape[0]} x {da.shape[1]}) ...")
118+
da = da.load()
119+
ds.close()
120+
121+
valid = int(np.isfinite(da.values).sum())
122+
total = da.size
123+
print(f" Valid cells: {valid:,} / {total:,} ({100*valid/total:.1f}%)")
124+
125+
return da
126+
127+
128+
def plot_global_spei(da: xr.DataArray, output_path: Path):
129+
"""Create publication-ready global SPEI map."""
130+
print("\nCreating map ...")
131+
132+
cmap, norm = build_colormap()
133+
134+
# Create figure — map fills most of the space, extra room at bottom for legend labels
135+
fig = plt.figure(figsize=FIGSIZE, facecolor='white')
136+
ax = fig.add_axes([0.02, 0.17, 0.96, 0.74], projection=PROJECTION)
137+
138+
# Background
139+
ax.set_global()
140+
ax.set_facecolor('#d9d9d9') # Light gray for ocean/background
141+
142+
# Plot SPEI data
143+
im = ax.pcolormesh(
144+
da.lon.values, da.lat.values, da.values,
145+
transform=DATA_CRS,
146+
cmap=cmap,
147+
norm=norm,
148+
shading='auto',
149+
rasterized=True, # Keeps file size reasonable
150+
)
151+
152+
# Country boundaries (50m resolution for better detail)
153+
ax.add_feature(
154+
cfeature.NaturalEarthFeature(
155+
'cultural', 'admin_0_boundary_lines_land', '50m',
156+
edgecolor='#404040', facecolor='none',
157+
),
158+
linewidth=0.3,
159+
alpha=0.7,
160+
)
161+
162+
# Coastlines (50m resolution)
163+
ax.add_feature(
164+
cfeature.NaturalEarthFeature(
165+
'physical', 'coastline', '50m',
166+
edgecolor='#202020', facecolor='none',
167+
),
168+
linewidth=0.4,
169+
)
170+
171+
# Gridlines
172+
gl = ax.gridlines(
173+
draw_labels=False,
174+
linewidth=0.3,
175+
color='gray',
176+
alpha=0.4,
177+
linestyle='--',
178+
)
179+
gl.xlocator = mticker.FixedLocator(range(-180, 181, 30))
180+
gl.ylocator = mticker.FixedLocator(range(-60, 81, 30))
181+
182+
# Title and subtitle with clear spacing
183+
date_label = f"{da.time.dt.strftime('%B %Y').values}"
184+
fig.text(
185+
0.5, 0.97, INDEX_LABEL,
186+
ha='center', va='top',
187+
fontsize=14, fontweight='bold',
188+
)
189+
fig.text(
190+
0.5, 0.93, f'as of {date_label}',
191+
ha='center', va='top',
192+
fontsize=11, color='#404040',
193+
)
194+
195+
# Colorbar — horizontal below the map, ticks at class boundaries
196+
cbar_left = 0.12
197+
cbar_width = 0.76
198+
cbar_ax = fig.add_axes([cbar_left, 0.08, cbar_width, 0.018])
199+
cbar = fig.colorbar(
200+
im, cax=cbar_ax, orientation='horizontal',
201+
extend='neither',
202+
ticks=[-2.0, -1.5, -1.2, -0.7, -0.5, 0.5, 0.7, 1.2, 1.5, 2.0],
203+
)
204+
cbar.ax.set_xticklabels(
205+
['-2.0', '-1.5', '-1.2', '-0.7', '-0.5', '0.5', '0.7', '1.2', '1.5', '2.0'],
206+
fontsize=7,
207+
)
208+
cbar.ax.set_xlabel(
209+
INDEX_LABEL,
210+
fontsize=9, labelpad=6,
211+
)
212+
213+
# Category labels above each color segment — use data coordinates on cbar axis
214+
label_fontsize = 5.5
215+
216+
for i, label in enumerate(SPI_LABELS):
217+
# Midpoint in data space — the colorbar x-axis IS in data space
218+
mid = (SPI_BOUNDS[i] + SPI_BOUNDS[i + 1]) / 2.0
219+
cbar.ax.text(mid, 1.3, label, transform=cbar.ax.get_xaxis_transform(),
220+
ha='center', va='bottom',
221+
fontsize=label_fontsize, color='#404040')
222+
223+
# Attribution text
224+
fig.text(0.02, 0.02, DATA_SOURCE, fontsize=8, color='#606060',
225+
ha='left', va='bottom')
226+
fig.text(0.98, 0.02, 'GOST/DEC Data Group/WBG', fontsize=8, color='#606060',
227+
ha='right', va='bottom')
228+
229+
# Save — use pad_inches to control whitespace
230+
plt.savefig(output_path, dpi=DPI, bbox_inches='tight', pad_inches=0.1,
231+
facecolor='white', edgecolor='none')
232+
plt.close()
233+
print(f"Saved: {output_path}")
234+
print(f" File size: {output_path.stat().st_size / 1e6:.1f} MB")
235+
236+
237+
def main():
238+
"""Main entry point."""
239+
print("=" * 60)
240+
print(" GLOBAL SPEI-12 MAP — DECEMBER 2024")
241+
print("=" * 60)
242+
243+
# Check input exists
244+
if not INPUT_FILE.exists():
245+
print(f"ERROR: Input file not found: {INPUT_FILE}")
246+
sys.exit(1)
247+
248+
# Ensure output directory exists
249+
OUTPUT_FILE.parent.mkdir(parents=True, exist_ok=True)
250+
251+
# Load data
252+
da = load_data(INPUT_FILE, TARGET_DATE, VAR_NAME)
253+
254+
# Plot
255+
plot_global_spei(da, OUTPUT_FILE)
256+
257+
print("\nDone!")
258+
259+
260+
if __name__ == '__main__':
261+
main()

‎tests/plot_global_spi.py‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -38,10 +38,10 @@
3838
# Input / output paths
3939
BASE_DIR = Path(__file__).parent.parent
4040
INPUT_FILE = BASE_DIR / 'global' / 'output' / 'netcdf' / 'wld_cli_chirps3_d3_spi_gamma_12_month.nc'
41-
OUTPUT_FILE = BASE_DIR / 'docs' / 'images' / 'global-spi12-202512.png'
41+
OUTPUT_FILE = BASE_DIR / 'docs' / 'images' / 'global-spi12-202412.png'
4242

4343
# Time slice to plot
44-
TARGET_DATE = '2025-12-21'
44+
TARGET_DATE = '2024-12-21'
4545

4646
# Variable name in the NetCDF
4747
VAR_NAME = 'spi_gamma_12_month'
@@ -234,7 +234,7 @@ def plot_global_spi(da: xr.DataArray, output_path: Path):
234234
def main():
235235
"""Main entry point."""
236236
print("=" * 60)
237-
print(" GLOBAL SPI-12 MAP — DECEMBER 2025")
237+
print(" GLOBAL SPI-12 MAP — DECEMBER 2024")
238238
print("=" * 60)
239239

240240
# Check input exists

0 commit comments

Comments
 (0)