第22章 OGGM教程全集

初级 中级 高级

本章提供了一套完整的教程，涵盖OGGM的每个主要方面，从安装到高级专题（如动态预热（dynamic spinup）和入海型冰川建模）。 教程中的代码分为三类：可以从空目录直接运行的完整脚本、需要沿用前面已经生成的gdir/gdirs对象的连续练习、以及用于解释思想的示意代码。每段关键代码前会说明依赖条件，读者不必把所有片段都视为可单独复制运行的脚本。

教程一：入门指南

1.1 安装OGGM

OGGM可通过conda（推荐）和pip安装：

# Option A: conda (recommended — handles all C dependencies)
conda create -n oggm_env python=3.11
conda activate oggm_env
conda install -c conda-forge oggm-deps
pip install oggm

# Option B: pip-only (requires system GDAL and netCDF4 libs)
pip install oggm

# Option C: development install (for modifying OGGM source)
git clone https://github.com/OGGM/oggm.git
cd oggm
pip install -e .

1.2 基础设置

import oggm
from oggm import cfg, workflow, utils

# Step 1 — Initialize OGGM (load configuration)
cfg.initialize()

# Step 2 — Set the working directory (where glacier data will be stored)
cfg.PATHS['working_dir'] = '/path/to/your/OGGM_data'

# Step 3 — Enable multiprocessing (optional but recommended)
cfg.PARAMS['use_multiprocessing'] = True

1.3 你的第一个冰川：Hintereisferner

# Define a single glacier to work with
rgi_ids = ['RGI60-11.00897']  # Hintereisferner, Austrian Alps

# Create GlacierDirectory object(s)
gdirs = workflow.init_glacier_directories(
    rgi_ids,
    from_prepro_level=0,  # start from scratch
    prepro_border=10       # map border in pixels
)

print(f"Created {len(gdirs)} glacier directory")
print(f"Glacier: {gdirs[0].rgi_id}")
print(f"Area: {gdirs[0].rgi_area_km2:.2f} km²")
print(f"Data directory: {gdirs[0].dir}")

预期输出：

Created 1 glacier directory
Glacier: RGI60-11.00897
Area: 8.45 km²
Data directory: /path/to/OGGM_data/per_glacier/RGI60-11/RGI60-11.00/RGI60-11.00897

1.4 运行预处理任务

# Run all GIS preprocessing tasks on the glacier
workflow.execute_entity_task(workflow.gis_prepro_tasks, gdirs)

# At this point, the glacier directory contains:
# - outlines.shp (glacier outline)
# - dem.tif (digital elevation model)
# - glacier_grid.json (grid definition)
# - gridded_data.nc (gridded topographic variables)
# - centerlines.pkl (computed centerlines)
# - model_flowlines.pkl (initialized flowlines)

1.5 理解输出文件

让我们检查一下生成了什么：

import os

gdir = gdirs[0]
for fname in os.listdir(gdir.dir):
    fpath = os.path.join(gdir.dir, fname)
    size_kb = os.path.getsize(fpath) / 1024
    print(f"  {fname:<40s} {size_kb:>8.1f} KB")

# Inspect the gridded data
import xarray as xr
with xr.open_dataset(gdir.get_filepath('gridded_data')) as ds:
    print("\nGridded data variables:")
    for vn in ds.data_vars:
        print(f"  {vn}: {ds[vn].shape}, {ds[vn].dtype}")

预期输出（近似值）：

  outlines.shp                                  12.3 KB
  dem.tif                                      156.7 KB
  glacier_grid.json                              0.5 KB
  gridded_data.nc                               45.2 KB
  centerlines.pkl                               23.1 KB
  model_flowlines.pkl                           18.4 KB
  diagnostics.json                               1.2 KB

Gridded data variables:
  topo: (107, 134), float64
  glacier_mask: (107, 134), uint8
  glacier_topo: (107, 134), float64
  aspect: (107, 134), float64
  slope: (107, 134), float64
  dis: (107, 134), float64
  catchment_mask: (107, 134), uint8

你现在已经完成了第一次OGGM预处理运行！冰川目录（GlacierDirectory）已准备好进行气候处理、冰厚反演和动力学模拟。 请参阅第5章（GIS与DEM处理）和第6章（中心线与几何参数化），深入了解每个任务的具体作用。

教程二：使用GlacierDirectory

2.1 深入理解GlacierDirectory

import oggm
from oggm import cfg, utils, workflow

cfg.initialize()
cfg.PATHS['working_dir'] = '/tmp/oggm_tutorial'

# Create a GlacierDirectory from an RGI shapefile
rgi_file = utils.get_demo_file('Hintereisferner_RGI5.shp')
gdir = utils.GlacierDirectory(rgi_file)

# Key attributes
print(f"RGI ID:       {gdir.rgi_id}")
print(f"RGI version:  {gdir.rgi_version}")
print(f"Name:         {gdir.rgi_name}")
print(f"Area (km²):   {gdir.rgi_area_km2:.2f}")
print(f"CenLat:       {gdir.cenlat:.3f}")
print(f"CenLon:       {gdir.cenlon:.3f}")
print(f"TerminusType: {gdir.terminus_type}")
print(f"Data dir:     {gdir.dir}")

2.2 BASENAMES系统

# BASENAMES provides standardized file names within a glacier directory
from oggm import cfg

print("Registered BASENAMES:")
for key, (filename, desc) in cfg.BASENAMES.items():
    print(f"  {key:<30s} → {filename:<30s} ({desc})")

# Get the full path to any registered file type
outline_path = gdir.get_filepath('outlines')
grid_path = gdir.get_filepath('gridded_data')
print(f"\nOutline file: {outline_path}")
print(f"Gridded data: {grid_path}")

# Check if a glacier has a particular file
if gdir.has_file('climate_historical'):
    print("Climate data is present")
else:
    print("Climate data not yet processed")

2.3 读取和写入冰川数据

# Reading gridded data (NetCDF)
import xarray as xr
ds = xr.open_dataset(gdir.get_filepath('gridded_data'))
glacier_topo = ds['glacier_topo'].values  # numpy array
mask = ds['glacier_mask'].values.astype(bool)
print(f"Glacier elevation range: {glacier_topo[mask].min():.0f} - "
      f"{glacier_topo[mask].max():.0f} m a.s.l.")

# Reading pickled objects
centerlines = gdir.read_pickle('centerlines')
print(f"Number of centerlines: {len(centerlines)}")

# Reading JSON
diagnostics = gdir.read_json('diagnostics')
print(f"DEM source: {diagnostics.get('dem_source', 'unknown')}")

# Writing custom data
import json
my_results = {'custom_metric': 42.0, 'rgi_id': gdir.rgi_id}
gdir.write_json(my_results, 'my_custom_results')

# Reading it back
loaded = gdir.read_json('my_custom_results')
print(f"Custom metric: {loaded['custom_metric']}")

2.4 查找GlacierDirectory

# Get all existing glacier directories in the working directory
all_gdirs = workflow.init_glacier_directories(
    rgi_ids=None,  # None = all that exist on disk
    from_prepro_level=0
)
print(f"Found {len(all_gdirs)} glacier directories")

# Filter by RGI region
alps_gdirs = workflow.init_glacier_directories(
    rgi_ids=None,
    from_prepro_level=0,
    reset=False
)
alps_only = [g for g in alps_gdirs if g.rgi_region == '11']
print(f"Alps glaciers: {len(alps_only)}")

教程三：冰川预处理

3.1 预处理工作流

函数 gis_prepro_tasks(gdir) 是一个便捷封装，可一次性运行所有预处理步骤。 为了深入理解，下面我们逐步运行每个步骤：

import oggm
from oggm import cfg, utils, workflow
from oggm.core import gis, centerlines

cfg.initialize()
cfg.PATHS['working_dir'] = '/tmp/oggm_tutorial_prepro'

# Start fresh
gdir = utils.GlacierDirectory(utils.get_demo_file('Hintereisferner_RGI5.shp'))

# Step 1: Define the glacier region on a local grid
gis.define_glacier_region(gdir)
print("1. Glacier region defined")

# Step 2: Process the DEM (reproject, crop, fill gaps)
gis.process_dem(gdir)
print("2. DEM processed")

# Step 3: Compute glacier masks
gis.glacier_masks(gdir)
print("3. Glacier masks computed")

# Step 4: Compute gridded attributes (slope, aspect, distance)
gis.gridded_attributes(gdir)
print("4. Gridded attributes computed")

3.2 中心线计算

# Step 5: Compute centerlines from the glacier mask
centerlines.compute_centerlines(gdir)
print("5. Centerlines computed")

# Visualize centerlines (requires matplotlib)
import matplotlib.pyplot as plt
centerlines.visualize_centerlines(gdir, figsize=(10, 8),
                                   use_flowlines=False)
plt.title(f"Centerlines for {gdir.rgi_id}")
plt.tight_layout()
plt.show()

# Step 6: Compute catchment areas for each centerline
centerlines.catchment_area(gdir)
print("6. Catchment areas computed")

# Step 7: Compute downstream centerline connections
centerlines.compute_downstream_line(gdir)
print("7. Downstream lines computed")

# Step 8: Determine flowline order and tributary structure
centerlines.compute_downstream_bedshape(gdir)
print("8. Downstream bed shape computed")

3.3 流线初始化

# Step 9: Initialize the flowlines
centerlines.initialize_flowlines(gdir)
print("9. Flowlines initialized")

# Inspect the flowlines
fls = gdir.read_pickle('model_flowlines')
for i, fl in enumerate(fls):
    print(f"\nFlowline {i}:")
    print(f"  Points: {len(fl.surface_h)}")
    print(f"  Length: {fl.dis_on_line[-1]:.0f} m")
    print(f"  Max width: {fl.widths.max():.0f} m")
    print(f"  Min elevation: {fl.surface_h.min():.0f} m")
    print(f"  Max elevation: {fl.surface_h.max():.0f} m")
    if hasattr(fl, 'is_tributary'):
        print(f"  Tributary: {fl.is_tributary}")
    if hasattr(fl, 'flows_to'):
        print(f"  Flows to: {fl.flows_to}")

# Visualize the flowlines
centerlines.visualize_centerlines(gdir, figsize=(10, 8),
                                   use_flowlines=True)
plt.title(f"Flowlines for {gdir.rgi_id}")
plt.tight_layout()
plt.show()

教程四：气候数据与物质平衡

4.1 处理气候数据

import oggm
from oggm import cfg, workflow
from oggm.core import climate

cfg.initialize()
cfg.PATHS['working_dir'] = '/tmp/oggm_tutorial_climate'

# Get a preprocessed glacier (skip the GIS steps for brevity)
gdirs = workflow.init_glacier_directories(
    ['RGI60-11.00897'],
    from_prepro_level=2,
    prepro_border=10
)

# Process climate data for the glacier
workflow.execute_entity_task(climate.process_climate_data, gdirs)
print("Climate data processed!")

# Inspect the climate file
import xarray as xr
gdir = gdirs[0]
with xr.open_dataset(gdir.get_filepath('climate_historical')) as ds:
    print(f"Climate data time span: {ds.time.values[0]} to {ds.time.values[-1]}")
    print(f"Variables: {list(ds.data_vars)}")
    print(f"Temperature range: {ds['temp'].min().values:.1f} to "
          f"{ds['temp'].max().values:.1f} °C")
    print(f"Precipitation range: {ds['prcp'].min().values:.1f} to "
          f"{ds['prcp'].max().values:.1f} kg m⁻²")
    print(f"Ref height: {ds.ref_hgt.values:.0f} m")
    print(f"Ref pixel area: {ds.ref_pix_area.values*1e-6:.2f} km²")

4.2 计算物质平衡

# Compute the monthly mass balance at reference heights
mb_df = gdir.get_ref_mb_data()
print("Mass balance at reference heights:")
print(mb_df.describe())

# Plot the mass balance profile
import matplotlib.pyplot as plt
import numpy as np

# Get a typical year's mass balance profile
heights = np.linspace(
    gdir.gridded_data['glacier_topo'].min(),
    gdir.gridded_data['glacier_topo'].max(),
    50
)

# Create a mass balance model
from oggm.core.massbalance import PastMassBalance
mb_model = PastMassBalance(gdir)

# Compute annual MB profile for a specific year
year = 2000
annual_mb = mb_model.get_annual_mb(heights, year=year)

fig, ax = plt.subplots(figsize=(8, 6))
ax.plot(annual_mb, heights, 'b-', linewidth=2)
ax.axvline(x=0, color='k', linestyle='--')
ax.axhline(y=mb_model.get_ela(year=year), color='r', linestyle='--',
           label=f"ELA = {mb_model.get_ela(year=year):.0f} m")
ax.set_xlabel('Annual Mass Balance (kg m⁻² yr⁻¹)')
ax.set_ylabel('Elevation (m a.s.l.)')
ax.set_title(f'Mass Balance Profile — {year}')
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

4.3 校准物质平衡模型

# OGGM automatically calibrates melt_f (temperature sensitivity)
# during process_climate_data, but you can inspect the results:

calib = gdir.read_json('mb_calib')
print("Mass balance calibration parameters:")
for key, val in calib.items():
    print(f"  {key}: {val}")

# The key parameters:
# - melt_f: temperature melt factor (kg m⁻² K⁻¹ day⁻¹)
# - temp_bias: temperature correction applied to match geodetic MB
# - prcp_fac: precipitation correction factor

# Compute the specific mass balance for a period
from oggm.core.massbalance import MultipleFlowlineMassBalance
mb_multi = MultipleFlowlineMassBalance(gdir)

years = np.arange(1980, 2010)
s_mb_list = []
for yr in years:
    ela = mb_multi.get_ela(year=yr)
    s_mb = mb_multi.get_specific_mb(year=yr)
    s_mb_list.append(s_mb)

print(f"\nMean specific MB (1980-2010): {np.mean(s_mb_list):.2f} mm w.e. yr⁻¹")
print(f"Mean ELA (1980-2010): {np.mean([mb_multi.get_ela(year=y) for y in years]):.0f} m")

教程五：冰厚反演

5.1 运行反演

import oggm
from oggm import cfg, workflow
from oggm.core import inversion

cfg.initialize()
cfg.PATHS['working_dir'] = '/tmp/oggm_tutorial_inversion'

# Start from a preprocessed glacier with climate data
gdirs = workflow.init_glacier_directories(
    ['RGI60-11.00897'],
    from_prepro_level=3,  # needs climate data
    prepro_border=10
)

# Run the inversion for each glacier
workflow.execute_entity_task(inversion.mass_conservation_inversion, gdirs)
print("Inversion complete!")

# Inspect the inversion output
gdir = gdirs[0]
inv_output = gdir.read_pickle('inversion_output')
print("\nInversion output keys:", list(inv_output.keys()))
print(f"Inverted volume: {inv_output['volume']*1e-9:.4f} km³")
print(f"Mean thickness: {inv_output['volume']/gdir.rgi_area_m2:.1f} m")

5.2 理解反演参数

# Key parameters affecting inversion results:
print("Inversion parameters:")
print(f"  Glen A:     {cfg.PARAMS['glen_a']:.2e} s⁻¹ Pa⁻³")
print(f"  Glen n:     {cfg.PARAMS['glen_n']}")
print(f"  Sliding fs: {cfg.PARAMS['fs']}")
print(f"  Inversion glen_a: {cfg.PARAMS['inversion_glen_a']:.2e}")
print(f"  Inversion fs:     {cfg.PARAMS['inversion_fs']}")
print(f"  mixed_min_shape:  {cfg.PARAMS['mixed_min_shape']}")

# The inversion solves for ice thickness H such that:
#   flux(H) = integrated_apparent_mass_balance
# where flux = u_bar * H * w and u_bar follows the SIA
# with deformation + sliding contributions.

5.3 可视化反演冰厚

import matplotlib.pyplot as plt
import numpy as np

# Plot the thickness profile along each flowline
fls = gdir.read_pickle('inversion_input')
inv = gdir.read_pickle('inversion_output')

fig, axes = plt.subplots(len(fls), 1, figsize=(12, 4 * len(fls)),
                          squeeze=False)
for i, fl in enumerate(fls):
    ax = axes[i, 0]
    x = fl.dis_on_line / 1000  # km
    thick = inv['thickness'][i]
    # Surface = bed + thickness
    bed = fl.bed_h
    surface = fl.surface_h
    computed_surface = bed + thick

    ax.fill_between(x, bed, surface, alpha=0.3, color='blue',
                     label='Measured surface')
    ax.fill_between(x, bed, computed_surface, alpha=0.5, color='red',
                     label='Inverted thickness')
    ax.plot(x, bed, 'k-', linewidth=1.5, label='Bed')
    ax.set_xlabel('Distance along flowline (km)')
    ax.set_ylabel('Elevation (m a.s.l.)')
    ax.set_title(f'Flowline {i}')
    ax.legend(loc='lower left')
    ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

5.4 过滤反演输出

# The raw inversion can have noise. Apply filtering:
inversion.filter_inversion_output(gdir)

# Re-read and compare
inv_filtered = gdir.read_pickle('inversion_output')
print(f"Volume before filter: {inv['volume']*1e-9:.4f} km³")
print(f"Volume after filter:  {inv_filtered['volume']*1e-9:.4f} km³")

# Sensitivity test: vary inversion parameters
import copy
original_glen_a = cfg.PARAMS['inversion_glen_a']

results = {}
for factor in [0.5, 1.0, 2.0]:
    cfg.PARAMS['inversion_glen_a'] = original_glen_a * factor
    # In a real sensitivity experiment, reset or overwrite the inversion
    # output before re-running mass_conservation_inversion.
    vol = inv_filtered['volume'] * 1e-9  # placeholder
    results[factor] = vol
    print(f"  glen_a factor {factor}: volume = {vol:.4f} km³")

cfg.PARAMS['inversion_glen_a'] = original_glen_a  # restore

教程六：冰川演化

6.1 运行历史模拟

import oggm
from oggm import cfg, workflow
from oggm.core import flowline, massbalance

cfg.initialize()
cfg.PATHS['working_dir'] = '/tmp/oggm_tutorial_evolution'

# Get a fully preprocessed glacier (prepro_level >= 3)
gdirs = workflow.init_glacier_directories(
    ['RGI60-11.00897'],
    from_prepro_level=3,
    prepro_border=10
)

gdir = gdirs[0]

# Create a mass balance model (historical climate)
mb_model = massbalance.PastMassBalance(gdir)

# Initialize the flowline model with the inverted bed
fls = gdir.read_pickle('model_flowlines')
model = flowline.FluxBasedModel(
    flowlines=fls,
    mb_model=mb_model,
    y0=2003,           # start year (RGI date)
    glen_a=cfg.PARAMS['glen_a'],
    fs=cfg.PARAMS['fs'],
)

# Run from 2003 to 2020
print("Running historical simulation (2003-2020)...")
model.run_until(2020)

# Inspect results
diag = model.get_diagnostics()
print(f"\nFinal year diagnostics ({model.yr:.0f}):")
print(f"  Volume:        {model.volume_km3:.4f} km³")
print(f"  Area:          {model.area_km2:.2f} km²")
print(f"  Length:        {model.length_m:.0f} m")
print(f"  ELA:           {diag.get('ela', 'N/A')}")
print(f"  Spec. MB:      {model.get_specific_mb():.2f} mm w.e. yr⁻¹")

6.2 存储模型诊断数据

# Run and store diagnostics at each year
from oggm.core.flowline import flowline_model_run

# Re-initialize
model = flowline.FluxBasedModel(
    flowlines=fls,
    mb_model=mb_model,
    y0=2003,
    glen_a=cfg.PARAMS['glen_a'],
    fs=cfg.PARAMS['fs'],
)

# Run with diagnostics
diag_path = gdir.get_filepath('model_diagnostics')
flowline_model_run(
    gdir, model,
    ys=2003, ye=2020,
    store_diagnostics=True,
    output_filesuffix='_historical'
)

# Read the diagnostics
import xarray as xr
with xr.open_dataset(
    gdir.get_filepath('model_diagnostics',
                      filesuffix='_historical')
) as ds_diag:
    print(f"Diagnostic variables: {list(ds_diag.data_vars)}")
    print(f"Time steps: {len(ds_diag.time)}")

    # Plot volume evolution
    import matplotlib.pyplot as plt
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))

    ax1.plot(ds_diag.time, ds_diag.volume_m3 * 1e-9)
    ax1.set_xlabel('Year')
    ax1.set_ylabel('Volume (km³)')
    ax1.set_title('Volume Evolution')
    ax1.grid(True, alpha=0.3)

    ax2.plot(ds_diag.time, ds_diag.area_m2 * 1e-6)
    ax2.set_xlabel('Year')
    ax2.set_ylabel('Area (km²)')
    ax2.set_title('Area Evolution')
    ax2.grid(True, alpha=0.3)

    plt.tight_layout()
    plt.show()

6.3 使用GCM数据进行未来预测

from oggm.core import climate

# Compute the climate delta (GCM - historical) for a future period
# This creates a "delta method" scenario
climate.historical_delta_method(gdir)

# Create a mass balance model for the future
mb_future = massbalance.ConstantMassBalance(gdir)

# Initialize a new model from the present-day state
model_future = flowline.FluxBasedModel(
    flowlines=fls,
    mb_model=mb_future,
    y0=2020,
    glen_a=cfg.PARAMS['glen_a'],
    fs=cfg.PARAMS['fs'],
    is_tidewater=gdir.is_tidewater,
)

# Run to 2100
print("Running future projection (2020-2100)...")
model_future.run_until(2100)

diag = model_future.get_diagnostics()
print(f"\nProjected state at 2100:")
print(f"  Volume:        {model_future.volume_km3:.4f} km³")
print(f"  Area:          {model_future.area_km2:.2f} km²")
print(f"  Length change: {model_future.length_m - gdir.rgi_length_m:.0f} m")

6.4 理解模型诊断数据

# Full list of available diagnostics from FluxBasedModel
print("Available diagnostics:")
common_diags = [
    'volume_m3', 'area_m2', 'length_m',
    'ela_m', 'calving_m3',
    'dmdtda_mmwea',  # specific mass balance
]
for d in common_diags:
    print(f"  {d}")

# Annual mass balance breakdown
print("\nModel.get_diagnostics() return values:")
full_diag = model_future.get_diagnostics()
for k, v in full_diag.items():
    if isinstance(v, (int, float)):
        print(f"  {k}: {v:.2f}")
    else:
        print(f"  {k}: {type(v).__name__}")

教程七：处理多个冰川

7.1 处理区域冰川目录（GlacierDirectory）

import oggm
from oggm import cfg, workflow, tasks

cfg.initialize()
cfg.PATHS['working_dir'] = '/tmp/oggm_tutorial_regional'

# Define a set of glaciers (e.g., a handful from the Alps)
rgi_ids = [
    'RGI60-11.00897',  # Hintereisferner
    'RGI60-11.00778',  # Vernagtferner
    'RGI60-11.00002',  # Pasterze
    'RGI60-11.02692',  # Findelengletscher
    'RGI60-11.02709',  # Gornergletscher
]

print(f"Processing {len(rgi_ids)} glaciers...")

# Initialize all glacier directories
gdirs = workflow.init_glacier_directories(
    rgi_ids,
    from_prepro_level=0,
    prepro_border=10,
    reset=True
)
print(f"Initialized {len(gdirs)} glacier directories")

# Run preprocessing pipeline (enable parallel processing for speed)
cfg.PARAMS['use_multiprocessing'] = True
workflow.execute_entity_task(workflow.gis_prepro_tasks, gdirs)
print("GIS preprocessing complete!")

# Quick status check
for gdir in gdirs:
    print(f"  {gdir.rgi_id}: {gdir.rgi_name}, "
          f"n_flowlines={len(gdir.read_pickle('model_flowlines'))}")

7.2 使用execute_entity_task并行处理

import time

# Run climate processing in parallel
print("Processing climate data (parallel)...")
t0 = time.time()
workflow.execute_entity_task(tasks.process_climate_data, gdirs)
print(f"Climate processing done in {time.time() - t0:.1f}s")

# Run inversion in parallel
print("Running mass-conservation inversion (parallel)...")
t0 = time.time()
workflow.execute_entity_task(tasks.mass_conservation_inversion, gdirs)
print(f"Inversion done in {time.time() - t0:.1f}s")

# Monitor progress
workflow.execute_entity_task(
    tasks.mass_conservation_inversion, gdirs,
    print_log=True  # print progress to stdout
)

7.3 汇总区域结果

import pandas as pd
import numpy as np

# Compile results from all glacier directories
results = []
for gdir in gdirs:
    try:
        inv = gdir.read_pickle('inversion_output')
        diag = gdir.read_json('diagnostics')
        calib = gdir.read_json('mb_calib')

        results.append({
            'rgi_id': gdir.rgi_id,
            'name': gdir.rgi_name,
            'area_km2': gdir.rgi_area_km2,
            'volume_km3': inv['volume'] * 1e-9,
            'mean_thick_m': inv['volume'] / gdir.rgi_area_m2,
            'n_flowlines': len(gdir.read_pickle('model_flowlines')),
            'melt_f': calib.get('melt_f', np.nan),
            'temp_bias': calib.get('temp_bias', np.nan),
            'prcp_fac': calib.get('prcp_fac', np.nan),
            'dem_source': diag.get('dem_source', 'unknown'),
        })
    except Exception as e:
        print(f"Failed for {gdir.rgi_id}: {e}")

df = pd.DataFrame(results)
print("\nRegional Summary:")
print(f"  Total glaciers:  {len(df)}")
print(f"  Total area:      {df['area_km2'].sum():.1f} km²")
print(f"  Total volume:    {df['volume_km3'].sum():.3f} km³")
print(f"  Mean thickness:  {df['mean_thick_m'].mean():.1f} m")
print(f"  Mean melt_f:    {df['melt_f'].mean():.2f}")
print("\nPer-glacier results:")
print(df.to_string(index=False))

7.4 错误处理和稳健处理

# Use "continue_on_error" for robust regional processing
cfg.PARAMS['continue_on_error'] = True

# This will log errors but not crash if one glacier fails
workflow.execute_entity_task(
    tasks.run_from_climate_data, gdirs,
    ys=2003, ye=2020,
)

# Check which glaciers succeeded
for gdir in gdirs:
    if gdir.has_file('model_diagnostics'):
        print(f"  [OK]  {gdir.rgi_id}")
    else:
        print(f"  [FAIL] {gdir.rgi_id} — check log for details")

# Reset the parameter
cfg.PARAMS['continue_on_error'] = False

教程八：高级专题

8.1 动态预热（dynamic spinup）

# Dynamic spinup searches for a temperature bias and recent-past
# initialization that matches the glacier area or volume near the RGI date.
from oggm.core import dynamic_spinup

gdir = gdirs[0]  # use a preprocessed gdir

# Run a short dynamic spinup and store the initialized model geometry
model = dynamic_spinup.run_dynamic_spinup(
    gdir,
    spinup_period=20,
    min_spinup_period=10,
    minimise_for='area',
    precision_percent=1,
    output_filesuffix='_dynamic_spinup',
)

print("Running spinup...")
# The spinup runs automatically before the main simulation
model.run_until(2020)

print(f"Post-spinup volume: {model.volume_km3:.4f} km³")
diag = gdir.get_diagnostics()
print(f"Spinup period: {diag['dynamic_spinup_period']} years")
print(f"Temperature bias: {diag['temp_bias_dynamic_spinup']:.2f} °C")

8.2 入海型冰川与崩解

import xarray as xr

# Tidewater glaciers require additional setup.
# Start from an already prepared list of GlacierDirectory objects, e.g.
# gdirs = workflow.init_glacier_directories(rgi_ids, from_prepro_level=3, ...)
tw_gdirs = [g for g in gdirs if g.is_tidewater]
if len(tw_gdirs) > 0:
    tw_gdir = tw_gdirs[0]
    print(f"Tidewater glacier: {tw_gdir.rgi_id} "
          f"({tw_gdir.rgi_name})")
    print(f"Terminus type: {tw_gdir.terminus_type}")

    # Configure k-calving. OGGM v1.6.3 uses 0.6 yr^-1 by default.
    cfg.PARAMS['use_kcalving_for_run'] = True
    cfg.PARAMS['calving_k'] = 0.6

    # Recompute calving-aware inversion diagnostics, initialize the model,
    # and run with the standard OGGM task interface.
    workflow.execute_entity_task(
        tasks.find_inversion_calving_from_any_mb, [tw_gdir]
    )
    workflow.execute_entity_task(
        tasks.init_present_time_glacier, [tw_gdir]
    )
    workflow.execute_entity_task(
        tasks.run_from_climate_data, [tw_gdir], ys=2000, ye=2020
    )

    with xr.open_dataset(tw_gdir.get_filepath('model_diagnostics')) as ds:
        print(float(ds.volume_m3.isel(time=-1)) / 1e9, "km3")
        print(float(ds.length_m.isel(time=-1)), "m")

8.3 合并支流冰川

from oggm.core import flowline

# Complex glaciers have multiple centerlines with tributaries
# You can control how tributaries are handled:

# Option 1: Independent flowlines (default)
cfg.PARAMS['use_multiple_flowlines'] = True

# Option 2: Merge tributaries into one main flowline
cfg.PARAMS['use_multiple_flowlines'] = False

# Inspect the tributary structure
fls = gdir.read_pickle('model_flowlines')
for i, fl in enumerate(fls):
    info = []
    if hasattr(fl, 'is_tributary') and fl.is_tributary:
        info.append('tributary')
    if hasattr(fl, 'flows_to'):
        info.append(f'flows_to={fl.flows_to}')
    if hasattr(fl, 'is_rectangular'):
        info.append(f'rectangular={fl.is_rectangular}')
    print(f"  Flowline {i}: {', '.join(info)}")

8.4 使用Delta方法进行气候预测

# The delta method computes climate change signals from GCMs
# and applies them as anomalies to the historical climate

from oggm.core import climate
import xarray as xr

# Step 1: Compute the delta for a given GCM and scenario
# (This is normally done during process_climate_data)
climate.historical_delta_method(gdir)

# Step 2: Inspect the delta
delta_file = gdir.get_filepath(
    'climate_historical',
    filesuffix='_delta'
)
print(f"Delta file: {delta_file}")

# Step 3: Create a future mass balance model
from oggm.core.massbalance import ConstantMassBalance
mb_future = ConstantMassBalance(
    gdir, y0=2003,
    halfsize=15,        # window half-size for climate averaging
    filename='climate_historical',
    input_filesuffix='_delta',
)

# The delta signal is already baked into the climate file,
# so the ConstantMassBalance model can use it directly.
# For GCM-driven runs, you compute deltas from raw GCM output:
# climate.process_gcm_data(gdir, filesuffix='_GCM_output')

8.5 交叉验证与敏感性分析

# Systematic parameter sensitivity test
import copy
from oggm.core.flowline import FluxBasedModel

baseline_params = {
    'glen_a': cfg.PARAMS['glen_a'],
    'fs': cfg.PARAMS['fs'],
}
fls = gdir.read_pickle('model_flowlines')
mb_model = massbalance.PastMassBalance(gdir)

sensitivity_results = {}
for glen_factor in [0.5, 1.0, 2.0]:
    for fs_value in [0.0, cfg.PARAMS['fs']]:
        # Clone flowlines to avoid in-place modification
        import copy
        fls_copy = copy.deepcopy(fls)

        model = FluxBasedModel(
            flowlines=fls_copy,
            mb_model=mb_model,
            y0=2003,
            glen_a=baseline_params['glen_a'] * glen_factor,
            fs=fs_value,
        )
        model.run_until(2020)
        key = f"glen={glen_factor:.1f}, fs={fs_value:.0f}"
        sensitivity_results[key] = model.volume_km3

print("\nSensitivity Analysis (final volume, km³):")
for key, vol in sensitivity_results.items():
    print(f"  {key}: {vol:.4f}")

教程九：常见错误与定位方法

# Minimal diagnostic checklist for one glacier directory
for name in ['gridded_data', 'climate_historical', 'mb_calib',
             'inversion_output', 'model_flowlines', 'model_diagnostics']:
    print(f"{name:25s}", gdir.has_file(name))

print(gdir.rgi_id, gdir.rgi_area_km2, gdir.rgi_region)
print(gdir.get_diagnostics())