Skip to content

special

download_era5_uncertainty

ERA5 Uncertainty Download (CLIMATE-22 gap-fill)

Downloads whole-year ERA5 reanalysis and ensemble_spread 2m temperature files from the Copernicus CDS, matching the layout Katrin Burkart produced by hand in /ihme/erf/ERA5/katrin_download_20240208/ for 2022/2023:

era5_{product_type}_{variable}_{year}.nc   (one file per whole year)

This fills the 2024/2025 gap that GBD temperature PAF work depends on. It is a standalone gap-fill: the ensemble_spread product is not (yet) consumed by the rest of this repo's pipeline, which only uses reanalysis mean temperature.

Notes / gotchas baked in below: * The ensemble products only exist on reanalysis-era5-single-levels (0.25° HRES / ~0.5° EDA), never on reanalysis-era5-land. * Both products are pulled 3-hourly to match the canonical GBD download and Katrin's 2022/2023 files. ensemble_spread (the EDA) is only 3-hourly; reanalysis (HRES) is subsampled from hourly to the same 3-hourly axis. * Credentials come from the caller's ~/.cdsapirc (cdsapi.Client() with no args), so no dependency on the shared per-user copernicus.yaml keyring (which does not carry every user).

download_era5_uncertainty(era5_variable: str, output_dir: str, queue: str, dry_run: bool) -> None

Download 2024/2025 ERA5 reanalysis + ensemble_spread (CLIMATE-22 gap-fill).

Source code in src/climate_data/special/download_era5_uncertainty.py 
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
@click.command()
@clio.with_era5_variable()
@clio.with_output_directory(DEFAULT_OUTPUT_DIR)
@clio.with_queue()
@clio.with_dry_run()
def download_era5_uncertainty(
    era5_variable: str,
    output_dir: str,
    queue: str,
    dry_run: bool,
) -> None:
    """Download 2024/2025 ERA5 reanalysis + ensemble_spread (CLIMATE-22 gap-fill)."""
    node_args = {
        "year": ERA5_UNCERTAINTY_YEARS,
        "product-type": ERA5_PRODUCT_TYPES,
    }

    run_parallel_maybe_dry_run(
        runner="cdtask special",
        task_name="download_era5_uncertainty",
        node_args=node_args,
        task_args={
            "era5-variable": era5_variable,
            "output-dir": output_dir,
        },
        task_resources={
            "queue": queue,
            "cores": 1,
            "memory": "20G",
            "runtime": "720m",
            "project": "proj_rapidresponse",
        },
        max_attempts=1,
        dry_run=dry_run,
    )

download_era5_uncertainty_task(year: str, product_type: str, era5_variable: str, output_dir: str) -> None

Download one whole-year ERA5 file (single year, single product type).

Source code in src/climate_data/special/download_era5_uncertainty.py 
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
@click.command()
@clio.with_year(years=cdc.HISTORY_YEARS)
@click.option(
    "--product-type",
    required=True,
    type=click.Choice(ERA5_PRODUCT_TYPES),
    help="ERA5 product type to download.",
)
@clio.with_era5_variable()
@clio.with_output_directory(DEFAULT_OUTPUT_DIR)
def download_era5_uncertainty_task(
    year: str,
    product_type: str,
    era5_variable: str,
    output_dir: str,
) -> None:
    """Download one whole-year ERA5 file (single year, single product type)."""
    download_era5_uncertainty_main(year, product_type, era5_variable, output_dir)

temperature_person_days

temperature_person_days_main(block_key: str, gcm_member: str, scenario: str, hierarchy: str, population_model_root: str, climate_data_root: str, output_dir: str, *, progress_bar: bool = False) -> None

Bin population into (temperature x temperature-zone) person-days per year.

The year span is derived from the temperature_zone actually on disk (not a hardcoded range): for the historical scenario this is the ERA5-only product spanning EXPOSURE_START_YEAR through the last year present (1990-2025 as delivered); for a forecast scenario, daily temperature is ERA5 before FORECAST_START_YEAR and the GCM scenario from then on, binned against that scenario's own zone. A year whose daily or population input is missing fails loudly rather than being skipped, so a gap in a product that must be square cannot slip through silently.

Source code in src/climate_data/special/temperature_person_days.py 
 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
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
def temperature_person_days_main(
    block_key: str,
    gcm_member: str,
    scenario: str,
    hierarchy: str,
    population_model_root: str,
    climate_data_root: str,
    output_dir: str,
    *,
    progress_bar: bool = False,
) -> None:
    """Bin population into (temperature x temperature-zone) person-days per year.

    The year span is derived from the ``temperature_zone`` actually on disk (not a
    hardcoded range): for the ``historical`` scenario this is the ERA5-only product
    spanning ``EXPOSURE_START_YEAR`` through the last year present (1990-2025 as
    delivered); for a forecast scenario, daily temperature is ERA5 before
    ``FORECAST_START_YEAR`` and the GCM scenario from then on, binned against that
    scenario's own zone. A year whose daily or population input is missing fails loudly
    rather than being skipped, so a gap in a product that must be square cannot slip
    through silently.
    """
    print(f"Aggregating {gcm_member} for {block_key}")
    pm_data = PopulationModelData(population_model_root)
    cd_data = ClimateData(climate_data_root, read_only=True)
    ca_data = ClimateAggregateData(Path(output_dir) / hierarchy)

    print("Building location masks")
    climate_slice, location_ids, location_idx = utils.build_location_index(
        hierarchy, block_key, pm_data
    )

    print("Building data index")
    temperature_bins = np.arange(-35, 45, 0.1)
    temperature_zone_bins = np.arange(-25, 35, 1)
    data_idx = pd.MultiIndex.from_product(
        [location_ids, temperature_bins, temperature_zone_bins]
    )
    out_template = np.zeros(
        (len(location_ids), len(temperature_bins), len(temperature_zone_bins)),
        dtype=np.float64,
    )

    print("Building historical temperature zone index")
    temperature_zone = cd_data.load_compiled_annual_results(
        scenario, "temperature_zone", gcm_member
    ).sel(**climate_slice)  # type: ignore[arg-type]
    historical_temperature_zone_idx = utils.to_idx(
        temperature_zone, temperature_zone_bins
    )

    print("Building temperature coordinates")
    temperature_coordinates = utils.get_temperature_coordinates(
        block_key, pm_data, temperature_zone
    )

    print("Aggregating temperature person days")
    # Drive the span from the zone actually on disk (historical: 1990..last present;
    # forecast: the compiled historical+forecast series). A year whose daily or
    # population input is missing now fails loudly rather than being skipped, so a gap
    # in a product that must be square can't slip through silently.
    zone_years = [int(y) for y in temperature_zone["year"].to_numpy()]
    zone_row_for_year = {y: i for i, y in enumerate(zone_years)}
    years = [y for y in zone_years if y >= cdc.EXPOSURE_START_YEAR]
    dfs = []
    for year in tqdm.tqdm(years, disable=not progress_bar):
        if scenario == "historical" or year < FORECAST_START_YEAR:
            temperature = cd_data.load_daily_results(
                "historical", "mean_temperature", year
            ).sel(**climate_slice)  # type: ignore[arg-type]
        else:
            temperature = cd_data.load_raw_daily_results(
                scenario, "mean_temperature", year, gcm_member
            ).sel(**climate_slice)  # type: ignore[arg-type]
        # Population nodata is stored as NaN. The aggregate stage tolerates this via
        # np.nansum, but here we accumulate into `out_arr` with `+=`, so a NaN pixel
        # would poison every output cell it touches (silently read back as 0, which
        # zeroed small locations like American Samoa). Zero-fill nodata first: no
        # modeled population contributes no person-days.
        pop_arr = pm_data.load_results(f"{year}q1", block_key)._ndarray.flatten()  # noqa: SLF001
        pop_arr = np.nan_to_num(pop_arr, nan=0.0)
        temperature_idx = utils.to_idx(temperature, temperature_bins)

        out_arr = out_template.copy()
        utils.compute_person_days(
            location_idx,
            temperature_idx,
            historical_temperature_zone_idx[zone_row_for_year[year]],
            pop_arr,
            temperature_coordinates,
            out_arr,
        )

        df = (
            pd.DataFrame({"person_days": out_arr.reshape(-1)}, index=data_idx)
            .assign(year=year)
            .set_index("year", append=True)
        )
        df.index.names = ["location_id", "temperature", "temperature_zone", "year"]
        df = df.reset_index()
        df["temperature"] = df["temperature"].round(1)
        df = df.set_index(["location_id", "year", "temperature_zone", "temperature"])[
            "person_days"
        ].unstack()
        df.columns.name = None
        dfs.append(df)

    if not dfs:
        msg = (
            f"No person-days produced for {scenario} {gcm_member} {block_key}: the "
            f"temperature zone has no year >= {cdc.EXPOSURE_START_YEAR}."
        )
        raise ValueError(msg)
    df = pd.concat(dfs)
    out_path = ca_data.person_days_path(block_key, scenario, gcm_member)
    mkdir(out_path.parent, parents=True, exist_ok=True)
    save_parquet(df, out_path)

temperature_zone

generate_temperature_zone_main(gcm_member: str, scenario: str, output_dir: str | Path) -> None

Generate the temperature zone for a given scenario and gcm member.

Parameters

gcm_member The gcm member to generate the temperature zone for. scenario The scenario to generate the temperature zone for. Pass historical (with gcm_member="era5") to build a pure-ERA5 zone from the raw historical annual mean temperature rather than the compiled historical+forecast series; the output then spans EXPOSURE_START_YEAR through the last historical year present on disk. output_dir The directory to save the temperature zone to (root for this run mode).

Source code in src/climate_data/special/temperature_zone.py 
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
51
52
53
54
55
def generate_temperature_zone_main(
    gcm_member: str,
    scenario: str,
    output_dir: str | Path,
) -> None:
    """Generate the temperature zone for a given scenario and gcm member.

    Parameters
    ----------
    gcm_member
        The gcm member to generate the temperature zone for.
    scenario
        The scenario to generate the temperature zone for.  Pass ``historical``
        (with ``gcm_member="era5"``) to build a pure-ERA5 zone from the raw
        historical annual mean temperature rather than the compiled
        historical+forecast series; the output then spans ``EXPOSURE_START_YEAR``
        through the last historical year present on disk.
    output_dir
        The directory to save the temperature zone to (root for this run mode).
    """
    print(f"Generating temperature zone for {scenario} {gcm_member}")
    cdata = ClimateData(output_dir)
    if scenario == "historical":
        ds = cdata.load_raw_annual_mfdataset("historical", "mean_temperature")
    else:
        ds = cdata.load_compiled_annual_results(
            scenario, "mean_temperature", gcm_member
        )
    temperature_zone = (
        ds.rolling(year=10)
        .mean()
        .sel(year=slice(cdc.EXPOSURE_START_YEAR, int(cdc.FORECAST_YEARS[-1])))
    )
    print(f"Saving temperature zone for {scenario} {gcm_member}")
    cdata.save_compiled_annual_results(
        temperature_zone,
        scenario=scenario,
        variable="temperature_zone",
        gcm_member=gcm_member,
        encoding_kwargs={"scale_factor": 0.01, "add_offset": 0.0},
    )

utils

aggregate_to_hierarchy(data: pd.DataFrame, hierarchy: pd.DataFrame) -> pd.DataFrame

Create all aggregate climate values for a given hierarchy from most-detailed data.

Parameters

data The most-detailed climate data to aggregate. hierarchy The hierarchy to aggregate the data to.

Returns

pd.DataFrame The climate data with values for all levels of the hierarchy.

Source code in src/climate_data/special/utils.py 
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
def aggregate_to_hierarchy(data: pd.DataFrame, hierarchy: pd.DataFrame) -> pd.DataFrame:
    """Create all aggregate climate values for a given hierarchy from most-detailed data.

    Parameters
    ----------
    data
        The most-detailed climate data to aggregate.
    hierarchy
        The hierarchy to aggregate the data to.

    Returns
    -------
    pd.DataFrame
        The climate data with values for all levels of the hierarchy.
    """
    agg_cols = sorted(
        set(data.columns) - {"location_id", "year_id", "temperature_zone"}
    )

    results = data.set_index("location_id")

    # Most detailed locations can be at multiple levels of the hierarchy,
    # so we loop over all levels from most detailed to global, aggregating
    # level by level and appending the results to the data.

    for level in reversed(list(range(1, hierarchy.level.max() + 1))):
        level_mask = hierarchy.level == level
        parent_map = hierarchy.loc[level_mask].set_index("location_id").parent_id

        subset = results.loc[results.index.intersection(parent_map.index)]
        subset["parent_id"] = parent_map

        parent_values = (
            subset.groupby(["year_id", "parent_id", "temperature_zone"])[agg_cols]
            .sum()
            .reset_index()
            .rename(columns={"parent_id": "location_id"})
            .set_index("location_id")
        )
        results = pd.concat([results, parent_values])
    results = (
        results.reset_index()
        .sort_values(["location_id", "year_id", "temperature_zone"])
        .reset_index(drop=True)
    )

    return results