Loading Seasonal Forecasts from CDS¶

This notebook demonstrates how to load seasonal forecasts from the Copernicus Climate Data Store (CDS):

Connect to the CDS API
Request seasonal forecast data for a specific country and time period
Visualize the data and compute zonal statistics

Note: For training purposes, we only load one lead time (24h) and one year to minimize runtime and data size. In practice, you can fetch all lead times and multiple years.

Additionally, you could compute Standardized Precipitation Index (SPI) or drought probabilities from the downloaded data.

1. Environment Setup¶

You only need to run this cell once. If the library is already installed in your environment, you can skip this step.

Install required packages. Ensure cdsapi is available for CDS API access.

In [ ]:

Copied!

!pixi add cdsapi
!pixi add cdsapi

2. Import Libraries¶

We use cdsapi for data retrieval, xarray for handling NetCDF files, and odc.geo for geospatial operations.

In [7]:

Copied!





import os
if os.getcwd().split("\\")[-1] != "anticipatory-action":
    os.chdir("..")
print(os.getcwd())

import cdsapi
import xarray as xr
import pandas as pd
from odc.geo.xr import xr_reproject
from config.params import Params
from hip.analysis.aoi.analysis_area import AnalysisArea
import os
if os.getcwd().split("\\")[-1] != "anticipatory-action":
    os.chdir("..")
print(os.getcwd())

import cdsapi
import xarray as xr
import pandas as pd
from odc.geo.xr import xr_reproject
from config.params import Params
from hip.analysis.aoi.analysis_area import AnalysisArea

C:\Users\amine.barkaoui\OneDrive - World Food Programme\Documents\GitHub\anticipatory-action

3. Define Parameters¶

country: ISO code of the country
issue: Month of forecast issuance
start_year, end_year: Time range (we use only one year for speed)
data_path, output_path: Paths for input/output

In [15]:

Copied!





country = "ISO"  # Replace with ISO code of the country
issue = 12       # Forecast issue month (December)
start_year = 1981
end_year = 2025  # For training, we only use the last year

data_path = "."  # anticipatory-action directory
output_path = "."
country = "ISO"  # Replace with ISO code of the country
issue = 12       # Forecast issue month (December)
start_year = 1981
end_year = 2025  # For training, we only use the last year

data_path = "."  # anticipatory-action directory
output_path = "."

4. Define Area of Interest (AOI)¶

We use AnalysisArea to get administrative boundaries and resolution.

In [3]:

Copied!





area = AnalysisArea.from_admin_boundaries(
    iso3=country,
    admin_level=0,
    resolution=0.25,
)
shp = area.get_dataset([area.BASE_AREA_DATASET])
area = AnalysisArea.from_admin_boundaries(
    iso3=country,
    admin_level=0,
    resolution=0.25,
)
shp = area.get_dataset([area.BASE_AREA_DATASET])

5. Connect to CDS API¶

Initialize the CDS API client.

In [4]:

Copied!

client = cdsapi.Client()
client = cdsapi.Client()

2025-12-15 15:04:19,744 INFO [2025-12-03T00:00:00Z] To improve our C3S service, we need to hear from you! Please complete this very short [survey](https://confluence.ecmwf.int/x/E7uBEQ/). Thank you.
INFO:ecmwf.datastores.legacy_client:[2025-12-03T00:00:00Z] To improve our C3S service, we need to hear from you! Please complete this very short [survey](https://confluence.ecmwf.int/x/E7uBEQ/). Thank you.

6. Prepare Data Request¶

We request seasonal-original-single-levels dataset from ECMWF system 51.

Important:

We only request one lead time (24h) and one year for speed.
You can extend this to all lead times and years.

In [8]:

Copied!





dataset = "seasonal-original-single-levels"
request = {
    "originating_centre": "ecmwf",
    "system": "51",
    "variable": ["total_precipitation"],
    "year": end_year,  # Uncomment this to load all years: [str(y) for y in range(start_year, end_year + 1)]
    "month": [str(issue)],
    "day": ["01"],
    "leadtime_hour": "24",  # Could be [str(lt) for lt in range(24, 5184, 24)] for all leadtimes
    "data_format": "netcdf",
    "area": [area.bbox[-1], area.bbox[0], area.bbox[1], area.bbox[2]]
}
target = f"{data_path}/data/zmb/zarr/2022/{issue}/ecmwf_cds_24h.nc"
dataset = "seasonal-original-single-levels"
request = {
    "originating_centre": "ecmwf",
    "system": "51",
    "variable": ["total_precipitation"],
    "year": end_year,  # Uncomment this to load all years: [str(y) for y in range(start_year, end_year + 1)]
    "month": [str(issue)],
    "day": ["01"],
    "leadtime_hour": "24",  # Could be [str(lt) for lt in range(24, 5184, 24)] for all leadtimes
    "data_format": "netcdf",
    "area": [area.bbox[-1], area.bbox[0], area.bbox[1], area.bbox[2]]
}
target = f"{data_path}/data/zmb/zarr/2022/{issue}/ecmwf_cds_24h.nc"

7. Retrieve Data¶

This step downloads the NetCDF file from CDS in {data_path}/data/zmb/zarr/{end_year}/{issue}/ecmwf_cds_24h.nc

In [9]:

Copied!

client.retrieve(dataset, request, target)
client.retrieve(dataset, request, target)

2025-12-15 15:05:59,332 INFO Request ID is ba89a364-9850-42bb-9bca-82bdd5a268f3
INFO:ecmwf.datastores.legacy_client:Request ID is ba89a364-9850-42bb-9bca-82bdd5a268f3
2025-12-15 15:05:59,400 INFO status has been updated to accepted
INFO:ecmwf.datastores.legacy_client:status has been updated to accepted
2025-12-15 15:06:09,240 INFO status has been updated to successful
INFO:ecmwf.datastores.legacy_client:status has been updated to successful

Out[9]:

'./data/zmb/zarr/2022/12/ecmwf_cds_24h.nc'

8. Load and Inspect Data¶

Use xarray to open the NetCDF file and attach CRS.

In [10]:

Copied!

da = xr.open_dataset(target).tp.rio.write_crs('epsg:4326')
date = pd.to_datetime(da.forecast_reference_time.values[0]).strftime('%Y-%m-%d')
da = xr.open_dataset(target).tp.rio.write_crs('epsg:4326')
date = pd.to_datetime(da.forecast_reference_time.values[0]).strftime('%Y-%m-%d')

C:\Users\amine.barkaoui\AppData\Local\Temp\ipykernel_20968\1413521792.py:1: FutureWarning: In a future version of xarray decode_timedelta will default to False rather than None. To silence this warning, set decode_timedelta to True, False, or a 'CFTimedeltaCoder' instance.
  da = xr.open_dataset(target).tp.rio.write_crs('epsg:4326')

In [14]:

Copied!

da
da

Out[14]:

<xarray.DataArray 'tp' (number: 51, forecast_reference_time: 1,
                        forecast_period: 1, latitude: 10, longitude: 12)> Size: 24kB
array([[[[[0.001673, ..., 0.001503],
          ...,
          [0.000623, ..., 0.004621]]]],



       ...,



       [[[[0.003464, ..., 0.000481],
          ...,
          [0.000352, ..., 0.003031]]]]], dtype=float32)
Coordinates:
  * number                   (number) int64 408B 0 1 2 3 4 5 ... 46 47 48 49 50
  * forecast_reference_time  (forecast_reference_time) datetime64[ns] 8B 2025...
  * forecast_period          (forecast_period) timedelta64[ns] 8B 1 days
  * latitude                 (latitude) float64 80B -9.079 -10.08 ... -18.08
  * longitude                (longitude) float64 96B 22.0 23.0 ... 32.0 32.99
    valid_time               datetime64[ns] 8B ...
    spatial_ref              int32 4B 0
Attributes: (12/33)
    GRIB_paramId:                             228
    GRIB_dataType:                            fc
    GRIB_numberOfPoints:                      120
    GRIB_typeOfLevel:                         surface
    GRIB_stepUnits:                           1
    GRIB_stepType:                            instant
    ...                                       ...
    GRIB_totalNumber:                         0
    GRIB_units:                               m
    long_name:                                Total precipitation
    units:                                    m
    standard_name:                            unknown
    GRIB_surface:                             0.0

xarray.DataArray

'tp'

number: 51
forecast_reference_time: 1
forecast_period: 1
latitude: 10
longitude: 12

0.001673 0.001986 0.003318 0.003499 ... 0.00193 0.002291 0.003031

array([[[[[0.001673, ..., 0.001503],
          ...,
          [0.000623, ..., 0.004621]]]],



       ...,



       [[[[0.003464, ..., 0.000481],
          ...,
          [0.000352, ..., 0.003031]]]]], dtype=float32)

Coordinates: (7)
- number
  (number)
  int64
  0 1 2 3 4 5 6 ... 45 46 47 48 49 50
  long_name :
  ensemble member numerical id
  units :
  1
  standard_name :
  realization
```
array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
       18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,
       36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50], dtype=int64)
```
- forecast_reference_time
  (forecast_reference_time)
  datetime64[ns]
  2025-12-01
  long_name :
  initial time of forecast
  standard_name :
  forecast_reference_time
```
array(['2025-12-01T00:00:00.000000000'], dtype='datetime64[ns]')
```
- forecast_period
  (forecast_period)
  timedelta64[ns]
  1 days
  long_name :
  time since forecast_reference_time
  standard_name :
  forecast_period
  dtype :
  timedelta64[ns]
```
array([86400000000000], dtype='timedelta64[ns]')
```
- latitude
  (latitude)
  float64
  -9.079 -10.08 ... -17.08 -18.08
  units :
  degrees_north
  standard_name :
  latitude
  long_name :
  latitude
  stored_direction :
  decreasing
```
array([ -9.079, -10.079, -11.079, -12.079, -13.079, -14.079, -15.079, -16.079,
       -17.079, -18.08 ])
```
- longitude
  (longitude)
  float64
  22.0 23.0 24.0 ... 31.0 32.0 32.99
  units :
  degrees_east
  standard_name :
  longitude
  long_name :
  longitude
```
array([21.995, 22.995, 23.995, 24.995, 25.995, 26.995, 27.995, 28.995, 29.995,
       30.995, 31.995, 32.995])
```
- valid_time
  ()
  datetime64[ns]
  ...
  standard_name :
  time
  long_name :
  time
```
[1 values with dtype=datetime64[ns]]
```
- spatial_ref
  ()
  int32
  0
  crs_wkt :
  GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]
  semi_major_axis :
  6378137.0
  semi_minor_axis :
  6356752.314245179
  inverse_flattening :
  298.257223563
  reference_ellipsoid_name :
  WGS 84
  longitude_of_prime_meridian :
  0.0
  prime_meridian_name :
  Greenwich
  geographic_crs_name :
  WGS 84
  horizontal_datum_name :
  World Geodetic System 1984
  grid_mapping_name :
  latitude_longitude
  spatial_ref :
  GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]
```
array(0)
```

Indexes: (5)

number

PandasIndex

PandasIndex(Index([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
       18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,
       36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50],
      dtype='int64', name='number'))

forecast_reference_time

PandasIndex

PandasIndex(DatetimeIndex(['2025-12-01'], dtype='datetime64[ns]', name='forecast_reference_time', freq=None))

forecast_period

PandasIndex

PandasIndex(TimedeltaIndex(['1 days'], dtype='timedelta64[ns]', name='forecast_period', freq=None))

latitude

PandasIndex

PandasIndex(Index([ -9.079, -10.079, -11.079, -12.079, -13.079, -14.079, -15.079, -16.079,
       -17.079,  -18.08],
      dtype='float64', name='latitude'))

longitude

PandasIndex

PandasIndex(Index([21.995, 22.995, 23.995, 24.995, 25.995, 26.995, 27.995, 28.995, 29.995,
       30.995, 31.995, 32.995],
      dtype='float64', name='longitude'))

Attributes: (33)
GRIB_paramId :
228
GRIB_dataType :
fc
GRIB_numberOfPoints :
120
GRIB_typeOfLevel :
surface
GRIB_stepUnits :
1
GRIB_stepType :
instant
GRIB_gridType :
regular_ll
GRIB_uvRelativeToGrid :
0
GRIB_NV :
0
GRIB_Nx :
12
GRIB_Ny :
10
GRIB_cfName :
unknown
GRIB_cfVarName :
tp
GRIB_gridDefinitionDescription :
Latitude/Longitude Grid
GRIB_iDirectionIncrementInDegrees :
1.0
GRIB_iScansNegatively :
0
GRIB_jDirectionIncrementInDegrees :
1.0
GRIB_jPointsAreConsecutive :
0
GRIB_jScansPositively :
0
GRIB_latitudeOfFirstGridPointInDegrees :
-9.079
GRIB_latitudeOfLastGridPointInDegrees :
-18.08
GRIB_longitudeOfFirstGridPointInDegrees :
21.995
GRIB_longitudeOfLastGridPointInDegrees :
32.995
GRIB_missingValue :
3.4028234663852886e+38
GRIB_name :
Total precipitation
GRIB_shortName :
tp
GRIB_system :
51
GRIB_totalNumber :
0
GRIB_units :
m
long_name :
Total precipitation
units :
m
standard_name :
unknown
GRIB_surface :
0.0

9. Visualize Forecast¶

Display the forecast for the AOI at different resolutions.

In [11]:

Copied!

da.isel(forecast_reference_time=0, forecast_period=0, number=0).rio.write_crs('EPSG:4326').rio.clip(shp.geometry).hip.viz.map(
    title=f"{date} - 24h leadtime - 1 degree"
)
da.isel(forecast_reference_time=0, forecast_period=0, number=0).rio.write_crs('EPSG:4326').rio.clip(shp.geometry).hip.viz.map(
    title=f"{date} - 24h leadtime - 1 degree"
)

Out[11]:

No description has been provided for this image

The following map shows the forecast at a 0.25 degree resolution. This forecast is obtained by reprojecting the native one using the bilinear interpolation method.

In [12]:

Copied!

xr_reproject(da, area.geobox, resampling='bilinear').isel(forecast_reference_time=0, forecast_period=0, number=0).dropna('latitude', how='all').dropna('longitude', how='all').rio.clip(shp.geometry).hip.viz.map(
    title=f"{date} - 24h leadtime - 0.25 degree"
)
xr_reproject(da, area.geobox, resampling='bilinear').isel(forecast_reference_time=0, forecast_period=0, number=0).dropna('latitude', how='all').dropna('longitude', how='all').rio.clip(shp.geometry).hip.viz.map(
    title=f"{date} - 24h leadtime - 0.25 degree"
)

Out[12]:

10. Compute Zonal Statistics¶

Aggregate forecast values by administrative boundaries (admin-2 level).

In [13]:

Copied!





area = AnalysisArea.from_admin_boundaries(iso3=country, admin_level=2, resolution=0.25)
admin_fc = area.zonal_stats(da.isel(forecast_reference_time=0, forecast_period=0, number=0), stats=['mean']).reset_index().assign(time=date).set_index(['zone', 'time'])
zonal_gdf = area.join_zonal_stats(admin_fc['mean'])
zonal_gdf.hip.viz.map(title=f"{date} - 24h leadtime - admin 2", column=date, annotate='Name', legend=True)
area = AnalysisArea.from_admin_boundaries(iso3=country, admin_level=2, resolution=0.25)
admin_fc = area.zonal_stats(da.isel(forecast_reference_time=0, forecast_period=0, number=0), stats=['mean']).reset_index().assign(time=date).set_index(['zone', 'time'])
zonal_gdf = area.join_zonal_stats(admin_fc['mean'])
zonal_gdf.hip.viz.map(title=f"{date} - 24h leadtime - admin 2", column=date, annotate='Name', legend=True)

WARNING:root:Mismatch between zone ids and unique zones detected in the data. zone ids length: 103, unique zones detected: 51. zone_ids will be adjusted to match unique zone values.

Out[13]:

Next Steps¶

Fetch all lead times and years for full analysis
Adjust the format: convert cumulative rainfall to daily amounts
Derive SPI and drought probabilities from complete dataset