Notebook 01: Data Acquisition

This notebook downloads and prepares the two core precipitation datasets required for the hybrid bias-correction workflow:

IMERG Late Run (GPM_3IMERGDL v07) – half-hourly satellite precipitation aggregated to daily totals at 0.1-degree resolution, accessed via the NASA GES DISC THREDDS NetCDF Subset Service.
CPC Unified Gauge-Based Analysis (CPC-UNI) – daily gauge-interpolated precipitation at 0.5-degree resolution from NOAA PSL, regridded to 0.1 degrees to match the IMERG grid.

Optionally, we also download the IMERG Final Run (GPM_3IMERGDF v07) so that the bias-corrected Late Run product can later be compared against the calibrated Final Run.

All downloads are spatially subset to the Area of Interest (AOI) defined in Notebook 00, keeping data volumes manageable even for multi-decade time spans.

1 Connect Google Drive (Colab only)

This section is only required when running in Google Colab and your project/data are stored in Google Drive.

Mounting Drive makes your repository and datasets accessible under /content/drive.
If you run this notebook locally (Jupyter / VS Code), skip this section.

Expected structure (Drive) After mounting, your project root should contain: - notebooks/ - src/ - config.yml (or config.yaml)

Proceed to the code cell below to mount Drive.

from google.colab import drive
import os

# Check if the drive is mounted
if os.path.exists('/content/drive'):
    # Try to unmount
    try:
        drive.flush_and_unmount()
        print("Successfully unmounted")
    except:
        print("Unmount failed, the drive might not be mounted or busy")

# Mount the drive
drive.mount('/content/drive')

Successfully unmounted
Mounted at /content/drive

Troubleshooting:
- If Colab becomes disconnected, Reconnect the runtime and rerun the mounting cell. - If we receive an error such as Mountpoint must not already contain files, delete all the sub-folders under “/content/drive” from the Files panel before retrying. We need to delete these one by one starting from the innermost folders, until the last “drive” folder is deleted.

2 Install packages (only if needed)

In most cases, Google Colab already includes the packages required for this workflow. The most common missing dependency is netCDF4 (NetCDF I/O support).

Check what is already installed (Colab)

Before installing anything, you can inspect the current environment by running !pip list.

If all required packages are present and only netCDF4 is missing, install only netCDF4.
If other required packages are missing from !pip list, install them together with netCDF4 in the code cell below.

Local Jupyter note

If you are running in a local environment (Jupyter / VS Code), assume all dependencies were installed when preparing the environment following the main repository README. In that case, you can skip this section.

Proceed to the code cell below only when installation is necessary.

# In Google Colab, almost all packages already available, except netCDF4
!pip install netCDF4 rasterio cartopy

Requirement already satisfied: netCDF4 in /usr/local/lib/python3.12/dist-packages (1.7.4)
Requirement already satisfied: rasterio in /usr/local/lib/python3.12/dist-packages (1.5.0)
Requirement already satisfied: cartopy in /usr/local/lib/python3.12/dist-packages (0.25.0)
Requirement already satisfied: cftime in /usr/local/lib/python3.12/dist-packages (from netCDF4) (1.6.5)
Requirement already satisfied: certifi in /usr/local/lib/python3.12/dist-packages (from netCDF4) (2026.2.25)
Requirement already satisfied: numpy>=1.21.2 in /usr/local/lib/python3.12/dist-packages (from netCDF4) (2.0.2)
Requirement already satisfied: affine in /usr/local/lib/python3.12/dist-packages (from rasterio) (2.4.0)
Requirement already satisfied: attrs in /usr/local/lib/python3.12/dist-packages (from rasterio) (25.4.0)
Requirement already satisfied: click!=8.2.*,>=4.0 in /usr/local/lib/python3.12/dist-packages (from rasterio) (8.3.1)
Requirement already satisfied: cligj>=0.5 in /usr/local/lib/python3.12/dist-packages (from rasterio) (0.7.2)
Requirement already satisfied: pyparsing in /usr/local/lib/python3.12/dist-packages (from rasterio) (3.3.2)
Requirement already satisfied: matplotlib>=3.6 in /usr/local/lib/python3.12/dist-packages (from cartopy) (3.10.0)
Requirement already satisfied: shapely>=2.0 in /usr/local/lib/python3.12/dist-packages (from cartopy) (2.1.2)
Requirement already satisfied: packaging>=21 in /usr/local/lib/python3.12/dist-packages (from cartopy) (26.0)
Requirement already satisfied: pyshp>=2.3 in /usr/local/lib/python3.12/dist-packages (from cartopy) (3.0.3)
Requirement already satisfied: pyproj>=3.3.1 in /usr/local/lib/python3.12/dist-packages (from cartopy) (3.7.2)
Requirement already satisfied: contourpy>=1.0.1 in /usr/local/lib/python3.12/dist-packages (from matplotlib>=3.6->cartopy) (1.3.3)
Requirement already satisfied: cycler>=0.10 in /usr/local/lib/python3.12/dist-packages (from matplotlib>=3.6->cartopy) (0.12.1)
Requirement already satisfied: fonttools>=4.22.0 in /usr/local/lib/python3.12/dist-packages (from matplotlib>=3.6->cartopy) (4.61.1)
Requirement already satisfied: kiwisolver>=1.3.1 in /usr/local/lib/python3.12/dist-packages (from matplotlib>=3.6->cartopy) (1.4.9)
Requirement already satisfied: pillow>=8 in /usr/local/lib/python3.12/dist-packages (from matplotlib>=3.6->cartopy) (11.3.0)
Requirement already satisfied: python-dateutil>=2.7 in /usr/local/lib/python3.12/dist-packages (from matplotlib>=3.6->cartopy) (2.9.0.post0)
Requirement already satisfied: six>=1.5 in /usr/local/lib/python3.12/dist-packages (from python-dateutil>=2.7->matplotlib>=3.6->cartopy) (1.17.0)

3 Project Folders

Next, we establish a clear folder structure on your mounted Drive so that raw inputs, downloads, intermediate masks, and final outputs each live in their own dedicated directory. The snippet below:

Declares a single main project folder
Builds subpaths for input data, downloaded files, mask outputs, and final results
Ensures each required directory exists, creating it if necessary

A well-organized directory layout makes it easy to trace where every file came from and prevents accidental overwrites—setting you up for a reproducible, scalable workflow.

# Import library
import os

# Configurable Directory
# Base folder on Google Drive for this project
main_dir = '/content/drive/MyDrive/hybrid-bias-correction'

# Define key subfolders
input_dir     = f'{main_dir}/data/input'      # RAW input data goes here
downloads_dir = f'{main_dir}/data/downloads'  # All downloaded files
mask_dir      = f'{main_dir}/data/mask'       # Generated land masks
output_dir    = f'{main_dir}/data/output'     # Final outputs & results

# List of directories to ensure exist
for d in [input_dir, downloads_dir, mask_dir, output_dir]:
    os.makedirs(d, exist_ok=True)
    print(f"✅ Ensured directory exists: {d}")

✅ Ensured directory exists: /content/drive/MyDrive/hybrid-bias-correction/data/input
✅ Ensured directory exists: /content/drive/MyDrive/hybrid-bias-correction/data/downloads
✅ Ensured directory exists: /content/drive/MyDrive/hybrid-bias-correction/data/mask
✅ Ensured directory exists: /content/drive/MyDrive/hybrid-bias-correction/data/output

What’s next?

With our workspace now structured, we’ll move on to configuring NASA Earthdata authentication so that we can programmatically fetch IMERG and CPC data.

4 Configuring NASA Earthdata Authentication

To programmatically download IMERG and CPC datasets from NASA’s GES DISC THREDDS server, you must authenticate with your Earthdata credentials. The snippet below:

Creates a .netrc file containing your username and password

machine urs.earthdata.nasa.gov login {username} password {password}

Moves that file into /root/.netrc so that both wget and Python’s requests library can automatically pick it up

Without this step, any attempt to access the NCSS endpoints will fail with an authentication error. Make sure you have:

Registered at https://urs.earthdata.nasa.gov/users/new
Approved the “NASA GESDISC DATA ARCHIVE” application via the link provided in your Earthdata account

import os


def create_netrc_file(username: str, password: str):
    """
    Write a .netrc file in the current directory with Earthdata credentials.
    """
    netrc_content = (
        f"machine urs.earthdata.nasa.gov "
        f"login {username} password {password}"
    )
    with open('.netrc', 'w') as f:
        f.write(netrc_content)
    print("✅ .netrc file created in working directory.")

def move_to_colab_root():
    """
    Move the .netrc file to /root so that wget and requests can use it.
    """
    target = '/root/.netrc'
    os.replace('.netrc', target)
    print(f"✅ .netrc file moved to {target}.")

def main():
    try:
        # TODO: replace these with your actual Earthdata credentials
        username = "your_username"
        password = "your_password"

        create_netrc_file(username, password)
        move_to_colab_root()
        print("✅ NASA Earthdata authentication is now configured.")

    except Exception as e:
        print(f"ERROR: ❌ Failed to configure .netrc: {e}")
        print("❌ An error occurred. See log for details.")

if __name__ == "__main__":
    main()

What’s next?

With authentication in place, we can safely issue HTTPS requests to the THREDDS NCSS endpoints (Section 3) and retrieve only the spatial–temporal slices of IMERG or CPC data that you need, without downloading full global archives.

5 Data Acquisition

In this section, we’ll obtain the two core datasets needed for our hybrid bias correction:

IMERG Late Run satellite rainfall estimates from NASA’s GES DISC https://gpm.nasa.gov/data/imerg
CPC-UNI gauge-based daily precipitation from NOAA PSL https://psl.noaa.gov/data/gridded/data.cpc.globalprecip.html

Optionally we could also download

IMERG Final Run - to compare how good the corrected IMERG Late Run from the official (Final) version.

Rather than downloading hundreds of gigabytes of global files (2001–2025), we use NASA’s THREDDS NetCDF Subset Service (NCSS) to pull only the spatial–temporal slice corresponding to our AOI. The THREDDS server exposes annual aggregation files at URLs of the form:

IMERG Late Run

https://gpm1.gesdisc.eosdis.nasa.gov/thredds/ncss/grid/aggregation/GPM_3IMERGDL.07/GPM_3IMERGDL.07_Aggregation_{year}.ncml.ncml

IMERG Final Run

https://gpm1.gesdisc.eosdis.nasa.gov/thredds/ncss/grid/aggregation/GPM_3IMERGDF.07/GPM_3IMERGDF.07_Aggregation_{year}.ncml.ncml

By supplying query parameters for the precipitation variable, north/west/south/east bounds, horizontal stride, and start/end timestamps, we instruct the server to return a much smaller NetCDF4‑classic subset rather than the full global field. This approach drastically reduces download volume and speeds up data access for any arbitrary bounding box.

The following code first computes a buffered bounding box around our idn_subset.nc so that every IMERG request is confined to exactly the region we care about.

We’ll calculate the buffered bounding box from our AOI mask. This script:

Loads the AOI land mask
Identifies the min/max latitude and longitude of all land pixels
Expands each side by a configurable buffer (e.g. 0.5°)
Prints out the north, south, west, and east values for use in subsequent download requests

Once we have these parameters, we’re ready to loop over each year and grab only the subset of IMERG data we need.

# Calculate buffered bounding box from AOI land mask
import os
import numpy as np
import xarray as xr

def load_aoi_mask(path: str) -> xr.DataArray:
    print(f"Loading AOI land mask from {path} ...")
    da = xr.open_dataarray(path)
    print("AOI land mask loaded")
    return da

def calculate_bbox(da: xr.DataArray, buffer: float = 0.5):
    """
    Compute the bounding box of non-NaN mask cells and expand each side
    by `buffer` degrees.  Returns (north, south, west, east).
    """
    lats, lons = np.meshgrid(da['lat'].values, da['lon'].values, indexing='ij')
    mask = ~np.isnan(da.values)
    if not mask.any():
        raise ValueError("AOI mask contains no land pixels.")

    lat_masked = lats[mask]
    lon_masked = lons[mask]

    south = max(lat_masked.min() - buffer, -90.0)
    north = min(lat_masked.max() + buffer,  90.0)
    west  = ((lon_masked.min() - buffer + 180) % 360) - 180
    east  = ((lon_masked.max() + buffer + 180) % 360) - 180

    return north, south, west, east

# --- Run and store as module-level variables ---
mask_path = os.path.join(main_dir, 'data/mask/iso3/idn_subset.nc')
da_aoi = load_aoi_mask(mask_path)
north, south, west, east = calculate_bbox(da_aoi, buffer=0.5)
da_aoi.close()

print(f"\n# ===== Bounding Box (used by all download cells below) =====")
print(f"north = {north}")
print(f"south = {south}")
print(f"west  = {west}")
print(f"east  = {east}")

5.1 Downloading IMERG Late Run (GPM_3IMERGDL)

The following code will download the IMERG-L from the GES DISC using predefined query. With our bounding box defined above, we construct a yearly loop that:

Builds the THREDDS NCSS URL for each calendar year
Supplies the precipitation variable, spatial bounds, and start/end timestamps
Streams a NetCDF4-classic file for each year into our /data/downloads/GPM_3IMERGDL_subset folder

NOTE

To understand correctly how the URL and Query are build, lets access this link NetCDF Subset Service for Grids with Dataset: /thredds/ncss/grid/aggregation/GPM_3IMERGDL.07/GPM_3IMERGDL.07_Aggregation_2000.ncml.ncml then choose and adjust some parameters: - Variables - Horizontal subset - Time subset - Vertical subset - Output format - CF compliance

We will get the full NCSS request URL at the bottom. Those URL are use to build the loop for downloading multi year data.

Step 1: Download process

This approach drastically reduces download volume and ensures that we only fetch the pixels and dates relevant to our AOI.

import requests
from datetime import datetime, timedelta
import os

# ===== User Parameters =====
# We may choose GPM_3IMERGDF for Final or GPM_3IMERGDL for Late Run
product_name = 'GPM_3IMERGDL'

# Bounding box — computed from AOI mask in the cell above
# (north, south, west, east are already defined)
print(f"Using bbox: north={north}, south={south}, west={west}, east={east}")

# Date range (inclusive)
start_date_str = "2001-01-01"
end_date_str   = "2025-12-31"

# Convert to datetime objects
start_date = datetime.strptime(start_date_str, "%Y-%m-%d")
end_date   = datetime.strptime(end_date_str,   "%Y-%m-%d")

# Optional: Horizontal stride (set to 1 for no skipping)
horiz_stride = 1

"""
IMERG v7, the latest version of the Integrated Multi-satellitE Retrievals for GPM,
offers two key variables for analyzing precipitation data from the THREDDS:
 - "precipitation" which is the merged satellite-gauge precipitation estimate.
 - "MWprecipitation" which represents microwave-only precipitation estimates.
"""
# Define the variable name (e.g., "precipitation")
variable = "precipitation"

# Output directory
imerg_download_dir = os.path.join(main_dir, 'data/downloads', f'{product_name}_subset')
os.makedirs(imerg_download_dir, exist_ok=True)

def year_bounds(year, start_date, end_date):
    """
    Return the (t0, t1) window for this calendar year slice.
    """
    t0 = start_date if year == start_date.year else datetime(year, 1, 1)
    t1 = end_date   if year == end_date.year   else datetime(year, 12, 31)
    return t0, t1

for yr in range(start_date.year, end_date.year + 1):
    t0, t1 = year_bounds(yr, start_date, end_date)

    # time_end must stay within the same calendar year's aggregation file;
    # use 23:29:59 on the last day (IMERG half-hourly last slot = 23:00-23:29)
    time_start = t0.strftime("%Y-%m-%dT00:00:00Z")
    time_end   = t1.strftime("%Y-%m-%dT23:29:59Z")

    # Per-year NCML aggregation URL
    base_url = (
        f"https://gpm1.gesdisc.eosdis.nasa.gov/thredds/ncss/"
        f"grid/aggregation/{product_name}.07/"
        f"{product_name}.07_Aggregation_{yr}.ncml.ncml"
    )

    # Build Query parameters
    params = {
        "var":         variable,
        "north":       f"{north:.3f}",
        "west":        f"{west:.3f}",
        "east":        f"{east:.3f}",
        "south":       f"{south:.3f}",
        "addLatLon":   "true",
        "horizStride": horiz_stride,
        "time_start":  time_start,
        "time_end":    time_end,
        "accept":      "netcdf4-classic"
    }

    fname = f"{product_name}_subset_{t0.strftime('%Y%m%d')}_{t1.strftime('%Y%m%d')}.nc4"
    outpath = os.path.join(imerg_download_dir, fname)

    if os.path.exists(outpath):
        print(f"  {yr}: already exists, skipping")
        continue

    print(f"Downloading {yr}: {time_start} -> {time_end} ...")
    r = requests.get(base_url, params=params, stream=True)
    if r.status_code == 200:
        with open(outpath, 'wb') as f:
            for chunk in r.iter_content(1024):
                if chunk:
                    f.write(chunk)
        print(f"  Saved -> {fname}")
    else:
        print(f"  Failed {yr}: HTTP {r.status_code} -> {r.url}")

print("All done.")

You can also use this exact same loop to grab the official IMERG Final Run—simply set product_name = 'GPM_3IMERGDF' (and adjust your output folder to /data/downloads/GPM_3IMERGDF_subset) so that later on we can compare our bias-corrected IMERG-L against the standard Final Run.

Step 2: Extract and Subset Variable

This script extract variable precipitation and clip using mask data, also add compression to the output to reduce file size.

import xarray as xr
import os
import glob
from tqdm import tqdm
import numpy as np

# Define the input and output directories
input_dir = os.path.join(main_dir, 'data/downloads/GPM_3IMERGDL_subset')
extract_dir = os.path.join(main_dir, 'data/downloads/GPM_3IMERGDL_extract')
mask_file = os.path.join(main_dir, 'data/mask/iso3/idn_subset.nc')

# Date range
start_year = 2001
end_year = 2025

os.makedirs(extract_dir, exist_ok=True)

# Load the mask dataset
mask_ds = xr.open_dataset(mask_file)
land_mask = mask_ds['land'].astype(float)

# List all downloaded subset files
file_list = sorted(glob.glob(os.path.join(input_dir, "GPM_3IMERGDL_subset_*.nc4")))
print(f"Found {len(file_list)} downloaded IMERG-L subset files")

# Process each yearly subset: extract precipitation, clip to AOI
for file_path in tqdm(file_list, desc="Processing IMERG-L subsets"):
    # Extract year from filename
    # e.g. GPM_3IMERGDL_subset_20010101_20011231.nc4
    #      [0]  [1]       [2]     [3]       [4]
    basename = os.path.basename(file_path)
    year_str = basename.split('_')[3][:4]
    year = int(year_str)

    if year < start_year or year > end_year:
        continue

    output_filename = f"idn_cli_imerg_ldd_{year}.nc4"
    output_filepath = os.path.join(extract_dir, output_filename)

    if os.path.exists(output_filepath):
        print(f"  {year}: already exists, skipping")
        continue

    try:
        ds = xr.open_dataset(file_path)
        ds = ds[['precipitation']]

        # Subset to mask extent
        ds = ds.sel(
            lat=slice(float(land_mask.lat.min()), float(land_mask.lat.max())),
            lon=slice(float(land_mask.lon.min()), float(land_mask.lon.max()))
        )

        # Align mask and apply
        mask_interp = land_mask.interp(lat=ds.lat, lon=ds.lon, method="nearest")
        ds = ds.where(mask_interp == 1, drop=True)

        # Ensure (time, lat, lon) dimension order
        ds = ds.transpose('time', 'lat', 'lon')

        # Save with compression
        comp = dict(zlib=True, complevel=5)
        encoding = {var: comp for var in ds.data_vars}
        ds.attrs['Conventions'] = 'CF-1.8'
        ds.to_netcdf(output_filepath, encoding=encoding, engine='netcdf4')
        ds.close()
        print(f"  Saved: {output_filename}")

    except Exception as e:
        print(f"  Error processing {year}: {e}")

mask_ds.close()
print("IMERG-L extraction complete.")

Step 3: Combine into a Single Multi-Year NetCDF

The bias-correction pipeline (notebook 02) expects a single NetCDF file containing the full daily time series across all years. This step concatenates the per-year extracted files along the time dimension and saves the result to data/input/.

import xarray as xr
import netCDF4
import os
import glob
import numpy as np

extract_dir = os.path.join(main_dir, 'data/downloads/GPM_3IMERGDL_extract')
output_file = os.path.join(main_dir, 'data/input/idn_cli_imergl_daily.nc4')

os.makedirs(os.path.dirname(output_file), exist_ok=True)

file_list = sorted(glob.glob(os.path.join(extract_dir, "idn_cli_imerg_ldd_*.nc4")))
print(f"Found {len(file_list)} yearly IMERG-L extract files")

if not file_list:
    print("No extract files found. Run Step 2 first.")
else:
    nc_out = None
    total_days = 0

    for i, fpath in enumerate(file_list):
        year = os.path.basename(fpath).split('_')[-1].replace('.nc4', '')
        ds = xr.open_dataset(fpath).sortby('time')
        ndays = len(ds.time)

        if i == 0:
            nc_out = netCDF4.Dataset(output_file, 'w', format='NETCDF4')
            nc_out.Conventions = 'CF-1.8'
            nc_out.title = 'IMERG Late Run daily precipitation subset for AOI'

            nc_out.createDimension('time', None)
            nc_out.createDimension('lat', len(ds.lat))
            nc_out.createDimension('lon', len(ds.lon))

            t_var = nc_out.createVariable('time', 'f8', ('time',))
            t_var.units = 'days since 1900-01-01'
            t_var.calendar = 'standard'
            t_var.standard_name = 'time'

            lat_var = nc_out.createVariable('lat', 'f4', ('lat',))
            lat_var[:] = ds.lat.values.astype('float32')
            lat_var.units = 'degrees_north'
            lat_var.standard_name = 'latitude'
            lat_var.axis = 'Y'

            lon_var = nc_out.createVariable('lon', 'f4', ('lon',))
            lon_var[:] = ds.lon.values.astype('float32')
            lon_var.units = 'degrees_east'
            lon_var.standard_name = 'longitude'
            lon_var.axis = 'X'

            p_var = nc_out.createVariable(
                'precipitation', 'f4', ('time', 'lat', 'lon'),
                zlib=True, complevel=5,
                chunksizes=(1, len(ds.lat), len(ds.lon)),
                fill_value=np.float32(9.96921e+36)
            )
            if 'units' in ds['precipitation'].attrs:
                p_var.units = ds['precipitation'].attrs['units']
            if 'long_name' in ds['precipitation'].attrs:
                p_var.long_name = ds['precipitation'].attrs['long_name']

        time_vals = netCDF4.date2num(
            ds.time.values.astype('datetime64[ms]').astype('O'),
            units='days since 1900-01-01',
            calendar='standard'
        )

        nc_out.variables['time'][total_days:total_days + ndays] = time_vals
        nc_out.variables['precipitation'][total_days:total_days + ndays, :, :] = \
            ds['precipitation'].values
        nc_out.sync()

        total_days += ndays
        print(f"  [{i+1}/{len(file_list)}] {year}: +{ndays} days "
              f"(total: {total_days})")
        ds.close()

    nc_out.close()

    with xr.open_dataset(output_file) as final:
        print(f"\nDone: {output_file}")
        print(f"  Time: {str(final.time.values[0])[:10]} to "
              f"{str(final.time.values[-1])[:10]}")
        print(f"  Shape: {dict(final.dims)}")

5.2 Downloading IMERG Final Run (GPM_3IMERGDF)

The IMERG Final Run incorporates monthly gauge calibration and is therefore considered the research-quality benchmark. Downloading it is optional but recommended: once the Late Run has been bias-corrected (notebook 02), the Final Run provides an independent reference to evaluate how closely the corrected product approaches the official calibrated estimate.

The workflow mirrors Section 5.1 exactly — download yearly subsets via THREDDS NCSS, extract and clip to the AOI, then concatenate into a single multi-year file — with product_name set to GPM_3IMERGDF.

Step 1: Download process

Same THREDDS NCSS approach as Section 5.1, but targeting the Final Run aggregation endpoint (GPM_3IMERGDF.07). Files are saved into data/downloads/GPM_3IMERGDF_subset/.

import requests
from datetime import datetime, timedelta
import os

# ===== User Parameters =====
product_name = 'GPM_3IMERGDF'

# Bounding box — reuses north/south/west/east from the AOI mask cell
print(f"Using bbox: north={north}, south={south}, west={west}, east={east}")

# Date range (inclusive)
start_date_str = "2001-01-01"
end_date_str   = "2025-12-31"

start_date = datetime.strptime(start_date_str, "%Y-%m-%d")
end_date   = datetime.strptime(end_date_str,   "%Y-%m-%d")

horiz_stride = 1
variable = "precipitation"

# Output directory
imergf_download_dir = os.path.join(main_dir, 'data/downloads', f'{product_name}_subset')
os.makedirs(imergf_download_dir, exist_ok=True)

def year_bounds(year, start_date, end_date):
    t0 = start_date if year == start_date.year else datetime(year, 1, 1)
    t1 = end_date   if year == end_date.year   else datetime(year, 12, 31)
    return t0, t1

for yr in range(start_date.year, end_date.year + 1):
    t0, t1 = year_bounds(yr, start_date, end_date)

    # time_end must stay within the same calendar year's aggregation file;
    # use 23:29:59 on the last day (IMERG half-hourly last slot = 23:00-23:29)
    time_start = t0.strftime("%Y-%m-%dT00:00:00Z")
    time_end   = t1.strftime("%Y-%m-%dT23:29:59Z")

    base_url = (
        f"https://gpm1.gesdisc.eosdis.nasa.gov/thredds/ncss/"
        f"grid/aggregation/{product_name}.07/"
        f"{product_name}.07_Aggregation_{yr}.ncml.ncml"
    )

    params = {
        "var":         variable,
        "north":       f"{north:.3f}",
        "west":        f"{west:.3f}",
        "east":        f"{east:.3f}",
        "south":       f"{south:.3f}",
        "addLatLon":   "true",
        "horizStride": horiz_stride,
        "time_start":  time_start,
        "time_end":    time_end,
        "accept":      "netcdf4-classic"
    }

    fname = f"{product_name}_subset_{t0.strftime('%Y%m%d')}_{t1.strftime('%Y%m%d')}.nc4"
    outpath = os.path.join(imergf_download_dir, fname)

    if os.path.exists(outpath):
        print(f"  {yr}: already exists, skipping")
        continue

    print(f"Downloading {yr}: {time_start} -> {time_end} ...")
    r = requests.get(base_url, params=params, stream=True)
    if r.status_code == 200:
        with open(outpath, 'wb') as f:
            for chunk in r.iter_content(1024):
                if chunk:
                    f.write(chunk)
        print(f"  Saved -> {fname}")
    else:
        print(f"  Failed {yr}: HTTP {r.status_code} -> {r.url}")

print("All done.")

Step 2: Extract and Subset Variable

Extract the precipitation variable, clip to the AOI land mask, and save compressed yearly files — identical to Section 5.1 Step 2 but reading from the Final Run download folder.

import xarray as xr
import os
import glob
from tqdm import tqdm
import numpy as np

# Define the input and output directories
input_dir_f = os.path.join(main_dir, 'data/downloads/GPM_3IMERGDF_subset')
extract_dir_f = os.path.join(main_dir, 'data/downloads/GPM_3IMERGDF_extract')
mask_file = os.path.join(main_dir, 'data/mask/iso3/idn_subset.nc')

start_year = 2001
end_year = 2025

os.makedirs(extract_dir_f, exist_ok=True)

# Load the mask dataset
mask_ds = xr.open_dataset(mask_file)
land_mask = mask_ds['land'].astype(float)

# List all downloaded subset files
file_list = sorted(glob.glob(os.path.join(input_dir_f, "GPM_3IMERGDF_subset_*.nc4")))
print(f"Found {len(file_list)} downloaded IMERG-F subset files")

for file_path in tqdm(file_list, desc="Processing IMERG-F subsets"):
    # Extract year from filename
    # e.g. GPM_3IMERGDF_subset_20010101_20011231.nc4
    #      [0]  [1]       [2]     [3]       [4]
    basename = os.path.basename(file_path)
    year_str = basename.split('_')[3][:4]
    year = int(year_str)

    if year < start_year or year > end_year:
        continue

    output_filename = f"idn_cli_imerg_fdd_{year}.nc4"
    output_filepath = os.path.join(extract_dir_f, output_filename)

    if os.path.exists(output_filepath):
        print(f"  {year}: already exists, skipping")
        continue

    try:
        ds = xr.open_dataset(file_path)
        ds = ds[['precipitation']]

        ds = ds.sel(
            lat=slice(float(land_mask.lat.min()), float(land_mask.lat.max())),
            lon=slice(float(land_mask.lon.min()), float(land_mask.lon.max()))
        )

        mask_interp = land_mask.interp(lat=ds.lat, lon=ds.lon, method="nearest")
        ds = ds.where(mask_interp == 1, drop=True)
        ds = ds.transpose('time', 'lat', 'lon')

        comp = dict(zlib=True, complevel=5)
        encoding = {var: comp for var in ds.data_vars}
        ds.attrs['Conventions'] = 'CF-1.8'
        ds.to_netcdf(output_filepath, encoding=encoding, engine='netcdf4')
        ds.close()
        print(f"  Saved: {output_filename}")

    except Exception as e:
        print(f"  Error processing {year}: {e}")

mask_ds.close()
print("IMERG-F extraction complete.")

Step 3: Combine into a Single Multi-Year NetCDF

Concatenate the per-year IMERG Final Run extracts into one file at data/input/idn_cli_imergf_daily.nc4, ready for use as a comparison baseline in downstream notebooks.

import xarray as xr
import netCDF4
import os
import glob
import numpy as np

extract_dir_f = os.path.join(main_dir, 'data/downloads/GPM_3IMERGDF_extract')
output_file_f = os.path.join(main_dir, 'data/input/idn_cli_imergf_daily.nc4')

os.makedirs(os.path.dirname(output_file_f), exist_ok=True)

file_list = sorted(glob.glob(os.path.join(extract_dir_f, "idn_cli_imerg_fdd_*.nc4")))
print(f"Found {len(file_list)} yearly IMERG-F extract files")

if not file_list:
    print("No extract files found. Run Step 2 first.")
else:
    nc_out = None
    total_days = 0

    for i, fpath in enumerate(file_list):
        year = os.path.basename(fpath).split('_')[-1].replace('.nc4', '')
        ds = xr.open_dataset(fpath).sortby('time')
        ndays = len(ds.time)

        if i == 0:
            nc_out = netCDF4.Dataset(output_file_f, 'w', format='NETCDF4')
            nc_out.Conventions = 'CF-1.8'
            nc_out.title = 'IMERG Final Run daily precipitation subset for AOI'

            nc_out.createDimension('time', None)
            nc_out.createDimension('lat', len(ds.lat))
            nc_out.createDimension('lon', len(ds.lon))

            t_var = nc_out.createVariable('time', 'f8', ('time',))
            t_var.units = 'days since 1900-01-01'
            t_var.calendar = 'standard'
            t_var.standard_name = 'time'

            lat_var = nc_out.createVariable('lat', 'f4', ('lat',))
            lat_var[:] = ds.lat.values.astype('float32')
            lat_var.units = 'degrees_north'
            lat_var.standard_name = 'latitude'
            lat_var.axis = 'Y'

            lon_var = nc_out.createVariable('lon', 'f4', ('lon',))
            lon_var[:] = ds.lon.values.astype('float32')
            lon_var.units = 'degrees_east'
            lon_var.standard_name = 'longitude'
            lon_var.axis = 'X'

            p_var = nc_out.createVariable(
                'precipitation', 'f4', ('time', 'lat', 'lon'),
                zlib=True, complevel=5,
                chunksizes=(1, len(ds.lat), len(ds.lon)),
                fill_value=np.float32(9.96921e+36)
            )
            if 'units' in ds['precipitation'].attrs:
                p_var.units = ds['precipitation'].attrs['units']
            if 'long_name' in ds['precipitation'].attrs:
                p_var.long_name = ds['precipitation'].attrs['long_name']

        time_vals = netCDF4.date2num(
            ds.time.values.astype('datetime64[ms]').astype('O'),
            units='days since 1900-01-01',
            calendar='standard'
        )

        nc_out.variables['time'][total_days:total_days + ndays] = time_vals
        nc_out.variables['precipitation'][total_days:total_days + ndays, :, :] = \
            ds['precipitation'].values
        nc_out.sync()

        total_days += ndays
        print(f"  [{i+1}/{len(file_list)}] {year}: +{ndays} days "
              f"(total: {total_days})")
        ds.close()

    nc_out.close()

    with xr.open_dataset(output_file_f) as final:
        print(f"\nDone: {output_file_f}")
        print(f"  Time: {str(final.time.values[0])[:10]} to "
              f"{str(final.time.values[-1])[:10]}")
        print(f"  Shape: {dict(final.dims)}")

5.3 Downloading CPC-UNI

The following code will download the CPC Global Unified Gauge-Based Analysis of Daily Precipitation data provided by the NOAA PSL, Boulder, Colorado, USA, from their website at https://psl.noaa.gov. The scripts are designed to prepare the data for further analysis by extracting the precipitation variable, fill missing value, regrid and clip using mask boundary. The processed data is structured to follow CF conventions, ensuring consistency and compatibility for subsequent climatological and hydrological studies. Each script is self-contained and builds upon the output from the previous step, maintaining a streamlined workflow from raw data acquisition to processed datasets ready for analysis.

If you need to modify the naming convention of the file output, input, and output folders, ensure these changes are consistently applied across all scripts to maintain dependencies on the output from the previous script. This approach ensures each step correctly processes data generated from the preceding step.

Step 1: Download process

This script download the Global CPC-UNI data. It will download data from 2001-2025, 25 files (1 nc file for 1 year daily data) with total size 1.5GB (as of Feb 2026). Please make sure you have unlimited data plan and high speed internet.

import os
import requests

# Define directory
cpcuni_dir = os.path.join(main_dir, 'data/downloads/CPC_UNI')
os.makedirs(cpcuni_dir, exist_ok=True)

# Download CPC-UNI data from 2001-2025
base_url = "https://downloads.psl.noaa.gov/Datasets/cpc_global_precip/precip."

for year in range(2001, 2026):
    url = f"{base_url}{year}.nc"
    local_path = os.path.join(cpcuni_dir, f"precip.{year}.nc")

    if os.path.exists(local_path):
        print(f"  {year}: already exists, skipping")
        continue

    print(f"Downloading {year} ...")
    response = requests.get(url)
    if response.status_code == 200:
        with open(local_path, 'wb') as f:
            f.write(response.content)
        print(f"  Saved: precip.{year}.nc")
    else:
        print(f"  Failed {year}: HTTP {response.status_code}")

print("All CPC-UNI files downloaded.")

Step 2: Clip, Fill, Regrid, and Mask

Processing the full global CPC grid (720 × 360 at 0.5°) through regridding to 0.1° is extremely slow and memory-intensive. Instead we apply a four-stage pipeline that operates only on the region of interest:

Coarse clip — subset the raw 0.5° data to the AOI bounding box plus a 1.0° buffer on each side (≥ 2 CPC cells beyond the AOI boundary). This “all-touch” margin ensures that every IMERG pixel has at least two neighbouring CPC cells in each direction, which is required for bilinear interpolation during native-resolution parameter fitting (Notebook 02).
Fill — iteratively spread non-NaN values into coastal/ocean NaN cells using a 3×3 nearest-neighbour mean (similar to CDO fillmiss). At 0.5° resolution some coastal cells carry no data; filling them before regridding avoids propagating gaps into the finer grid.
- 2b. Save filled 0.5° intermediate — the filled-but-not-regridded data is saved separately for use by the BCSD parameter fitting path in Notebook 02 (Option B), which fits distribution parameters at native CPC resolution and bilinearly interpolates them to the IMERG grid.
Regrid — linearly interpolate from 0.5° to 0.1° to match the IMERG grid. Because the input is already clipped, this step processes only a few hundred pixels instead of millions.
Mask clip — apply the AOI land mask to produce the final output with exactly the same spatial footprint as the IMERG extracts.

import os
import xarray as xr
import numpy as np
from scipy.ndimage import generic_filter

# Define directories
cpcuni_dir = os.path.join(main_dir, 'data/downloads/CPC_UNI')
clip_dir   = os.path.join(main_dir, 'data/downloads/CPC_UNI_extract/clip')
filled_05deg_dir = os.path.join(main_dir, 'data/downloads/CPC_UNI_extract/filled_05deg')
os.makedirs(clip_dir, exist_ok=True)
os.makedirs(filled_05deg_dir, exist_ok=True)

# Compression settings
comp = dict(zlib=True, complevel=5)

# Load AOI land mask once
mask_file = os.path.join(main_dir, 'data/mask/iso3/idn_subset.nc')
mask_ds = xr.open_dataset(mask_file)
land_mask = mask_ds['land'].astype(float)

# Buffered bounding box for coarse clip.
# Use 1.0° buffer (>= 2 CPC 0.5° cells) so that every IMERG pixel has
# sufficient CPC neighbours for bilinear interpolation during native-
# resolution parameter fitting (BCSD, Notebook 02).  This "all-touch"
# margin prevents NaN at the AOI boundary after interpolation.
buf = 1.0
lat_min = float(land_mask.lat.min()) - buf
lat_max = float(land_mask.lat.max()) + buf
lon_min = float(land_mask.lon.min()) - buf
lon_max = float(land_mask.lon.max()) + buf
print(f"Coarse clip bbox: lat [{lat_min:.2f}, {lat_max:.2f}], "
      f"lon [{lon_min:.2f}, {lon_max:.2f}]")

# Build a regional 0.1° target grid (not global)
target_lats = np.arange(
    np.floor(lat_min * 10) / 10 + 0.05,
    np.ceil(lat_max * 10) / 10,
    0.1
)
target_lons = np.arange(
    np.floor(lon_min * 10) / 10 + 0.05,
    np.ceil(lon_max * 10) / 10,
    0.1
)
print(f"Regional 0.1° grid: {len(target_lats)} lat × {len(target_lons)} lon")


def fillmiss_2d(arr, max_iter=50):
    """
    Iteratively fill NaN cells from their nearest non-NaN neighbours,
    spreading outward one cell per iteration (similar to CDO fillmiss).

    Each iteration replaces every NaN that has at least one non-NaN
    neighbour with the mean of its valid neighbours (3×3 window).
    Repeats until no NaN remain or max_iter is reached.
    """
    filled = arr.copy()
    for _ in range(max_iter):
        nans = np.isnan(filled)
        if not nans.any():
            break

        # Mean of valid neighbours (ignoring NaN) in a 3×3 window
        def _nanmean_or_nan(values):
            valid = values[~np.isnan(values)]
            return np.mean(valid) if valid.size > 0 else np.nan

        smoothed = generic_filter(filled, _nanmean_or_nan, size=3,
                                  mode='constant', cval=np.nan)
        filled[nans] = smoothed[nans]

    return filled


def process_year(year):
    """Coarse clip → spatial fill → save 0.5° → regrid → mask clip for one year."""
    input_path  = os.path.join(cpcuni_dir, f"precip.{year}.nc")
    output_path = os.path.join(clip_dir, f"idn_cpc_{year}.nc4")

    # Check whether BOTH outputs already exist
    filled_05_path = os.path.join(filled_05deg_dir, f"idn_cpc_filled05_{year}.nc4")
    both_exist = os.path.exists(output_path) and os.path.exists(filled_05_path)
    if both_exist:
        print(f"  {year}: both outputs already exist, skipping")
        return

    # --- 1) Coarse clip to buffered bbox ---
    ds = xr.open_dataset(input_path)
    if ds['lat'].values[0] > ds['lat'].values[-1]:
        ds = ds.sortby('lat')
    ds = ds.assign_coords(lon=((ds['lon'] + 180) % 360) - 180).sortby('lon')
    ds = ds.sel(lat=slice(lat_min, lat_max), lon=slice(lon_min, lon_max))

    # --- 2) Spatial fill (per timestep, 2D nearest-neighbour spread) ---
    #     Fills ocean/coastal NaN so every cell in the bbox has a value.
    precip = ds['precip'].values  # shape: (time, lat, lon)
    for t in range(precip.shape[0]):
        precip[t] = fillmiss_2d(precip[t])
    ds['precip'].values = precip

    # --- 2b) Save filled 0.5° intermediate for native-resolution fitting ---
    #     Used by Notebook 02 (Option B / BCSD) to fit CPC distribution
    #     parameters at the native 0.5° grid, then bilinearly interpolate
    #     to 0.1° — eliminating block boundary artefacts.
    if not os.path.exists(filled_05_path):
        ds_05 = ds.copy(deep=True)
        ds_05['precip'] = ds_05['precip'].transpose('time', 'lat', 'lon')
        ds_05['lat'].attrs.update({
            'long_name': 'Latitude', 'units': 'degrees_north',
            'axis': 'Y', 'standard_name': 'latitude'
        })
        ds_05['lon'].attrs.update({
            'long_name': 'Longitude', 'units': 'degrees_east',
            'axis': 'X', 'standard_name': 'longitude'
        })
        enc_05 = {var: comp for var in ds_05.data_vars}
        enc_05['lat'] = {'dtype': 'float32'}
        enc_05['lon'] = {'dtype': 'float32'}
        ds_05.attrs['Conventions'] = 'CF-1.8'
        ds_05.to_netcdf(filled_05_path, encoding=enc_05, engine='netcdf4')
        ds_05.close()
        print(f"  {year}: saved filled 0.5° intermediate")

    # --- 3) Regrid 0.5° → 0.1° (regional grid only) ---
    if not os.path.exists(output_path):
        ds = ds.interp(lat=target_lats, lon=target_lons, method='linear')
        ds['precip'] = ds['precip'].transpose('time', 'lat', 'lon')

        # --- 4) Mask clip to exact AOI ---
        mask_interp = land_mask.interp(lat=ds['lat'], lon=ds['lon'],
                                       method='nearest')
        ds = ds.where(mask_interp == 1, drop=True)

        # Coordinate attributes
        ds['lat'].attrs.update({
            'long_name': 'Latitude', 'units': 'degrees_north',
            'axis': 'Y', 'standard_name': 'latitude'
        })
        ds['lon'].attrs.update({
            'long_name': 'Longitude', 'units': 'degrees_east',
            'axis': 'X', 'standard_name': 'longitude'
        })

        encoding = {var: comp for var in ds.data_vars}
        encoding['lat'] = {'dtype': 'float32'}
        encoding['lon'] = {'dtype': 'float32'}
        ds.attrs['Conventions'] = 'CF-1.8'
        ds.to_netcdf(output_path, encoding=encoding, engine='netcdf4')
        print(f"  {year}: saved 0.1° regridded")

    ds.close()


# --- Process all years ---
for year in range(2001, 2026):
    print(f"Processing {year} ...")
    try:
        process_year(year)
    except Exception as e:
        print(f"  {year}: ERROR — {e}")

mask_ds.close()
print("CPC-UNI processing complete.")

Coarse clip bbox: lat [-12.05, 7.05], lon [93.95, 142.05]
Regional 0.1° grid: 192 lat × 482 lon
Processing 2001 ...
  2001: saved filled 0.5° intermediate
  2001: saved 0.1° regridded
Processing 2002 ...
  2002: saved filled 0.5° intermediate
  2002: saved 0.1° regridded
Processing 2003 ...
  2003: saved filled 0.5° intermediate
  2003: saved 0.1° regridded
Processing 2004 ...
  2004: saved filled 0.5° intermediate
  2004: saved 0.1° regridded
Processing 2005 ...
  2005: saved filled 0.5° intermediate
  2005: saved 0.1° regridded
Processing 2006 ...
  2006: saved filled 0.5° intermediate
  2006: saved 0.1° regridded
Processing 2007 ...
  2007: saved filled 0.5° intermediate
  2007: saved 0.1° regridded
Processing 2008 ...
  2008: saved filled 0.5° intermediate
  2008: saved 0.1° regridded
Processing 2009 ...
  2009: saved filled 0.5° intermediate
  2009: saved 0.1° regridded
Processing 2010 ...
  2010: saved filled 0.5° intermediate
  2010: saved 0.1° regridded
Processing 2011 ...
  2011: saved filled 0.5° intermediate
  2011: saved 0.1° regridded
Processing 2012 ...
  2012: saved filled 0.5° intermediate
  2012: saved 0.1° regridded
Processing 2013 ...
  2013: saved filled 0.5° intermediate
  2013: saved 0.1° regridded
Processing 2014 ...
  2014: saved filled 0.5° intermediate
  2014: saved 0.1° regridded
Processing 2015 ...
  2015: saved filled 0.5° intermediate
  2015: saved 0.1° regridded
Processing 2016 ...
  2016: saved filled 0.5° intermediate
  2016: saved 0.1° regridded
Processing 2017 ...
  2017: saved filled 0.5° intermediate
  2017: saved 0.1° regridded
Processing 2018 ...
  2018: saved filled 0.5° intermediate
  2018: saved 0.1° regridded
Processing 2019 ...
  2019: saved filled 0.5° intermediate
  2019: saved 0.1° regridded
Processing 2020 ...
  2020: saved filled 0.5° intermediate
  2020: saved 0.1° regridded
Processing 2021 ...
  2021: saved filled 0.5° intermediate
  2021: saved 0.1° regridded
Processing 2022 ...
  2022: saved filled 0.5° intermediate
  2022: saved 0.1° regridded
Processing 2023 ...
  2023: saved filled 0.5° intermediate
  2023: saved 0.1° regridded
Processing 2024 ...
  2024: saved filled 0.5° intermediate
  2024: saved 0.1° regridded
Processing 2025 ...
  2025: saved filled 0.5° intermediate
  2025: saved 0.1° regridded
CPC-UNI processing complete.

Step 3: Combine into a Single Multi-Year NetCDF (0.1°)

The bias-correction pipeline (notebook 02) expects one CPC-UNI file covering the full daily time series. This step concatenates all per-year clipped files along the time dimension and writes the result to data/input/idn_cli_cpc_daily.nc4.

import xarray as xr
import netCDF4
import os
import glob
import numpy as np

clip_dir = os.path.join(main_dir, 'data/downloads/CPC_UNI_extract/clip')
output_file_cpc = os.path.join(main_dir, 'data/input/cpcuni/idn_cpcuni.nc4')

os.makedirs(os.path.dirname(output_file_cpc), exist_ok=True)

file_list = sorted(glob.glob(os.path.join(clip_dir, "idn_cpc_*.nc4")))
print(f"Found {len(file_list)} yearly CPC-UNI clipped files")

if not file_list:
    print("No clipped files found. Run Step 2 first.")
else:
    nc_out = None
    total_days = 0

    for i, fpath in enumerate(file_list):
        year = os.path.basename(fpath).replace('idn_cpc_', '').replace('.nc4', '')
        ds = xr.open_dataset(fpath).sortby('time')
        ndays = len(ds.time)

        if i == 0:
            # Create the output file with unlimited time dimension
            nc_out = netCDF4.Dataset(output_file_cpc, 'w', format='NETCDF4')
            nc_out.Conventions = 'CF-1.8'
            nc_out.title = 'CPC Unified daily precipitation (0.1 deg, AOI clipped)'

            # Dimensions
            nc_out.createDimension('time', None)  # unlimited
            nc_out.createDimension('lat', len(ds.lat))
            nc_out.createDimension('lon', len(ds.lon))

            # Time variable
            t_var = nc_out.createVariable('time', 'f8', ('time',))
            t_var.units = 'days since 1900-01-01'
            t_var.calendar = 'standard'
            t_var.standard_name = 'time'

            # Lat/lon variables
            lat_var = nc_out.createVariable('lat', 'f4', ('lat',))
            lat_var[:] = ds.lat.values.astype('float32')
            lat_var.units = 'degrees_north'
            lat_var.standard_name = 'latitude'
            lat_var.axis = 'Y'

            lon_var = nc_out.createVariable('lon', 'f4', ('lon',))
            lon_var[:] = ds.lon.values.astype('float32')
            lon_var.units = 'degrees_east'
            lon_var.standard_name = 'longitude'
            lon_var.axis = 'X'

            # Precip variable with compression
            p_var = nc_out.createVariable(
                'precip', 'f4', ('time', 'lat', 'lon'),
                zlib=True, complevel=5,
                chunksizes=(1, len(ds.lat), len(ds.lon)),
                fill_value=np.float32(9.96921e+36)
            )
            if 'units' in ds['precip'].attrs:
                p_var.units = ds['precip'].attrs['units']
            if 'long_name' in ds['precip'].attrs:
                p_var.long_name = ds['precip'].attrs['long_name']

        # Convert time to numeric values
        time_vals = netCDF4.date2num(
            ds.time.values.astype('datetime64[ms]').astype('O'),
            units='days since 1900-01-01',
            calendar='standard'
        )

        # Append this year's data
        nc_out.variables['time'][total_days:total_days + ndays] = time_vals
        nc_out.variables['precip'][total_days:total_days + ndays, :, :] = \
            ds['precip'].values
        nc_out.sync()

        total_days += ndays
        print(f"  [{i+1}/{len(file_list)}] {year}: +{ndays} days "
              f"(total: {total_days})")
        ds.close()

    nc_out.close()

    # Final verification
    with xr.open_dataset(output_file_cpc) as final:
        print(f"\nDone: {output_file_cpc}")
        print(f"  Time: {str(final.time.values[0])[:10]} to "
              f"{str(final.time.values[-1])[:10]}")
        print(f"  Shape: {dict(final.sizes)}")

Found 25 yearly CPC-UNI clipped files
  [1/25] 2001: +365 days (total: 365)
  [2/25] 2002: +365 days (total: 730)
  [3/25] 2003: +365 days (total: 1095)
  [4/25] 2004: +366 days (total: 1461)
  [5/25] 2005: +365 days (total: 1826)
  [6/25] 2006: +365 days (total: 2191)
  [7/25] 2007: +365 days (total: 2556)
  [8/25] 2008: +366 days (total: 2922)
  [9/25] 2009: +365 days (total: 3287)
  [10/25] 2010: +365 days (total: 3652)
  [11/25] 2011: +365 days (total: 4017)
  [12/25] 2012: +366 days (total: 4383)
  [13/25] 2013: +365 days (total: 4748)
  [14/25] 2014: +365 days (total: 5113)
  [15/25] 2015: +365 days (total: 5478)
  [16/25] 2016: +366 days (total: 5844)
  [17/25] 2017: +365 days (total: 6209)
  [18/25] 2018: +365 days (total: 6574)
  [19/25] 2019: +365 days (total: 6939)
  [20/25] 2020: +366 days (total: 7305)
  [21/25] 2021: +365 days (total: 7670)
  [22/25] 2022: +365 days (total: 8035)
  [23/25] 2023: +365 days (total: 8400)
  [24/25] 2024: +366 days (total: 8766)
  [25/25] 2025: +365 days (total: 9131)

Done: /content/drive/MyDrive/hybrid-bias-correction/data/input/cpcuni/idn_cpcuni.nc4
  Time: 2001-01-01 to 2025-12-31
  Shape: {'time': 9131, 'lat': 171, 'lon': 461}

Step 3b: Combine Filled 0.5° Files into a Single Multi-Year NetCDF

The BCSD parameter fitting path in Notebook 02 (Option B) needs the CPC data at its native 0.5° resolution — gap-filled but not regridded. This step concatenates the per-year filled 0.5° files produced in Step 2b into a single multi-year file at data/input/cpcuni/idn_cpcuni_native05.nc4.

This file is referenced by cpc_native_file in config.yml.

import xarray as xr
import netCDF4
import os
import glob
import numpy as np

filled_05deg_dir = os.path.join(main_dir, 'data/downloads/CPC_UNI_extract/filled_05deg')
output_file_cpc_05 = os.path.join(main_dir, 'data/input/cpcuni/idn_cpcuni_native05.nc4')

os.makedirs(os.path.dirname(output_file_cpc_05), exist_ok=True)

file_list = sorted(glob.glob(os.path.join(filled_05deg_dir, "idn_cpc_filled05_*.nc4")))
print(f"Found {len(file_list)} yearly CPC-UNI filled 0.5° files")

if not file_list:
    print("No filled 0.5° files found. Run Step 2 first.")
else:
    nc_out = None
    total_days = 0

    for i, fpath in enumerate(file_list):
        year = os.path.basename(fpath).replace('idn_cpc_filled05_', '').replace('.nc4', '')
        ds = xr.open_dataset(fpath).sortby('time')
        ndays = len(ds.time)

        if i == 0:
            nc_out = netCDF4.Dataset(output_file_cpc_05, 'w', format='NETCDF4')
            nc_out.Conventions = 'CF-1.8'
            nc_out.title = 'CPC Unified daily precipitation (native 0.5 deg, gap-filled, AOI buffered)'

            nc_out.createDimension('time', None)
            nc_out.createDimension('lat', len(ds.lat))
            nc_out.createDimension('lon', len(ds.lon))

            t_var = nc_out.createVariable('time', 'f8', ('time',))
            t_var.units = 'days since 1900-01-01'
            t_var.calendar = 'standard'
            t_var.standard_name = 'time'

            lat_var = nc_out.createVariable('lat', 'f4', ('lat',))
            lat_var[:] = ds.lat.values.astype('float32')
            lat_var.units = 'degrees_north'
            lat_var.standard_name = 'latitude'
            lat_var.axis = 'Y'

            lon_var = nc_out.createVariable('lon', 'f4', ('lon',))
            lon_var[:] = ds.lon.values.astype('float32')
            lon_var.units = 'degrees_east'
            lon_var.standard_name = 'longitude'
            lon_var.axis = 'X'

            p_var = nc_out.createVariable(
                'precip', 'f4', ('time', 'lat', 'lon'),
                zlib=True, complevel=5,
                chunksizes=(1, len(ds.lat), len(ds.lon)),
                fill_value=np.float32(9.96921e+36)
            )
            if 'units' in ds['precip'].attrs:
                p_var.units = ds['precip'].attrs['units']
            if 'long_name' in ds['precip'].attrs:
                p_var.long_name = ds['precip'].attrs['long_name']

        time_vals = netCDF4.date2num(
            ds.time.values.astype('datetime64[ms]').astype('O'),
            units='days since 1900-01-01',
            calendar='standard'
        )

        nc_out.variables['time'][total_days:total_days + ndays] = time_vals
        nc_out.variables['precip'][total_days:total_days + ndays, :, :] = \
            ds['precip'].values
        nc_out.sync()

        total_days += ndays
        print(f"  [{i+1}/{len(file_list)}] {year}: +{ndays} days "
              f"(total: {total_days})")
        ds.close()

    nc_out.close()

    with xr.open_dataset(output_file_cpc_05) as final:
        print(f"\nDone: {output_file_cpc_05}")
        print(f"  Time: {str(final.time.values[0])[:10]} to "
              f"{str(final.time.values[-1])[:10]}")
        print(f"  Shape: {dict(final.sizes)}")
        res = abs(float(final.lat[1] - final.lat[0]))
        print(f"  Resolution: ~{res:.2f}°")

Found 25 yearly CPC-UNI filled 0.5° files
  [1/25] 2001: +365 days (total: 365)
  [2/25] 2002: +365 days (total: 730)
  [3/25] 2003: +365 days (total: 1095)
  [4/25] 2004: +366 days (total: 1461)
  [5/25] 2005: +365 days (total: 1826)
  [6/25] 2006: +365 days (total: 2191)
  [7/25] 2007: +365 days (total: 2556)
  [8/25] 2008: +366 days (total: 2922)
  [9/25] 2009: +365 days (total: 3287)
  [10/25] 2010: +365 days (total: 3652)
  [11/25] 2011: +365 days (total: 4017)
  [12/25] 2012: +366 days (total: 4383)
  [13/25] 2013: +365 days (total: 4748)
  [14/25] 2014: +365 days (total: 5113)
  [15/25] 2015: +365 days (total: 5478)
  [16/25] 2016: +366 days (total: 5844)
  [17/25] 2017: +365 days (total: 6209)
  [18/25] 2018: +365 days (total: 6574)
  [19/25] 2019: +365 days (total: 6939)
  [20/25] 2020: +366 days (total: 7305)
  [21/25] 2021: +365 days (total: 7670)
  [22/25] 2022: +365 days (total: 8035)
  [23/25] 2023: +365 days (total: 8400)
  [24/25] 2024: +366 days (total: 8766)
  [25/25] 2025: +365 days (total: 9131)

Done: /content/drive/MyDrive/hybrid-bias-correction/data/input/cpcuni/idn_cpcuni_native05.nc4
  Time: 2001-01-01 to 2025-12-31
  Shape: {'time': 9131, 'lat': 38, 'lon': 96}
  Resolution: ~0.50°

6 Summary

This notebook acquired and prepared the precipitation datasets that drive the remainder of the hybrid bias-correction workflow:

Dataset	Source	Native resolution	Final file
IMERG Late Run (v07)	NASA GES DISC THREDDS	0.1° half-hourly → daily	`data/input/idn_cli_imergl_daily.nc4`
IMERG Final Run (v07)	NASA GES DISC THREDDS	0.1° half-hourly → daily	`data/input/idn_cli_imergf_daily.nc4`
CPC Unified Gauge (0.1°)	NOAA PSL	0.5° daily → regridded 0.1°	`data/input/idn_cli_cpc_daily.nc4`
CPC Unified Gauge (native 0.5°)	NOAA PSL	0.5° daily, gap-filled	`data/input/cpcuni/idn_cpcuni_native05.nc4`

For each dataset the pipeline followed the same three steps: (1) download raw yearly subsets, (2) extract the precipitation variable and clip to the AOI land mask produced in Notebook 00, and (3) concatenate all years into a single multi-year NetCDF file with CF-1.8 conventions and zlib compression.

For CPC-UNI, an additional intermediate output is saved at Step 2b: the gap-filled 0.5° data before regridding. This native-resolution file is used by the BCSD parameter fitting path in Notebook 02 to fit distribution parameters at the reference resolution and bilinearly interpolate to 0.1°, eliminating the 5×5 block boundary artefact that arises from nearest-neighbour regridding.

Next step → Open Notebook 02 to run the LS → LSEQM → LSEQM+DL bias-correction chain on the IMERG Late Run using CPC-UNI as the gauge reference.