In [ ]:

# Uncomment to install climdata in Google Colab or other environments
# !pip install climdata

# Uncomment to install climdata in Google Colab or other environments # !pip install climdata

Bias Correction & Statistical Downscaling (BCSD)¶

This notebook walks through the climdata.sdba BCSD pipeline: fine-resolution MSWX observations are used to bias-correct and downscale coarse CMIP6 output over Europe.

Steps: regrid observations to the GCM grid → bias-correct at coarse resolution → statistically downscale back to the fine grid.

In [ ]:

import xarray as xr
import numpy as np
# Import climdata BCSD modules
from climdata.sdba import BCSD, BiasCorrection, StatisticalDownscaling, regrid_to_coarse

print("✓ All imports successful!")

import xarray as xr import numpy as np # Import climdata BCSD modules from climdata.sdba import BCSD, BiasCorrection, StatisticalDownscaling, regrid_to_coarse print("✓ All imports successful!")

In [ ]:

# Time periods
hist_start = '2004-01-01'
hist_end = '2014-12-31'
fut_start = '2015-01-01'
fut_end = '2050-12-31'

print(f"Historical period: {hist_start} to {hist_end}")
print(f"Future period: {fut_start} to {fut_end}")

# Time periods hist_start = '2004-01-01' hist_end = '2014-12-31' fut_start = '2015-01-01' fut_end = '2050-12-31' print(f"Historical period: {hist_start} to {hist_end}") print(f"Future period: {fut_start} to {fut_end}")

Load Fine-Resolution Reference Data (MSWX)¶

MSWX provides global weather data at 0.1° resolution, blending multiple data sources.

In [ ]:

from climdata import ClimData

options = [
    "dataset=mswx",  # Select the MSWX dataset for extraction
    "region=europe", # Select the region
    "variables=[tas]",
    f"time_range.start_date={hist_start}",  # Start date of extraction
    f"time_range.end_date={hist_end}",    # End date of extraction
    "data_dir=./data",    # Local directory to store downloaded/intermediate files
]
mswx = ClimData(overrides=options)

obs_fine_real = mswx.extract()

print("\n✓ MSWX data loaded!")
print(f"  Dimensions: {obs_fine_real.dims}")
print(f"  Variables: {list(obs_fine_real.data_vars)}")
print("  Resolution: ~0.1° (~10 km)")
print(f"  Time range: {obs_fine_real.time.values[0]} to {obs_fine_real.time.values[-1]}")

from climdata import ClimData options = [ "dataset=mswx", # Select the MSWX dataset for extraction "region=europe", # Select the region "variables=[tas]", f"time_range.start_date={hist_start}", # Start date of extraction f"time_range.end_date={hist_end}", # End date of extraction "data_dir=./data", # Local directory to store downloaded/intermediate files ] mswx = ClimData(overrides=options) obs_fine_real = mswx.extract() print("\n✓ MSWX data loaded!") print(f" Dimensions: {obs_fine_real.dims}") print(f" Variables: {list(obs_fine_real.data_vars)}") print(" Resolution: ~0.1° (~10 km)") print(f" Time range: {obs_fine_real.time.values[0]} to {obs_fine_real.time.values[-1]}")

Load Coarse GCM Data (CMIP6)¶

Load historical and future CMIP6 data for the same region.

In [ ]:

# Choose a CMIP6 model (using MPI-ESM1-2-HR as example)
cmip_model = 'MPI-ESM1-2-HR'

options = [
    "dataset=cmip",  # Select the CMIP6 dataset for extraction
    "region=europe", # Select the region
    "variables=[tas]",
    f"time_range.start_date={hist_start}",  # Start date of extraction
    f"time_range.end_date={hist_end}",    # End date of extraction
    "data_dir=./data",    # Local directory to store downloaded/intermediate files
    f"source_id={cmip_model}",
    "experiment_id=historical",
]
cmip_hist = ClimData(overrides=options)

sim_hist_coarse_real = cmip_hist.extract()

print("\n✓ Historical CMIP6 data loaded!")
print(f"  Model: {cmip_model}")
print(f"  Dimensions: {sim_hist_coarse_real.dims}")
print("  Resolution: ~1-2° (model-dependent)")
print(f"  Time range: {sim_hist_coarse_real.time.values[0]} to {sim_hist_coarse_real.time.values[-1]}")

# Choose a CMIP6 model (using MPI-ESM1-2-HR as example) cmip_model = 'MPI-ESM1-2-HR' options = [ "dataset=cmip", # Select the CMIP6 dataset for extraction "region=europe", # Select the region "variables=[tas]", f"time_range.start_date={hist_start}", # Start date of extraction f"time_range.end_date={hist_end}", # End date of extraction "data_dir=./data", # Local directory to store downloaded/intermediate files f"source_id={cmip_model}", "experiment_id=historical", ] cmip_hist = ClimData(overrides=options) sim_hist_coarse_real = cmip_hist.extract() print("\n✓ Historical CMIP6 data loaded!") print(f" Model: {cmip_model}") print(f" Dimensions: {sim_hist_coarse_real.dims}") print(" Resolution: ~1-2° (model-dependent)") print(f" Time range: {sim_hist_coarse_real.time.values[0]} to {sim_hist_coarse_real.time.values[-1]}")

In [ ]:

# Same model, future scenario
options = [
    "dataset=cmip",  # Select the CMIP6 dataset for extraction
    "region=europe", # Select the region
    "variables=[tas]",
    f"time_range.start_date={fut_start}",  # Start date of extraction
    f"time_range.end_date={fut_end}",    # End date of extraction
    "data_dir=./data",    # Local directory to store downloaded/intermediate files
    f"source_id={cmip_model}",
    "experiment_id=ssp585",
]
cmip_fut = ClimData(overrides=options)

sim_fut_coarse_real = cmip_fut.extract()

print("\n✓ Future CMIP6 data loaded!")
print("  Scenario: SSP5-8.5")
print(f"  Dimensions: {sim_fut_coarse_real.dims}")
print(f"  Time range: {sim_fut_coarse_real.time.values[0]} to {sim_fut_coarse_real.time.values[-1]}")

# Same model, future scenario options = [ "dataset=cmip", # Select the CMIP6 dataset for extraction "region=europe", # Select the region "variables=[tas]", f"time_range.start_date={fut_start}", # Start date of extraction f"time_range.end_date={fut_end}", # End date of extraction "data_dir=./data", # Local directory to store downloaded/intermediate files f"source_id={cmip_model}", "experiment_id=ssp585", ] cmip_fut = ClimData(overrides=options) sim_fut_coarse_real = cmip_fut.extract() print("\n✓ Future CMIP6 data loaded!") print(" Scenario: SSP5-8.5") print(f" Dimensions: {sim_fut_coarse_real.dims}") print(f" Time range: {sim_fut_coarse_real.time.values[0]} to {sim_fut_coarse_real.time.values[-1]}")

Inspect the Data¶

Check shapes, resolutions and whether the grids are compatible with ISIMIP3BASD downscaling (which needs integer downscaling factors).

In [ ]:

# Check data summary
print("Data Summary:")
print("\nObservations (MSWX):")
print("  Variable: tas (temperature)")
print(f"  Shape: {obs_fine_real['tas'].shape}")
print(f"  Resolution: ~{(obs_fine_real.lat.values[1] - obs_fine_real.lat.values[0]):.2f}°")

print("\nHistorical GCM (CMIP6):")
print("  Variable: tas")
print(f"  Shape: {sim_hist_coarse_real['tas'].shape}")
print(f"  Resolution: ~{(sim_hist_coarse_real.lat.values[1] - sim_hist_coarse_real.lat.values[0]):.2f}°")

print("\nFuture GCM (CMIP6):")
print("  Variable: tas")
print(f"  Shape: {sim_fut_coarse_real['tas'].shape}")

# Check if grids are compatible for downscaling
lat_ratio = len(obs_fine_real.lat) / len(sim_hist_coarse_real.lat)
lon_ratio = len(obs_fine_real.lon) / len(sim_hist_coarse_real.lon)

print("\n📊 Downscaling factors:")
print(f"  Latitude: {lat_ratio:.2f}x")
print(f"  Longitude: {lon_ratio:.2f}x")

if lat_ratio == int(lat_ratio) and lon_ratio == int(lon_ratio):
    print("  ✓ Grids are compatible for ISIMIP3BASD downscaling!")
else:
    print("  ⚠ Warning: Grids are not exact integer multiples.")
    print("    ISIMIP3BASD downscaling requires integer downscaling factors.")
    print("    Consider regridding to compatible resolutions first.")

# Check data summary print("Data Summary:") print("\nObservations (MSWX):") print(" Variable: tas (temperature)") print(f" Shape: {obs_fine_real['tas'].shape}") print(f" Resolution: ~{(obs_fine_real.lat.values[1] - obs_fine_real.lat.values[0]):.2f}°") print("\nHistorical GCM (CMIP6):") print(" Variable: tas") print(f" Shape: {sim_hist_coarse_real['tas'].shape}") print(f" Resolution: ~{(sim_hist_coarse_real.lat.values[1] - sim_hist_coarse_real.lat.values[0]):.2f}°") print("\nFuture GCM (CMIP6):") print(" Variable: tas") print(f" Shape: {sim_fut_coarse_real['tas'].shape}") # Check if grids are compatible for downscaling lat_ratio = len(obs_fine_real.lat) / len(sim_hist_coarse_real.lat) lon_ratio = len(obs_fine_real.lon) / len(sim_hist_coarse_real.lon) print("\n📊 Downscaling factors:") print(f" Latitude: {lat_ratio:.2f}x") print(f" Longitude: {lon_ratio:.2f}x") if lat_ratio == int(lat_ratio) and lon_ratio == int(lon_ratio): print(" ✓ Grids are compatible for ISIMIP3BASD downscaling!") else: print(" ⚠ Warning: Grids are not exact integer multiples.") print(" ISIMIP3BASD downscaling requires integer downscaling factors.") print(" Consider regridding to compatible resolutions first.")

In [ ]:

# Initialize BCSD for temperature
bcsd = BCSD(
    variable='tas',
    regridding_tool='xesmf',  # or 'cdo' if xESMF not available
    regridding_method='conservative',  # area-weighted conservative regridding
    bias_correction_kwargs={
        'n_processes': 4  # Use parallel processing
    },
    downscaling_kwargs={
        'n_iterations': 20,  # MBCn iterations
        'n_processes': 4
    }
)

print("BCSD pipeline configured!")

# Initialize BCSD for temperature bcsd = BCSD( variable='tas', regridding_tool='xesmf', # or 'cdo' if xESMF not available regridding_method='conservative', # area-weighted conservative regridding bias_correction_kwargs={ 'n_processes': 4 # Use parallel processing }, downscaling_kwargs={ 'n_iterations': 20, # MBCn iterations 'n_processes': 4 } ) print("BCSD pipeline configured!")

Run the Workflow Step by Step¶

The pipeline is:

1. Regrid the fine observations onto the coarse GCM grid
2. Bias-correct at coarse resolution
3. Downscale the corrected data back to fine resolution

In [ ]:

# Step 1: Regrid fine observations to coarse GCM grid
print("Step 1: Regridding observations to coarse resolution...")
obs_hist_coarse = regrid_to_coarse(
    obs_fine_real,
    sim_hist_coarse_real,
    method='conservative',
    regridding_tool='xesmf'
)
print(f"  ✓ Coarse observations shape: {obs_hist_coarse['tas'].shape}")

# Ensure all datasets have exactly the same spatial grid
print("\nStep 1b: Aligning spatial grids...")

# Remove extra dimensions from CMIP data (source_id dimension)
if 'source_id' in sim_hist_coarse_real.dims:
    sim_hist_coarse_real = sim_hist_coarse_real.squeeze('source_id', drop=True)
    print("  ✓ Removed source_id dimension from historical data")
if 'source_id' in sim_fut_coarse_real.dims:
    sim_fut_coarse_real = sim_fut_coarse_real.squeeze('source_id', drop=True)
    print("  ✓ Removed source_id dimension from future data")

# Use sim_hist as the reference grid - ensure obs and sim_fut match it exactly
obs_hist_coarse = obs_hist_coarse.sel(
    lat=sim_hist_coarse_real.lat,
    lon=sim_hist_coarse_real.lon,
    method='nearest'
)
# Explicitly assign coordinates to ensure perfect match
obs_hist_coarse = obs_hist_coarse.assign_coords({
    'lat': sim_hist_coarse_real.lat,
    'lon': sim_hist_coarse_real.lon
})

sim_fut_coarse_aligned = sim_fut_coarse_real.sel(
    lat=sim_hist_coarse_real.lat,
    lon=sim_hist_coarse_real.lon,
    method='nearest'
)
# Explicitly assign coordinates
sim_fut_coarse_aligned = sim_fut_coarse_aligned.assign_coords({
    'lat': sim_hist_coarse_real.lat,
    'lon': sim_hist_coarse_real.lon
})

print(f"  ✓ Aligned shapes - obs: {obs_hist_coarse['tas'].shape}, hist: {sim_hist_coarse_real['tas'].shape}, fut: {sim_fut_coarse_aligned['tas'].shape}")
print(f"  ✓ Coordinate match - lat: {(obs_hist_coarse.lat == sim_hist_coarse_real.lat).all().values}, lon: {(obs_hist_coarse.lon == sim_hist_coarse_real.lon).all().values}")

# Step 2: Bias correction at coarse resolution
print("\nStep 2: Bias correction...")
bc = BiasCorrection(
    variable='tas'
)
sim_fut_ba = bc.correct(
    obs_hist=obs_hist_coarse,
    sim_hist=sim_hist_coarse_real,
    sim_fut=sim_fut_coarse_aligned
)
print(f"  ✓ Bias-corrected shape: {sim_fut_ba['tas'].shape}")

# Step 1: Regrid fine observations to coarse GCM grid print("Step 1: Regridding observations to coarse resolution...") obs_hist_coarse = regrid_to_coarse( obs_fine_real, sim_hist_coarse_real, method='conservative', regridding_tool='xesmf' ) print(f" ✓ Coarse observations shape: {obs_hist_coarse['tas'].shape}") # Ensure all datasets have exactly the same spatial grid print("\nStep 1b: Aligning spatial grids...") # Remove extra dimensions from CMIP data (source_id dimension) if 'source_id' in sim_hist_coarse_real.dims: sim_hist_coarse_real = sim_hist_coarse_real.squeeze('source_id', drop=True) print(" ✓ Removed source_id dimension from historical data") if 'source_id' in sim_fut_coarse_real.dims: sim_fut_coarse_real = sim_fut_coarse_real.squeeze('source_id', drop=True) print(" ✓ Removed source_id dimension from future data") # Use sim_hist as the reference grid - ensure obs and sim_fut match it exactly obs_hist_coarse = obs_hist_coarse.sel( lat=sim_hist_coarse_real.lat, lon=sim_hist_coarse_real.lon, method='nearest' ) # Explicitly assign coordinates to ensure perfect match obs_hist_coarse = obs_hist_coarse.assign_coords({ 'lat': sim_hist_coarse_real.lat, 'lon': sim_hist_coarse_real.lon }) sim_fut_coarse_aligned = sim_fut_coarse_real.sel( lat=sim_hist_coarse_real.lat, lon=sim_hist_coarse_real.lon, method='nearest' ) # Explicitly assign coordinates sim_fut_coarse_aligned = sim_fut_coarse_aligned.assign_coords({ 'lat': sim_hist_coarse_real.lat, 'lon': sim_hist_coarse_real.lon }) print(f" ✓ Aligned shapes - obs: {obs_hist_coarse['tas'].shape}, hist: {sim_hist_coarse_real['tas'].shape}, fut: {sim_fut_coarse_aligned['tas'].shape}") print(f" ✓ Coordinate match - lat: {(obs_hist_coarse.lat == sim_hist_coarse_real.lat).all().values}, lon: {(obs_hist_coarse.lon == sim_hist_coarse_real.lon).all().values}") # Step 2: Bias correction at coarse resolution print("\nStep 2: Bias correction...") bc = BiasCorrection( variable='tas' ) sim_fut_ba = bc.correct( obs_hist=obs_hist_coarse, sim_hist=sim_hist_coarse_real, sim_fut=sim_fut_coarse_aligned ) print(f" ✓ Bias-corrected shape: {sim_fut_ba['tas'].shape}")

In [ ]:

# Step 3: Statistical downscaling
print("\nStep 3: Statistical downscaling...")

# Check if grids are compatible for downscaling
lat_factor = len(obs_fine_real.lat) / len(sim_fut_ba.lat)
lon_factor = len(obs_fine_real.lon) / len(sim_fut_ba.lon)

print(f"  Current downscaling factors: lat={lat_factor:.2f}x, lon={lon_factor:.2f}x")

if lat_factor != int(lat_factor) or lon_factor != int(lon_factor):
    print("  ⚠ Grids not compatible - creating compatible fine grid...")
    
    # Create a compatible fine grid (integer multiples)
    target_lat_factor = round(lat_factor)
    target_lon_factor = round(lon_factor)
    
    # Create target coordinates
    coarse_lat = sim_fut_ba.lat.values
    coarse_lon = sim_fut_ba.lon.values
    
    # Generate fine grid with exact integer spacing
    fine_lat = np.linspace(coarse_lat[0], coarse_lat[-1], len(coarse_lat) * target_lat_factor)
    fine_lon = np.linspace(coarse_lon[0], coarse_lon[-1], len(coarse_lon) * target_lon_factor)
    
    # Create a dummy dataset with the target fine grid
    target_fine_grid = xr.Dataset({
        'dummy': (('lat', 'lon'), np.zeros((len(fine_lat), len(fine_lon))))
    }, coords={'lat': fine_lat, 'lon': fine_lon})
    
    # Regrid obs_fine to the compatible grid
    print(f"  Regridding to {len(fine_lat)}×{len(fine_lon)} ({target_lat_factor}×{target_lon_factor} factors)...")
    
    import xesmf as xe
    regridder = xe.Regridder(obs_fine_real, target_fine_grid, 'bilinear')
    obs_fine_compatible = regridder(obs_fine_real)
    # regridder.clean_weight_file()
    
    # Ensure coordinates have proper metadata for iris
    obs_fine_compatible['lat'].attrs.update({
        'standard_name': 'latitude',
        'long_name': 'latitude',
        'units': 'degrees_north',
        'axis': 'Y'
    })
    obs_fine_compatible['lon'].attrs.update({
        'standard_name': 'longitude',
        'long_name': 'longitude',
        'units': 'degrees_east',
        'axis': 'X'
    })
    obs_fine_compatible['time'].attrs.update({
        'standard_name': 'time',
        'long_name': 'time',
        'axis': 'T'
    })
    
    print(f"  ✓ Compatible fine grid created: {obs_fine_compatible['tas'].shape}")
else:
    obs_fine_compatible = obs_fine_real
    print("  ✓ Grids already compatible!")

sd = StatisticalDownscaling(
    variable='tas',
    n_iterations=20
)
result = sd.downscale(
    obs_fine=obs_fine_compatible,
    sim_coarse=sim_fut_ba
)
print(f"  ✓ Downscaled shape: {result['tas'].shape}")

print("\n✅ BCSD workflow completed successfully!")
print(f"Output dimensions: {result.dims}")
print(f"Output resolution: {(result.lat.values[1] - result.lat.values[0]):.3f}°")

# Step 3: Statistical downscaling print("\nStep 3: Statistical downscaling...") # Check if grids are compatible for downscaling lat_factor = len(obs_fine_real.lat) / len(sim_fut_ba.lat) lon_factor = len(obs_fine_real.lon) / len(sim_fut_ba.lon) print(f" Current downscaling factors: lat={lat_factor:.2f}x, lon={lon_factor:.2f}x") if lat_factor != int(lat_factor) or lon_factor != int(lon_factor): print(" ⚠ Grids not compatible - creating compatible fine grid...") # Create a compatible fine grid (integer multiples) target_lat_factor = round(lat_factor) target_lon_factor = round(lon_factor) # Create target coordinates coarse_lat = sim_fut_ba.lat.values coarse_lon = sim_fut_ba.lon.values # Generate fine grid with exact integer spacing fine_lat = np.linspace(coarse_lat[0], coarse_lat[-1], len(coarse_lat) * target_lat_factor) fine_lon = np.linspace(coarse_lon[0], coarse_lon[-1], len(coarse_lon) * target_lon_factor) # Create a dummy dataset with the target fine grid target_fine_grid = xr.Dataset({ 'dummy': (('lat', 'lon'), np.zeros((len(fine_lat), len(fine_lon)))) }, coords={'lat': fine_lat, 'lon': fine_lon}) # Regrid obs_fine to the compatible grid print(f" Regridding to {len(fine_lat)}×{len(fine_lon)} ({target_lat_factor}×{target_lon_factor} factors)...") import xesmf as xe regridder = xe.Regridder(obs_fine_real, target_fine_grid, 'bilinear') obs_fine_compatible = regridder(obs_fine_real) # regridder.clean_weight_file() # Ensure coordinates have proper metadata for iris obs_fine_compatible['lat'].attrs.update({ 'standard_name': 'latitude', 'long_name': 'latitude', 'units': 'degrees_north', 'axis': 'Y' }) obs_fine_compatible['lon'].attrs.update({ 'standard_name': 'longitude', 'long_name': 'longitude', 'units': 'degrees_east', 'axis': 'X' }) obs_fine_compatible['time'].attrs.update({ 'standard_name': 'time', 'long_name': 'time', 'axis': 'T' }) print(f" ✓ Compatible fine grid created: {obs_fine_compatible['tas'].shape}") else: obs_fine_compatible = obs_fine_real print(" ✓ Grids already compatible!") sd = StatisticalDownscaling( variable='tas', n_iterations=20 ) result = sd.downscale( obs_fine=obs_fine_compatible, sim_coarse=sim_fut_ba ) print(f" ✓ Downscaled shape: {result['tas'].shape}") print("\n✅ BCSD workflow completed successfully!") print(f"Output dimensions: {result.dims}") print(f"Output resolution: {(result.lat.values[1] - result.lat.values[0]):.3f}°")