diff --git a/.github/workflows/test_suite.yml b/.github/workflows/test_suite.yml index 30f9dfeb..2a15a264 100644 --- a/.github/workflows/test_suite.yml +++ b/.github/workflows/test_suite.yml @@ -57,11 +57,14 @@ jobs: - name: 'Clone the PyGEM-notebooks repo' run: | - # Use PyGEM-notebook:main for master branch and PyGEM-notebooks:dev otherwise BRANCH=${GITHUB_REF#refs/heads/} - git clone --depth 1 --branch $([[ "$BRANCH" == "master" ]] && echo "main" || echo "dev") \ - https://github.com/pygem-community/PyGEM-notebooks.git - echo "PYGEM_NOTEBOOKS_DIRPATH=$(pwd)/PyGEM-notebooks" >> $GITHUB_ENV + if [ "$BRANCH" = "master" ]; then + NOTEBOOK_BRANCH="main" + else + NOTEBOOK_BRANCH="dev" + fi + git clone --depth 1 --branch "$NOTEBOOK_BRANCH" https://github.com/pygem-community/PyGEM-notebooks.git + echo "PYGEM_NOTEBOOKS_DIRPATH=$(pwd)/PyGEM-notebooks" >> "$GITHUB_ENV" - name: 'Run tests' run: | diff --git a/README.md b/README.md index 94efdbc4..a77a2fd1 100755 --- a/README.md +++ b/README.md @@ -27,6 +27,8 @@ To support model testing and demonstration, a suite of Jupyter notebooks can be Citation +   + DOI diff --git a/docs/conf.py b/docs/conf.py index db2f2529..f881f30e 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -56,11 +56,11 @@ html_theme = 'sphinx_book_theme' -html_static_path = ['_static'] html_theme_options = { 'repository_url': 'https://github.com/PyGEM-Community/PyGEM', 'use_repository_button': True, - 'show_nav_level': 2, - 'navigation_depth': 3, + 'show_nav_level': 1, + 'navigation_depth': 4, + 'toc_title': 'On this page', } diff --git a/docs/configuration.md b/docs/configuration.md index a5ec8df0..14ba06f5 100644 --- a/docs/configuration.md +++ b/docs/configuration.md @@ -59,8 +59,8 @@ The user information is strictly for bookkeeping purposes. This information is s | Variable | Type | Comment/Note | | :--- | :--- | :--- | -| `setup.rgi_region01` | `list` of `integers` | List of RGI region numbers to include | -| `setup.rgi_region02` | `list` of `integers` or `"all"` | List of RGI region numbers or `"all"` to include all regions | +| `setup.rgi_region01` | `list` of `integers` | List of RGI order 01 region numbers to include | +| `setup.rgi_region02` | `list` of `integers` or `'all'` | List of RGI order 02 subregion numbers to include, or `'all'` to include all subregions | | `setup.glac_no_skip` | `null` or `list` of `strings` | Glacier numbers to skip | | `setup.glac_no` | `null` or `list` of `strings` | List of RGI glacier numbers (e.g., `1.00570`) | | `setup.min_glac_area_km2` | `float` | Minimum glacier area (km$^{2}$) threshold | @@ -79,7 +79,7 @@ Set glac_no will always overwrite the rgi regions, so if you want to use the rgi | Variable | Type | Comment/Note | | :--- | :--- | :--- | | `oggm.base_url` | `string` | URL for OGGM glacier directories | -| `oggm.logging_level` | `string` | Log level (`DEBUG`, `INFO`, `WARNING`, `ERROR`, `WORKFLOW`, `CRITICAL`) | +| `oggm.logging_level` | `string` | Log level (`'DEBUG'`, `'INFO'`, `'WARNING'`, `'ERROR'`, `'WORKFLOW'`, `'CRITICAL'`) | | `oggm.border` | `integer` | Border size (options: `10`, `80`, `160`, `240`) | | `oggm.oggm_gdir_relpath` | `string` | Relative path to OGGM glacier directories | | `oggm.overwrite_gdirs` | `boolean` | Overwrite glacier directories if they exist | @@ -91,22 +91,19 @@ Set glac_no will always overwrite the rgi regions, so if you want to use the rgi | Variable | Type | Comment/Note | | :--- | :--- | :--- | -| `climate.ref_gcm_name` | `string` | Reference climate dataset (`ERA5`, etc.) | -| `climate.ref_startyear` | `integer` | Start year for reference dataset | -| `climate.ref_endyear` | `integer` | End year for reference dataset | -| `climate.ref_wateryear` | `string` | Type of year (`calendar`, `hydro`, `custom`) | -| `climate.ref_spinupyears` | `integer` | Number of spin-up years | -| `climate.gcm_name` | `string` | GCM dataset used for simulations | -| `climate.scenario` | `null` or `string` | Climate scenario | -| `climate.gcm_startyear` | `integer` | Start year for GCM dataset | -| `climate.gcm_endyear` | `integer` | End year for GCM dataset | -| `climate.gcm_wateryear` | `string` | Year type (`calendar`, `hydro`, `custom`) | -| `climate.constantarea_years` | `integer` | Years to keep glacier area constant | -| `climate.gcm_spinupyears` | `integer` | Number of spin-up years for simulation | -| `climate.hindcast` | `boolean` | Enable or disable hindcasting | -| `climate.paths` | `string` | Relative filepaths, filenames and variable names for climate datasets | +| `climate.ref_climate_name` | `string` | Historical reference climate dataset (e.g.`'ERA5'`) | +| `climate.ref_startyear` | `integer` | Reference period start year (inclusive, default is `2000`) | +| `climate.ref_endyear` | `integer` | Reference period end year (inclusive, default is `2019`) | +| `climate.ref_wateryear` | `string` | Reference period year type (`'calendar'`, `'hydro'`, `'custom'`) | +| `climate.sim_climate_name` | `string` | Climate dataset used for simulations (can be reference climate or GCM, e.g. `ERA5` or `CESM2`) | +| `climate.sim_climate_scenario` | `null` or `string` | Climate emission scenario (e.g. `'ssp245', 'rcp65'`) | +| `climate.sim_startyear` | `integer` | Start year for model simulation (inclusive, e.g. `2000`) | +| `climate.sim_endyear` | `integer` | End year for model simulation (inclusinve, e.g. `2100`) | +| `climate.sim_wateryear` | `string` | Year type (`'calendar'`, `'hydro'`, `'custom'`) | +| `climate.constantarea_years` | `integer` | Number of years from `climate.ref_startyear` to keep glacier area constant | +| `climate.paths` | `string` | Relative filepaths, filenames and variable names for climate datasets (relative to `root`) | ```{note} -The hindcast option will flip the climate data array so 1960-2000 would run 2000-1960 ensuring that glacier area at 2000 is correct; however, due to nonlinearities a run starting at 1960 would not provide the same mass change and area at 2000.. +`ref` (short for 'reference') refers to the calibration period (will typically vary between 1980-present). `sim` (short for 'simulation') refers to simulation period. ``` --- @@ -115,32 +112,32 @@ The hindcast option will flip the climate data array so 1960-2000 would run 2000 | Variable | Type | Comment/Note | | :--- | :--- | :--- | -| `calib.option_calibration` | `string` | Calibration option ('emulator', 'MCMC', 'HH2015', 'HH2015mod', 'null') | -| `calib.priors_reg_fn` | `string` | Prior distribution (specify filename, relative to `root/Output/calibration/`, or set to `null`) | +| `calib.option_calibration` | `string` | Calibration option (`'emulator'`, `'MCMC'`, `'HH2015'`, `'HH2015mod'`, `'null'`) | +| `calib.priors_reg_fn` | `string` | Prior distribution of regional model parameters (specify filename, relative to `root/Output/calibration/`, or set to `'null'`) | ### HH2015 Parameters | Variable | Type | Comment/Note | | :--- | :--- | :--- | -| `calib.HH2015_params.tbias_init` | `float` | Initial temperature bias | -| `calib.HH2015_params.tbias_step` | `float` | Temperature bias step size | -| `calib.HH2015_params.kp_init` | `float` | Initial precipitation factor | -| `calib.HH2015_params.kp_bndlow` | `float` | Lower bound for precipitation factor | -| `calib.HH2015_params.kp_bndhigh` | `float` | Upper bound for precipitation factor | -| `calib.HH2015_params.ddfsnow_init` | `float` | Initial degree-day factor for snow | -| `calib.HH2015_params.ddfsnow_bndlow` | `float` | Lower bound for degree-day factor for snow | -| `calib.HH2015_params.ddfsnow_bndhigh` | `float` | Upper bound for degree-day factor for snow | +| `calib.HH2015_params.tbias_init` | `float` | Initial temperature bias [K] | +| `calib.HH2015_params.tbias_step` | `float` | Temperature bias step size [K] | +| `calib.HH2015_params.kp_init` | `float` | Initial precipitation factor [-] | +| `calib.HH2015_params.kp_bndlow` | `float` | Lower bound for precipitation factor [-] | +| `calib.HH2015_params.kp_bndhigh` | `float` | Upper bound for precipitation factor [-] | +| `calib.HH2015_params.ddfsnow_init` | `float` | Initial degree-day factor for snow [mm w.e. d$^{-1}$ K$^{-1}$] | +| `calib.HH2015_params.ddfsnow_bndlow` | `float` | Lower bound for degree-day factor for snow [mm w.e. d$^{-1}$ K$^{-1}$] | +| `calib.HH2015_params.ddfsnow_bndhigh` | `float` | Upper bound for degree-day factor for snow [mm w.e. d$^{-1}$ K$^{-1}$] | ### HH2015mod Parameters | Variable | Type | Comment/Note | | :--- | :--- | :--- | -| `calib.HH2015mod_params.tbias_init` | `float` | Initial temperature bias | -| `calib.HH2015mod_params.tbias_step` | `float` | Temperature bias step size | -| `calib.HH2015mod_params.kp_init` | `float` | Initial precipitation factor | -| `calib.HH2015mod_params.kp_bndlow` | `float` | Lower bound for precipitation factor | -| `calib.HH2015mod_params.kp_bndhigh` | `float` | Upper bound for precipitation factor | -| `calib.HH2015mod_params.ddfsnow_init`| `float` | Initial degree-day factor for snow | +| `calib.HH2015mod_params.tbias_init` | `float` | Initial temperature bias [K] | +| `calib.HH2015mod_params.tbias_step` | `float` | Temperature bias step size [K] | +| `calib.HH2015mod_params.kp_init` | `float` | Initial precipitation factor [-] | +| `calib.HH2015mod_params.kp_bndlow` | `float` | Lower bound for precipitation factor [-] | +| `calib.HH2015mod_params.kp_bndhigh` | `float` | Upper bound for precipitation factor [-] | +| `calib.HH2015mod_params.ddfsnow_init`| `float` | Initial degree-day factor for snow [mm w.e. d$^{-1}$ K$^{-1}$] | | `calib.HH2015mod_params.method_opt` | `string` | SciPy optimization scheme ('SLSQP' or 'L-BFGS-B') | | `calib.HH2015mod_params.params2opt` | `list` of `strings` | Parameters to optimize (e.g., 'tbias', 'kp') | | `calib.HH2015mod_params.ftol_opt` | `float` | Tolerance for SciPy optimization scheme | @@ -153,22 +150,22 @@ The hindcast option will flip the climate data array so 1960-2000 would run 2000 | `calib.emulator_params.emulator_sims` | `integer` | Number of simulations to develop the emulator | | `calib.emulator_params.overwrite_em_sims` | `boolean` | Whether to overwrite emulator simulations | | `calib.emulator_params.opt_hh2015_mod` | `boolean` | Option to also perform HH2015_mod calibration using emulator | -| `calib.emulator_params.tbias_step` | `float` | Temperature bias step size | -| `calib.emulator_params.tbias_init` | `float` | Initial temperature bias | -| `calib.emulator_params.kp_init` | `float` | Initial precipitation factor | -| `calib.emulator_params.kp_bndlow` | `float` | Lower bound for precipitation factor | -| `calib.emulator_params.kp_bndhigh` | `float` | Upper bound for precipitation factor | -| `calib.emulator_params.ddfsnow_init` | `float` | Initial degree-day factor for snow | +| `calib.emulator_params.tbias_step` | `float` | Temperature bias step size [K] | +| `calib.emulator_params.tbias_init` | `float` | Initial temperature bias [K] | +| `calib.emulator_params.kp_init` | `float` | Initial precipitation factor [-] | +| `calib.emulator_params.kp_bndlow` | `float` | Lower bound for precipitation factor [-] | +| `calib.emulator_params.kp_bndhigh` | `float` | Upper bound for precipitation factor [-] | +| `calib.emulator_params.ddfsnow_init` | `float` | Initial degree-day factor for snow [mm w.e. d$^{-1}$ K$^{-1}$] | | `calib.emulator_params.option_areaconstant` | `boolean` | Option to keep area constant or evolve | | `calib.emulator_params.tbias_disttype` | `string` | Temperature bias distribution type ('truncnormal', 'uniform') | -| `calib.emulator_params.tbias_sigma` | `float` | Temperature bias standard deviation | -| `calib.emulator_params.kp_gamma_alpha` | `float` | Precipitation factor gamma distribution alpha | -| `calib.emulator_params.kp_gamma_beta` | `float` | Precipitation factor gamma distribution beta | +| `calib.emulator_params.tbias_sigma` | `float` | Temperature bias standard deviation [K] | +| `calib.emulator_params.kp_gamma_alpha` | `float` | Precipitation factor gamma distribution alpha [-] | +| `calib.emulator_params.kp_gamma_beta` | `float` | Precipitation factor gamma distribution beta [-] | | `calib.emulator_params.ddfsnow_disttype` | `string` | Degree-day factor for snow distribution type ('truncnormal') | -| `calib.emulator_params.ddfsnow_mu` | `float` | Degree-day factor for snow mean | -| `calib.emulator_params.ddfsnow_sigma` | `float` | Degree-day factor for snow standard deviation | -| `calib.emulator_params.ddfsnow_bndlow` | `float` | Lower bound for degree-day factor for snow | -| `calib.emulator_params.ddfsnow_bndhigh` | `float` | Upper bound for degree-day factor for snow | +| `calib.emulator_params.ddfsnow_mu` | `float` | Degree-day factor for snow mean [mm w.e. d$^{-1}$ K$^{-1}$] | +| `calib.emulator_params.ddfsnow_sigma` | `float` | Degree-day factor for snow standard deviation [mm w.e. d$^{-1}$ K$^{-1}$] | +| `calib.emulator_params.ddfsnow_bndlow` | `float` | Lower bound for degree-day factor for snow [mm w.e. d$^{-1}$ K$^{-1}$] | +| `calib.emulator_params.ddfsnow_bndhigh` | `float` | Upper bound for degree-day factor for snow [mm w.e. d$^{-1}$ K$^{-1}$] | | `calib.emulator_params.method_opt` | `string` | SciPy optimization scheme ('SLSQP' or 'L-BFGS-B') | | `calib.emulator_params.params2opt` | `list` of `strings` | Parameters to optimize (e.g., 'tbias', 'kp') | | `calib.emulator_params.ftol_opt` | `float` | Tolerance for SciPy optimization scheme | @@ -180,47 +177,47 @@ The hindcast option will flip the climate data array so 1960-2000 would run 2000 | :--- | :--- | :--- | | `calib.MCMC_params.option_use_emulator` | `boolean` | Switch to emulator instead of full mass balance model | | `calib.MCMC_params.emulator_sims` | `integer` | Number of emulator simulations | -| `calib.MCMC_params.tbias_step` | `float` | Temperature bias step size | -| `calib.MCMC_params.tbias_stepsmall` | `float` | Small temperature bias step size | +| `calib.MCMC_params.tbias_step` | `float` | Temperature bias step size [K] | +| `calib.MCMC_params.tbias_stepsmall` | `float` | Small temperature bias step size [K] | | `calib.MCMC_params.option_areaconstant` | `boolean` | Option to keep area constant or evolve | -| `calib.MCMC_params.mcmc_step` | `float` | MCMC step size (in terms of standard deviation) | +| `calib.MCMC_params.mcmc_step` | `float` | MCMC step size (in terms of standard deviation of each parameters prior distribution) | | `calib.MCMC_params.n_chains` | `integer` | Number of MCMC chains (min 1, max 3) | | `calib.MCMC_params.mcmc_sample_no` | `integer` | Number of MCMC steps | | `calib.MCMC_params.mcmc_burn_pct` | `float` | Percentage of steps to burn-in (0 records all steps) | | `calib.MCMC_params.thin_interval` | `integer` | Thin interval for reducing file size | | `calib.MCMC_params.ddfsnow_disttype` | `string` | Degree-day factor for snow distribution type ('truncnormal', 'uniform') | -| `calib.MCMC_params.ddfsnow_mu` | `float` | Degree-day factor for snow mean | -| `calib.MCMC_params.ddfsnow_sigma` | `float` | Degree-day factor for snow standard deviation | -| `calib.MCMC_params.ddfsnow_bndlow` | `float` | Lower bound for degree-day factor for snow | -| `calib.MCMC_params.ddfsnow_bndhigh` | `float` | Upper bound for degree-day factor for snow | +| `calib.MCMC_params.ddfsnow_mu` | `float` | Degree-day factor for snow mean [mm w.e. d$^{-1}$ K$^{-1}$] | +| `calib.MCMC_params.ddfsnow_sigma` | `float` | Degree-day factor for snow standard deviation [mm w.e. d$^{-1}$ K$^{-1}$] | +| `calib.MCMC_params.ddfsnow_bndlow` | `float` | Lower bound for degree-day factor for snow [mm w.e. d$^{-1}$ K$^{-1}$] | +| `calib.MCMC_params.ddfsnow_bndhigh` | `float` | Upper bound for degree-day factor for snow [mm w.e. d$^{-1}$ K$^{-1}$] | | `calib.MCMC_params.kp_disttype` | `string` | Precipitation factor distribution type ('gamma', 'lognormal', 'uniform') | | `calib.MCMC_params.tbias_disttype` | `string` | Temperature bias distribution type ('normal', 'truncnormal', 'uniform') | -| `calib.MCMC_params.tbias_mu` | `float` | Temperature bias mean | -| `calib.MCMC_params.tbias_sigma` | `float` | Temperature bias standard deviation | -| `calib.MCMC_params.tbias_bndlow` | `float` | Lower bound for temperature bias | -| `calib.MCMC_params.tbias_bndhigh` | `float` | Upper bound for temperature bias | -| `calib.MCMC_params.kp_gamma_alpha` | `float` | Precipitation factor gamma distribution alpha | -| `calib.MCMC_params.kp_gamma_beta` | `float` | Precipitation factor gamma distribution beta | -| `calib.MCMC_params.kp_lognorm_mu` | `float` | Precipitation factor lognormal distribution mean | -| `calib.MCMC_params.kp_lognorm_tau` | `float` | Precipitation factor lognormal distribution tau | -| `calib.MCMC_params.kp_mu` | `float` | Precipitation factor normal distribution mean | -| `calib.MCMC_params.kp_sigma` | `float` | Precipitation factor normal distribution standard deviation | -| `calib.MCMC_params.kp_bndlow` | `float` | Lower bound for precipitation factor | -| `calib.MCMC_params.kp_bndhigh` | `float` | Upper bound for precipitation factor | +| `calib.MCMC_params.tbias_mu` | `float` | Temperature bias mean [K] | +| `calib.MCMC_params.tbias_sigma` | `float` | Temperature bias standard deviation [K] | +| `calib.MCMC_params.tbias_bndlow` | `float` | Lower bound for temperature bias [K] | +| `calib.MCMC_params.tbias_bndhigh` | `float` | Upper bound for temperature bias [K] | +| `calib.MCMC_params.kp_gamma_alpha` | `float` | Precipitation factor gamma distribution alpha [-] | +| `calib.MCMC_params.kp_gamma_beta` | `float` | Precipitation factor gamma distribution beta [-] | +| `calib.MCMC_params.kp_lognorm_mu` | `float` | Precipitation factor lognormal distribution mean [-] | +| `calib.MCMC_params.kp_lognorm_tau` | `float` | Precipitation factor lognormal distribution tau [-] | +| `calib.MCMC_params.kp_mu` | `float` | Precipitation factor normal distribution mean [-] | +| `calib.MCMC_params.kp_sigma` | `float` | Precipitation factor normal distribution standard deviation [-]| +| `calib.MCMC_params.kp_bndlow` | `float` | Lower bound for precipitation factor [-] | +| `calib.MCMC_params.kp_bndhigh` | `float` | Upper bound for precipitation factor [-] | ### Calibration Datasets | Variable | Type | Comment/Note | | :--- | :--- | :--- | -| `calib.data.massbalance.hugonnet2021_relpath` | `string` | Path to Hugonnet (2021) mass balance data | +| `calib.data.massbalance.hugonnet2021_relpath` | `string` | Path to Hugonnet (2021) mass balance data (relative to `root`) | | `calib.data.massbalance.hugonnet2021_fn` | `string` | Filename for Hugonnet (2021) mass balance data | | `calib.data.massbalance.hugonnet2021_facorrected_fn` | `string` | Filename for corrected Hugonnet (2021) mass balance data | -| `calib.data.oib.oib_relpath` | `string` | Path to OIB lidar data | -| `calib.data.oib.oib_rebin` | `integer` | Elevation rebinning in meters | -| `calib.data.oib.oib_filter_pctl` | `float` | Pixel count percentile filter | -| `calib.data.frontalablation.frontalablation_relpath` | `string` | Path to frontal ablation data | +| `calib.data.oib.oib_relpath` | `string` | Path to OIB lidar surface elevation data (relative to `root`) | +| `calib.data.oib.oib_rebin` | `integer` | Elevation rebinning interval [meters] | +| `calib.data.oib.oib_filter_pctl` | `float` | Pixel count percentile filter (filters out any data from each survey where the pixel count for a given bin is below the specified threshold) | +| `calib.data.frontalablation.frontalablation_relpath` | `string` | Path to frontal ablation data (relative to `root`) | | `calib.data.frontalablation.frontalablation_cal_fn` | `string` | Filename for frontal ablation calibration data | -| `calib.data.icethickness.h_consensus_relpath` | `string` | Path to ice thickness data | +| `calib.data.icethickness.h_ref_relpath` | `string` | Path to reference ice thickness data (relative to `root`) | ### Ice Thickness Calibration @@ -234,29 +231,41 @@ The hindcast option will flip the climate data array so 1960-2000 would run 2000 | Variable | Type | Comment/Note | | :--- | :--- | :--- | -| `sim.option_dynamics` | `null` or `string` | Glacier dynamics scheme (`OGGM`, `MassRedistributionCurves`, `null`) | -| `sim.option_bias_adjustment` | `integer` | Bias adjustment option (`0`, `1`, `2`, `3`) | -| `sim.nsims` | `integer` | Number of simulations | +| `sim.option_dynamics` | `null` or `string` | Glacier dynamics scheme (`'OGGM'`, `'MassRedistributionCurves'`, `'null'`) | +| `sim.option_bias_adjustment` | `integer` | Bias adjustment option (`0`: no adjustment, `1`: precipitation factor and temp bias scaling (builds off HH2015 methods), `2`: HH2015 methods, `3`: quantile delta mapping) | +| `sim.nsims` | `integer` | Number of simulations to run. Default is 1. More than 1 simulation can be run for glaciers that have been calibrated using the `'MCMC'` option. (Running 1 simulation with `sim.option_dynamics == 'MCMC` will take the **median** parameter values) | ### Output Options | Variable | Type | Comment/Note | | :--- | :--- | :--- | -| `sim.out.sim_stats` | `list` | Output statistics (`mean`, `std`, `median`, etc.) | -| `sim.out.export_all_simiters` | `boolean` | Export all simulation results | -| `sim.out.export_extra_vars` | `boolean` | Export extra variables (temp, prec, melt, etc.) | -| `sim.out.export_binned_data` | `boolean` | Export binned ice thickness | +| `sim.out.sim_stats` | `list` of `strings` | Output statistics (`'mean'`, `'std'`, `'median'`, etc.) | +| `sim.out.export_all_simiters` | `boolean` | Export all simulation results (for `sim.nsims` > 1), otherwise calculate statistics specified by `sim.out.sim_stats` | +| `sim.out.export_extra_vars` | `boolean` | Export extra variables (temperature, precipitation, melt, etc.) | +| `sim.out.export_binned_data` | `boolean` | Export binned ice data | | `sim.out.export_binned_components` | `boolean` | Export mass balance components | +### Dynamics +| Variable | Type | Comment/Note | +| :--- | :--- | :--- | +| `sim.oggm_dynamics.cfl_number` | `float` | Time step threshold [seconds] | +| `sim.oggm_dynamics.cfl_number_calving` | `float` | Time step threshold for marine-terimating glaciers [seconds] | +| `sim.oggm_dynamics.use_reg_glena` | `boolean` | Use regionally averaged glen_a creep paremeter | +| `sim.oggm_dynamics.glena_reg_relpath` | `string` | Path to comma-separated value file containing regionally averaged `glen_a` and `fs` parameter values (relative to `root`) | +| `sim.oggm_dynamics.fs` | `float` | Oerlemans sliding parameter | +| `sim.oggm_dynamics.glen_a_multiplier` | `float` | glen_a multiplier if not using regionally calibrated glen_a | +| `sim.icethickness_advancethreshold` | `float` | advancing glacier ice thickness change threshold (if not running OGGM dynamics) | +| `sim.terminus_percentage` | `float` | glacier elevation range considered terminus, used to size advancing new bin (if not running OGGM dynamics) [%] | + ### Model Parameters | Variable | Type | Comment/Note | | :--- | :--- | :--- | -| `params.use_constant_lapserate` | `boolean` | Use constant or variable lapse rate | -| `params.kp` | `float` | Precipitation factor | -| `params.tbias` | `float` | Temperature bias ($^\circ$C) | -| `params.ddfsnow` | `float` | Snow degree-day factor (m w.e. d$^{-1}$ $^\circ$C$^{-1}$) | -| `params.ddfsnow_iceratio` | `float` | Ratio of snow to ice degree-day factor | +| `sim.params.use_constant_lapserate` | `boolean` | Use constant or variable lapse rate | +| `sim.params.kp` | `float` | Precipitation factor [-] | +| `sim.params.tbias` | `float` | Temperature bias [K] | +| `sim.params.ddfsnow` | `float` | Snow degree-day factor [mm w.e. d$^{-1}$ K$^{-1}$] | +| `sim.params.ddfsnow_iceratio` | `float` | Ratio of snow to ice degree-day factor [-] | --- (input_mb_model_options_target)= @@ -264,27 +273,25 @@ The hindcast option will flip the climate data array so 1960-2000 would run 2000 | Variable | Type | Comment/Note | | :--- | :--- | :--- | -| `mb.option_surfacetype_initial` | `integer` | 1: median elevation (default), 2: mean elevation | +| `mb.option_surfacetype_initial` | `integer` | `1`: median elevation (default), `2`: mean elevation | | `mb.include_firn` | `boolean` | `true`: firn included, `false`: firn is modeled as snow | | `mb.include_debris` | `boolean` | `true`: account for debris with melt factors, `false`: do not account for debris | -| `mb.debris_relpath` | `string` | Path to debris dataset | +| `mb.debris_relpath` | `string` | Path to debris dataset (relative to `root`) | | `mb.option_elev_ref_downscale` | `string` | `'Zmed'`, `'Zmax'`, or `'Zmin'` for median, max, or min glacier elevations | -| `mb.option_temp2bins` | `integer` | `1`: adjust temp from GCM to glacier bins using lapse rates | -| `mb.option_adjusttemp_surfelev` | `integer` | `1`: adjust temps based on surface elevation changes, `0`: no adjustment | -| `mb.option_prec2bins` | `integer` | `1`: adjust precipitation from GCM to glacier bins | -| `mb.option_preclimit` | `integer` | `1`: limit uppermost 25% using an exponential function | +| `mb.option_adjusttemp_surfelev` | `integer` | `1`: adjust temperatures based on surface elevation changes, `0`: no adjustment | +| `mb.option_preclimit` | `integer` | `1`: limit uppermost 25% using an exponential function (currently the only scheme in place - any value other than 1 will result in no precipitation limit to be applied) | | `mb.option_accumulation` | `integer` | `1`: single threshold, `2`: threshold ± 1$^\circ$C using linear interpolation | -| `mb.option_ablation` | `integer` | `1`: monthly temp, `2`: superimposed daily temps enabling melt near 0$^\circ$C | +| `mb.option_ablation` | `integer` | `1`: monthly temperatures, `2`: superimpose daily temperatures to allow melt during months with negative mean temperatures but occasional positive daily values | | `mb.option_ddf_firn` | `integer` | `0`: ddf_firn = ddf_snow, `1`: ddf_firn = mean of ddf_snow and ddf_ice | | `mb.option_refreezing` | `string` | `'Woodward'`: annual air temp (Woodward et al. 1997), `'HH2015'`: heat conduction (Huss & Hock 2015) | -| `mb.Woodard_rf_opts.rf_month` | `integer` | Refreezing month | -| `mb.HH2015_rf_opts.rf_layers` | `integer` | Number of refreezing layers (default: 8) | -| `mb.HH2015_rf_opts.rf_dz` | `integer` | Layer thickness (m) | +| `mb.Woodard_rf_opts.rf_month` | `integer` | Month number in which refreezing is reset to zero each year (`1`-`12`, dafult is 10) | +| `mb.HH2015_rf_opts.rf_layers` | `integer` | Number of refreezing layers (default is 8) | +| `mb.HH2015_rf_opts.rf_dz` | `integer` | Layer thickness [m] | | `mb.HH2015_rf_opts.rf_dsc` | `integer` | Number of time steps for numerical stability | -| `mb.HH2015_rf_opts.rf_meltcrit` | `float` | Critical melt amount (m w.e.) to initialize refreezing | +| `mb.HH2015_rf_opts.rf_meltcrit` | `float` | Critical melt amount to initialize refreezing [m w.e.] | | `mb.HH2015_rf_opts.pp` | `float` | Additional refreeze water fraction for bare ice | -| `mb.HH2015_rf_opts.rf_dens_top` | `integer` | Snow density at surface (kg m$^{-3}$) | -| `mb.HH2015_rf_opts.rf_dens_bot` | `integer` | Snow density at bottom refreezing layer (kg m$^{-3}$) | +| `mb.HH2015_rf_opts.rf_dens_top` | `integer` | Snow density at surface [kg m$^{-3}$] | +| `mb.HH2015_rf_opts.rf_dens_bot` | `integer` | Snow density at bottom refreezing layer [kg m$^{-3}$] | | `mb.HH2015_rf_opts.option_rf_limit_meltsnow` | `integer` | Refreezing limit option | --- @@ -293,7 +300,7 @@ The hindcast option will flip the climate data array so 1960-2000 would run 2000 | Variable | Type | Comment/Note | | :--- | :--- | :--- | -| `rgi.rgi_relpath` | `string` | Filepath for RGI files | +| `rgi.rgi_relpath` | `string` | Filepath for RGI files (relative to `root`) | | `rgi.rgi_lat_colname` | `string` | Name of latitude column | | `rgi.rgi_lon_colname` | `string` | Name of longitude column (360-based) | | `rgi.elev_colname` | `string` | Name of elevation column | @@ -308,13 +315,12 @@ The hindcast option will flip the climate data array so 1960-2000 would run 2000 | Variable | Type | Comment/Note | | :--- | :--- | :--- | -| `time.option_leapyear` | `integer` | `1`: include leap years, `0`: exclude leap years (February always 28 days) | +| `time.option_leapyear` | `integer` | `1`: include 29FEB in leap years, `0`: exclude 29FEB during leap years | | `time.startmonthday` | `string` | Start date (MM-DD), used with custom calendars | | `time.endmonthday` | `string` | End date (MM-DD), used with custom calendars | | `time.wateryear_month_start` | `integer` | Water year start month | | `time.winter_month_start` | `integer` | First month of winter | | `time.summer_month_start` | `integer` | First month of summer | -| `time.option_dates` | `integer` | `1`: Use dates from date table (1st of each month), `2`: Use dates from climate data | | `time.timestep` | `string` | Time step (currently only `'monthly'` is supported) | --- @@ -323,20 +329,14 @@ The hindcast option will flip the climate data array so 1960-2000 would run 2000 | Variable | Type | Comment/Note | | :--- | :--- | :--- | -| `constants.density_ice` | `integer` | Density of ice (kg m$^{-3}$) | -| `constants.density_water` | `integer` | Density of water (kg m$^{-3}$) | -| `constants.area_ocean` | `float` | Ocean surface area (m$^{2}$) | -| `constants.k_ice` | `float` | Thermal conductivity of ice (J s$^{-1}$ K$^{-1}$ m$^{-1}$) | -| `constants.k_air` | `float` | Thermal conductivity of air (J s$^{-1}$ K$^{-1}$ m$^{-1}$) | -| `constants.ch_ice` | `integer` | Volumetric heat capacity of ice (J K$^{-1}$ m$^{-3}$) | -| `constants.ch_air` | `integer` | Volumetric heat capacity of air (J K$^{-1}$ m$^{-3}$) | -| `constants.Lh_rf` | `integer` | Latent heat of fusion (J kg$^{-1}$) | +| `constants.density_ice` | `integer` | Density of ice [900 kg m$^{-3}$] | +| `constants.density_water` | `integer` | Density of water [1000 kg m$^{-3}$] | +| `constants.k_ice` | `float` | Thermal conductivity of ice [2.330 J s$^{-1}$ K$^{-1}$ m$^{-1}$] | +| `constants.k_air` | `float` | Thermal conductivity of air [0.023 J s$^{-1}$ K$^{-1}$ m$^{-1}$] | +| `constants.ch_ice` | `integer` | Volumetric heat capacity of ice [1890000 J K$^{-1}$ m$^{-3}$] | +| `constants.ch_air` | `integer` | Volumetric heat capacity of air [1297 J K$^{-1}$ m$^{-3}$] | +| `constants.Lh_rf` | `integer` | Latent heat of fusion [333550 J kg$^{-1}$] | | `constants.tolerance` | `float` | Model tolerance (used to remove rounding errors) | -| `constants.gravity` | `float` | Gravity (m s$^{-2}$) | -| `constants.pressure_std` | `integer` | Standard atmospheric pressure (Pa) | -| `constants.temp_std` | `float` | Standard temperature (K) | -| `constants.R_gas` | `float` | Universal gas constant (J mol$^{-1}$ K$^{-1}$) | -| `constants.molarmass_air` | `float` | Molar mass of Earth's air (kg/mol) | --- (input_debugging_options)= diff --git a/docs/contributing.md b/docs/contributing.md index 06845764..36dfce8f 100644 --- a/docs/contributing.md +++ b/docs/contributing.md @@ -12,7 +12,7 @@ Next, clone PyGEM. This will place the code at your current directory, so you ma ``` git clone https://github.com/PyGEM-Community/PyGEM.git ``` -If you opted to create your own fork, clone using appropriate repo URL: `git clone https://github.com/YOUR-USERNAME/PyGEM.git` +If you opted to create your own fork, clone using appropriate repo URL: `git clone https://github.com//PyGEM.git` Navigate to root project directory: ``` diff --git a/docs/index.rst b/docs/index.rst index 7fb88910..34c90935 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -48,8 +48,9 @@ with respect to modeling glacier dynamics and ice thickness inversions. :caption: Getting Started: install_pygem - test_pygem + configuration scripts_overview + test_pygem model_output .. toctree:: @@ -57,8 +58,12 @@ with respect to modeling glacier dynamics and ice thickness inversions. :caption: Contributing: contributing - citing +.. toctree:: + :maxdepth: 1 + :caption: Citing: + + citing .. Indices and tables .. ================== diff --git a/docs/install_pygem.md b/docs/install_pygem.md index f1e28dbd..531a5b9a 100644 --- a/docs/install_pygem.md +++ b/docs/install_pygem.md @@ -1,6 +1,8 @@ (install_pygem_target)= # Installation -The Python Glacier Evolution Model has been packaged using Poetry and is hosted on the Python Package Index ([PyPI](https://pypi.org/project/pygem/)), to ensure that all dependencies install seamlessly. It is recommended that users create a [conda](https://docs.conda.io/projects/conda/en/latest/user-guide/index.html) environment from which to install the model dependencies and core code. If you do not yet have conda installed, see [conda's documentation](https://docs.conda.io/projects/conda/en/latest/user-guide/install) for instructions. +PyGEM has been tested successfully on Linux and macOS systems. For Windows users (Windows 10/11), we recommend installing the [Windows Subsystem for Linux](https://docs.microsoft.com/en-us/windows/wsl) (WSL), and then installing and running PyGEM from there. + +PyGEM has been packaged using [Poetry](https://python-poetry.org/) and is hosted on the Python Package Index ([PyPI](https://pypi.org/project/pygem/)), to ensure that all dependencies install seamlessly. It is recommended that users create a [conda](https://docs.conda.io/projects/conda/en/latest/user-guide/index.html) environment from which to install the model dependencies and core code. If you do not yet have conda installed, see [conda's documentation](https://docs.conda.io/projects/conda/en/latest/user-guide/install) for instructions. Next, choose your preferred PyGEM installation option:
- [**stable**](stable_install_target): this is the latest version that has been officially released to PyPI, with a fixed version number (e.g. v1.0.1). It is intended for general use. diff --git a/docs/mb_ablation.md b/docs/mb_ablation.md index afce142b..568a9b79 100644 --- a/docs/mb_ablation.md +++ b/docs/mb_ablation.md @@ -3,13 +3,13 @@ There are currently two model options for ablation. Both model options use a deg ### Option 1: monthly temperatures The first calculates ablation ($a$) using the mean monthly temperature: -$$a=f_{snow/ice/firn/debris} \cdot T_{m}^{+} \cdot n$$ +$a=f_{snow/ice/firn/debris} \cdot T_{m}^{+} \cdot n$ -where $f$ is the degree-day factor of snow, ice, firn, or debris (m w.e. d-1 °C$^{-1}$), $T_{m}^{+}$ is the positive monthly mean temperature (°C), and $n$ is the number of days per month. +where $f$ is the degree-day factor of snow, ice, firn, or debris (m w.e. d$^{-1}$ °C$^{-1}$), $T_{m}^{+}$ is the positive monthly mean temperature (°C), and $n$ is the number of days per month. ### Option 2: monthly temperatures with daily variance The second option incorporates the daily variance associated with the temperature for each month according to [Huss and Hock (2015)](https://www.frontiersin.org/articles/10.3389/feart.2015.00054/full): -$$a=f_{snow/ice/firn/debris} \cdot \sum_{i=1}^{ndays} T_{d,i}^{+} $$ +$a=f_{snow/ice/firn/debris} \cdot \sum_{i=1}^{ndays} T_{d,i}^{+} $ where $T_{d}$ is the daily positive mean air temperature and is estimated by superimposing random variability from the standard deviation of the daily temperature for each month. @@ -17,6 +17,6 @@ The degree-day factors for snow, ice, firn, and debris depend on the surface typ ### Temperature at elevation bins Temperature for each elevation bin ($T_{bin}$) is determined by selecting the temperature from the gridded climate data ($T_{gcm}$) based on the nearest neighbor, which is then downscaled to the elevation bins on the glacier according to: -$$T_{bin} = T_{gcm} + lr_{gcm} \cdot (z_{ref} - z_{gcm}) + lr_{glac} \cdot (z_{bin} - z_{ref}) + T_{bias}$$ +$T_{bin} = T_{gcm} + lr_{gcm} \cdot (z_{ref} - z_{gcm}) + lr_{glac} \cdot (z_{bin} - z_{ref}) + T_{bias}$ where $lr_{gcm}$ and $lr_{glac}$ are lapse rates (°C m-1) associated with downscaling the climate data to the glacier and then over the glacier elevation bins, respectively; $z_{ref}$, $z_{gcm}$, and $z_{bin}$ are the elevations from the glacier’s reference point (median or mean elevation), the climate data, and the elevation bin, respectively; and $T_{bias}$ is the temperature bias. The temperature bias is one of three model parameters that is calibrated and serves to account for any biases resulting from the use of coarse climate data that is unable to capture local topographic variations. By default, the $lr_{gcm}$ and $lr_{glac}$ are assumed to be equal. \ No newline at end of file diff --git a/docs/mb_accumulation.md b/docs/mb_accumulation.md index b329569a..cb09d3ed 100644 --- a/docs/mb_accumulation.md +++ b/docs/mb_accumulation.md @@ -3,7 +3,7 @@ Accumulation ($c$) is calculated for each elevation bin as a function of the pre ### Option 1: Threshold +/- 1$^{\circ}$C The first (default) option is to estimate the ratio of liquid and solid precipitation based on the air temperature: -$$c = \delta \cdot P_{bin}$$ +$c = \delta \cdot P_{bin}$ where $\delta=1$; if $T_{bin} \leq T_{snow}-1$ @@ -18,9 +18,9 @@ The alternative option is to classify precipitation as snow or rain based on a s ### Precipitation at elevation bins Precipitation at each elevation bin of the glacier ($P_{bin}$) is determined by selecting the precipitation from the gridded climate data ($P_{gcm}$) based on the nearest neighbor, which is then downscaled to the elevation bins on the glacier: -$$P_{bin} = P_{GCM} \cdot k_{p} \cdot (1 + d_{prec} \cdot (z_{bin} - z_{ref}))$$ +$P_{bin} = P_{GCM} \cdot k_{p} \cdot (1 + d_{prec} \cdot (z_{bin} - z_{ref}))$
where $k_{p}$ is the precipitation factor and $d_{prec}$ is the precipitation gradient. The precipitation factor is a model parameter that is used to adjust from the climate data to the glacier, which could be caused by local topographic effects due to differences in elevation, rain shadow effects, etc. The precipitation gradient is another model parameter, which is used to redistribute the precipitation along the glacier and can be thought of as a precipitation lapse rate. Typical values for the precipitation gradient vary from 0.01 – 0.025% m$^{-1}$([Huss and Hock, 2015](https://www.frontiersin.org/articles/10.3389/feart.2015.00054/full) who cited [WGMS, 2012](https://wgms.ch/products_fog/)). The default assumes a precipitation gradient of 0.01% m$^{-1}$ to reduce the number of model parameters. Additionally, for glaciers with high relief (> 1000 m), the precipitation in the uppermost 25% of the glacier’s elevation is reduced using an exponential function ([Huss and Hock, 2015](https://www.frontiersin.org/articles/10.3389/feart.2015.00054/full)): -$$P_{bin,exp} = P_{bin} \cdot exp(\frac{z_{bin} - z_{75\%}}{z_{max} - z_{75\%}}) $$ +$P_{bin,exp} = P_{bin} \cdot exp(\frac{z_{bin} - z_{75\%}}{z_{max} - z_{75\%}}) $ where $P_{bin,exp}$ is the adjusted precipitation, and $z_{max}$ and $z_{75\%}$ are the elevation of the glacier maximum and the glacier’s 75th percentile elevation, respectively. The adjusted precipitation cannot be lower than 87.5% of the maximum precipitation on the glacier. This adjustment accounts for the reduced air moisture and increased wind erosion at higher elevations ([Benn and Lehmkuhl, 2000](https://risweb.st-andrews.ac.uk/portal/en/researchoutput/mass-balance-and-equilibriumline-altitudes-of-glaciers-in-highmountain-environments(080f17fc-33dd-4805-bc97-a5aaa018a457)/export.html)). \ No newline at end of file diff --git a/docs/mb_refreezing.md b/docs/mb_refreezing.md index 74e3dfb9..67ba3c82 100644 --- a/docs/mb_refreezing.md +++ b/docs/mb_refreezing.md @@ -3,7 +3,7 @@ There are two model options for computing refreezing. The default option estima ### Option 1: Mean annual air temperature For the default option, refreezing (R) is calculated for each elevation bin as a function of its weighted annual mean air temperature (Ta) following [Woodward et al. (1997)](https://www.cambridge.org/core/journals/annals-of-glaciology/article/influence-of-superimposedice-formation-on-the-sensitivity-of-glacier-mass-balance-to-climate-change/84DFA08E9CC8F28BE0729F1EBF4DA4E1)): -$$R = -0.0069 \cdot T_{a} + 0.000096$$ +$R = -0.0069 \cdot T_{a} + 0.000096$ The weighted annual mean accounts for the number of days in each month. Refreezing cannot be negative. The model assumes that refreezing occurs in the snowpack as opposed to being superimposed ice, so in the ablation zone the refreezing cannot exceed the snow depth. Since this option estimates annual refreezing, the equation above provides the maximum amount of potential refreezing. Each year, the potential refreezing is reset in October as this is the transition season from predominantly summer melt to winter accumulation. Therefore, it is possible that accumulated snow may melt and refreeze diurnally. The model replicates this physical behavior by first determining the amount of snow that has melted in that month. If the refreezing potential is greater than zero, the model assumes the snow melts and refreezes in the snow pack. Refreezing cannot exceed the amount of snow melt in any given month. The amount of refreezing is then subtracted from the potential refreezing until the potential refreezing is fully depleted or reset. Once the snow and refreezing completely melts, the model can melt the underlying ice/firn/snow. This model is used as the default because it is computationally cheap compared to the alternative. ### Option 2: Heat conduction diff --git a/docs/model_output.md b/docs/model_output.md index da20a91a..b24f2be3 100644 --- a/docs/model_output.md +++ b/docs/model_output.md @@ -1,6 +1,8 @@ (model_output_overview_target)= # Model Output -The model outputs a variety of data including monthly mass balance and its components (accumulation, melt, refreezing, frontal ablation, glacier runoff), and annual mass, mass below sea level, and area. Results are written as a netcdf file (.nc) for each glacier. If multiple simulations are performed (e.g., for Monte Carlo simulations), then statistics related to the median and median absolute deviation are output for each parameter. In addition to this standard glacier-wide output, binned output is also available, which include the bins surface elevation, volume, thickness, and climatic mass balance annually. +The model outputs a variety of data including monthly mass balance and its components (accumulation, melt, refreezing, frontal ablation, glacier runoff), and annual mass, mass below sea level, and area. Results are written as a netcdf file (.nc) for each glacier. If multiple simulations are performed (e.g., for Monte Carlo simulations), then statistics related to the median and median absolute deviation are output for each parameter. + +In addition to this standard glacier-wide output, binned outputs are also available (by setting `sim["out"]["export_binned_data"]` to `true` in the [PyGEM configuration file](contributing_pygem_target)), which include each bins initial surface elevation, annual area, annual mass, annual thickness, annual climatic mass balance, and monthly climatic mass balance. Monthly climatic mass balance components can also be stored by setting `sim["out"]["export_binned_components"]` to `true`. ## Post-processing Data PyGEM simulations are output for each glacier individually. For most analyses, it is useful to aggregate or merge analyses to a regional or global scale. PyGEM's *postproc.compile_simulations.py* is designed to do just so. diff --git a/docs/model_structure.md b/docs/model_structure.md index f0df984d..086f9326 100644 --- a/docs/model_structure.md +++ b/docs/model_structure.md @@ -15,7 +15,7 @@ Currently, the model does not have a “required” set of directories. The rela ~/pygem_data/
├ [DEMs](input_mb_data_target): geodetic mass balance data derived from DEM differencing that is used for calibration.
-├ [IceThickness_Farinotti](input_thickness_data_target): consensus ice thickness estimates (Farinotti et al. 2019). The directory is optional unless you want to match the consensus ice thickness estimates
+├ [IceThickness_Farinotti](input_thickness_data_target): reference ice thickness estimates (Farinotti et al. 2019). The directory is optional unless you want to match the reference ice thickness estimates
├ [oggm_gdirs](input_glacier_data_target): glacier directories downloaded from OGGM. This directory will be automatically generated by the pre-processing steps from OGGM.
├ Output: model output data
├ [RGI](input_glacier_inventory_target): Randolph Glacier Inventory information
@@ -33,7 +33,7 @@ The model code itself is heavily commented with the hope that the code is easy t * [Pre-process data](preprocessing_target) (optional if including more data) * [Set up configuration file](config_workflow_target) * [Calibrate frontal ablation parameter](workflow_cal_frontalablation_target) (optional for marine-terimating glaciers) -* [Calibrate mass balance parameters](workflow_cal_prms_target) +* [Calibrate climatic mass balance parameters](workflow_cal_prms_target) * [Calibrate ice viscosity parameter](workflow_cal_glena_target) * [Run model simulation](workflow_sim_target) * [Post-process output](workflow_post_target) @@ -57,7 +57,7 @@ The only critical component of the configuration file which must be set prior to For more details, see the [configuration file overview](pygem_config_overview_target). ```{note} -The first time you process glacier directories (e.g., happens automatically with the run_calibration.py), you want to set overwrite_gdirs=True in *~/PyGEM/config.yaml*. This will add the mass balance and consensus ice thickness data to the glacier directories. To avoid redownloading/reprocessing the glacier directories every time, after performed once set overwrite_gdirs=False. +The first time you process glacier directories (e.g., happens automatically with the run_calibration.py), you want to set overwrite_gdirs=True in *~/PyGEM/config.yaml*. This will add the mass balance and reference ice thickness data to the glacier directories. To avoid redownloading/reprocessing the glacier directories every time, after performed once set overwrite_gdirs=False. ``` (workflow_cal_prms_target)= @@ -94,7 +94,7 @@ Circularity issues exist in calibrating the frontal ablation parameter as the ma (workflow_cal_glena_target)= ### Calibrate ice viscosity model parameter -The ice viscosity ("Glen A") model parameter is calibrated such that the ice volume estimated using the calibrated mass balance gradients are consistent with the consensus ice volume estimates ([Farinotti et al. (2019)](https://www.nature.com/articles/s41561-019-0300-3)) for each RGI region. This is done by running the following: +The ice viscosity ("Glen A") model parameter is calibrated such that the ice volume estimated using the calibrated mass balance gradients are consistent with the reference ice volume estimates ([Farinotti et al. (2019)](https://www.nature.com/articles/s41561-019-0300-3)) for each RGI region. This is done by running the following: ``` run_calibration_reg_glena ``` diff --git a/docs/publications.md b/docs/publications.md index 582682b4..d3c2e57c 100644 --- a/docs/publications.md +++ b/docs/publications.md @@ -1,5 +1,10 @@ (publications_target)= -# Select Publications +# Publications +## Publications describing PyGEM * Rounce, D.R., Hock, R., Maussion, F., Hugonnet, R., Kochtitzky, W., Huss, M., Berthier, E., Brinkerhoff, D., Compagno, L., Copland, L., Farinotti, D., Menounos, B., and McNabb, R.W. (2023). “[Global glacier change in the 21st century: Every increase in temperature matters](https://www.science.org/doi/10.1126/science.abo1324)”, Science, 379(6627), pp. 78-83, doi:10.1126/science.abo1324 * Rounce, D.R., Hock, R., and Shean, D.E. (2020). “[Glacier mass change in High Mountain Asia through 2100 using the open-source Python Glacier Evolution Model (PyGEM)](https://www.frontiersin.org/articles/10.3389/feart.2019.00331/full)”, Frontiers in Earth Science, 7(331), pp. 1-20, doi:10.3389/feart.2019.00331 -* Rounce, D.R., Khurana, T., Short, M.B., Hock, R., Shean, D.E., and Brinkherhoff, D.J. (2020). “[Quantifying parameter uncertainty in a large-scale glacier evolution model using Bayesian inference – Application to High Mountain Asia](https://www.cambridge.org/core/journals/journal-of-glaciology/article/quantifying-parameter-uncertainty-in-a-largescale-glacier-evolution-model-using-bayesian-inference-application-to-high-mountain-asia/61D8956E9A6C27CC1A5AEBFCDADC0432)”, Journal of Glaciology, 66(256), pp. 175-187, doi:10.1017/jog.2019.91 \ No newline at end of file +* Rounce, D.R., Khurana, T., Short, M.B., Hock, R., Shean, D.E., and Brinkherhoff, D.J. (2020). “[Quantifying parameter uncertainty in a large-scale glacier evolution model using Bayesian inference – Application to High Mountain Asia](https://www.cambridge.org/core/journals/journal-of-glaciology/article/quantifying-parameter-uncertainty-in-a-largescale-glacier-evolution-model-using-bayesian-inference-application-to-high-mountain-asia/61D8956E9A6C27CC1A5AEBFCDADC0432)”, Journal of Glaciology, 66(256), pp. 175-187, doi:10.1017/jog.2019.91 + +## Publications using PyGEM +* Yang, R, Hock, R., Rounce, D., Kang, S. (2024). "[Regional-scale response of glacier speed to seasonal runoff variations on the Kenai Peninsula, Alaska.](https://doi.org/10.1029/2025GL115248)", Geophysical Research Letters. 52, e2025GL115248, doi:10.1029/2025GL115248. +* Kochtitzky W., Copland, L., Van Wychen, W., Hock, R., Rounce, D.R., Jiskoot, H., Scambos, T.A., Morlighem, M., King, M., Cha, L., Gould, L., Merrill,P-M., Glazovsky, A., Hugonnet, R., Strozzi, T., Noël, B., Navarro, F., Millan, R., Dowdeswell, J.A., Cook, A., Dalton, A., Khan, S., Jania, J. (2023). "[Progress towards globally complete frontal ablation estimates of marine-terminating glaciers](https://doi.org/10.1017/aog.2023.35)". Annals of Glaciology. doi:10.1017/aog.2023.35. \ No newline at end of file diff --git a/docs/scripts_overview.md b/docs/scripts_overview.md index 8f78bca6..06118aa3 100644 --- a/docs/scripts_overview.md +++ b/docs/scripts_overview.md @@ -12,14 +12,4 @@ run_calibration_frontalablation_overview run_calibration_overview run_calibration_reg_glena_overview run_simulation_overview -``` - - -```{toctree} ---- -caption: Other: -maxdepth: 1 ---- - -config.yaml ``` \ No newline at end of file diff --git a/pygem/bin/op/list_failed_simulations.py b/pygem/bin/op/list_failed_simulations.py index 7ed08a24..29e98b22 100644 --- a/pygem/bin/op/list_failed_simulations.py +++ b/pygem/bin/op/list_failed_simulations.py @@ -27,7 +27,16 @@ import pygem.pygem_modelsetup as modelsetup -def run(reg, simpath, gcm, scenario, calib_opt, bias_adj, gcm_startyear, gcm_endyear): +def run( + reg, + simpath, + gcm, + sim_climate_scenario, + calib_opt, + bias_adj, + sim_startyear, + sim_endyear, +): # define base directory base_dir = simpath + '/' + str(reg).zfill(2) + '/' @@ -43,19 +52,19 @@ def run(reg, simpath, gcm, scenario, calib_opt, bias_adj, gcm_startyear, gcm_end glacno_list_all = list(main_glac_rgi_all['rgino_str'].values) # get list of glacier simulation files - if scenario: - sim_dir = base_dir + gcm + '/' + scenario + '/stats/' + if sim_climate_scenario: + sim_dir = base_dir + gcm + '/' + sim_climate_scenario + '/stats/' else: sim_dir = base_dir + gcm + '/stats/' - # check if gcm has given scenario + # check if gcm has given sim_climate_scenario assert os.path.isdir(sim_dir), f'Error: simulation path not found, {sim_dir}' # instantiate list of galcnos that are not in sim_dir failed_glacnos = [] fps = glob.glob( - sim_dir + f'*_{calib_opt}_ba{bias_adj}_*_{gcm_startyear}_{gcm_endyear}_all.nc' + sim_dir + f'*_{calib_opt}_ba{bias_adj}_*_{sim_startyear}_{sim_endyear}_all.nc' ) # Glaciers with successful runs to process @@ -67,7 +76,7 @@ def run(reg, simpath, gcm, scenario, calib_opt, bias_adj, gcm_startyear, gcm_end main_glac_rgi_all.apply(lambda x: x.rgino_str in glacno_ran, axis=1) ] print( - f'{gcm} {str(scenario).replace("None", "")} glaciers successfully simulated:\n - {main_glac_rgi.shape[0]} of {main_glac_rgi_all.shape[0]} glaciers ({np.round(main_glac_rgi.shape[0] / main_glac_rgi_all.shape[0] * 100, 3)}%)' + f'{gcm} {str(sim_climate_scenario).replace("None", "")} glaciers successfully simulated:\n - {main_glac_rgi.shape[0]} of {main_glac_rgi_all.shape[0]} glaciers ({np.round(main_glac_rgi.shape[0] / main_glac_rgi_all.shape[0] * 100, 3)}%)' ) print( f' - {np.round(main_glac_rgi.Area.sum(), 0)} km2 of {np.round(main_glac_rgi_all.Area.sum(), 0)} km2 ({np.round(main_glac_rgi.Area.sum() / main_glac_rgi_all.Area.sum() * 100, 3)}%)' @@ -88,7 +97,7 @@ def run(reg, simpath, gcm, scenario, calib_opt, bias_adj, gcm_startyear, gcm_end def main(): # Set up CLI parser = argparse.ArgumentParser( - description="""description: script to check for failed PyGEM glacier simulations\n\nexample call: $python list_failed_simulations.py -rgi_region01=1 -gcm_name=CanESM5 -scenrio=ssp585 -outdir=/path/to/output/failed/glaciers/""", + description="""description: script to check for failed PyGEM glacier simulations\n\nexample call: $python list_failed_simulations.py -rgi_region01=1 -sim_climate_name=CanESM5 -scenrio=ssp585 -outdir=/path/to/output/failed/glaciers/""", formatter_class=argparse.RawTextHelpFormatter, ) requiredNamed = parser.add_argument_group('required named arguments') @@ -100,30 +109,30 @@ def main(): nargs='+', ) parser.add_argument( - '-gcm_name', + '-sim_climate_name', type=str, default=None, help='GCM name to compile results from (ex. ERA5 or CESM2)', ) parser.add_argument( - '-scenario', + '-sim_climate_scenario', action='store', type=str, default=None, - help='rcp or ssp scenario used for model run (ex. rcp26 or ssp585)', + help='rcp or ssp sim_climate_scenario used for model run (ex. rcp26 or ssp585)', ) parser.add_argument( - '-gcm_startyear', + '-sim_startyear', action='store', type=int, - default=pygem_prms['climate']['gcm_startyear'], + default=pygem_prms['climate']['sim_startyear'], help='start year for the model run', ) parser.add_argument( - '-gcm_endyear', + '-sim_endyear', action='store', type=int, - default=pygem_prms['climate']['gcm_endyear'], + default=pygem_prms['climate']['sim_endyear'], help='start year for the model run', ) parser.add_argument( @@ -152,13 +161,13 @@ def main(): args = parser.parse_args() region = args.rgi_region01 - scenario = args.scenario - gcm_name = args.gcm_name + sim_climate_scenario = args.sim_climate_scenario + sim_climate_name = args.sim_climate_name bias_adj = args.option_bias_adjustment simpath = pygem_prms['root'] + '/Output/simulations/' - if gcm_name in ['ERA5', 'ERA-Interim', 'COAWST']: - scenario = None + if sim_climate_name in ['ERA5', 'ERA-Interim', 'COAWST']: + sim_climate_scenario = None bias_adj = 0 if not isinstance(region, list): @@ -172,31 +181,31 @@ def main(): failed_glacs = run( reg, simpath, - args.gcm_name, - scenario, + args.sim_climate_name, + sim_climate_scenario, args.option_calibration, bias_adj, - args.gcm_startyear, - args.gcm_endyear, + args.sim_startyear, + args.sim_endyear, ) if len(failed_glacs) > 0: if args.outdir: fout = os.path.join( args.outdir, - f'R{str(reg).zfill(2)}_{args.gcm_name}_{scenario}_{args.gcm_startyear}_{args.gcm_endyear}_failed_rgiids.json', + f'R{str(reg).zfill(2)}_{args.sim_climate_name}_{sim_climate_scenario}_{args.sim_startyear}_{args.sim_endyear}_failed_rgiids.json', ).replace('None_', '') with open(fout, 'w') as f: json.dump(failed_glacs, f) print( - f'List of failed glaciers for {gcm_name} {str(scenario).replace("None", "")} exported to: {fout}' + f'List of failed glaciers for {sim_climate_name} {str(sim_climate_scenario).replace("None", "")} exported to: {fout}' ) if args.verbose: print( - f'Failed glaciers for RGI region R{str(reg).zfill(2)} {args.gcm_name} {str(scenario).replace("None", "")} {args.gcm_startyear}-{args.gcm_endyear}:' + f'Failed glaciers for RGI region R{str(reg).zfill(2)} {args.sim_climate_name} {str(sim_climate_scenario).replace("None", "")} {args.sim_startyear}-{args.sim_endyear}:' ) print(failed_glacs) else: print( - f'No glaciers failed from R{region}, for {gcm_name} {scenario.replace("None", "")}' + f'No glaciers failed from R{region}, for {sim_climate_name} {sim_climate_scenario.replace("None", "")}' ) diff --git a/pygem/bin/postproc/postproc_compile_simulations.py b/pygem/bin/postproc/postproc_compile_simulations.py index cbaa4c96..72a19247 100644 --- a/pygem/bin/postproc/postproc_compile_simulations.py +++ b/pygem/bin/postproc/postproc_compile_simulations.py @@ -63,11 +63,11 @@ def run(args): simpath, gcms, realizations, - scenario, + sim_climate_scenario, calibration, bias_adj, - gcm_startyear, - gcm_endyear, + sim_startyear, + sim_endyear, vars, ) = args print(f'RGI region {reg}') @@ -110,34 +110,41 @@ def run(args): ############################################################ ### get time values - should be the same across all sims ### + ### also ensure that specified GCM/sim_climate_scenario pair was run ### ############################################################ - if scenario: - # ensure scenario has been run for each gcm - for gcm in gcms: - if scenario not in os.listdir(base_dir + '/' + gcm): - # remove the gcm from our gcm list if the desired scenario is not contained + for gcm in gcms: + gcm_path = os.path.join(base_dir, gcm) + if not os.path.exists(gcm_path): + print(f'{gcm} not found, skipping') + gcms.remove(gcm) + + if sim_climate_scenario: + sim_climate_scenario_path = os.path.join(gcm_path, sim_climate_scenario) + if not os.path.exists(sim_climate_scenario_path): + print(f'{sim_climate_scenario} not found for {gcm}, skipping') gcms.remove(gcm) - print(f'scenario {scenario} not found for {gcm}, skipping') - fn = glob.glob( - base_dir - + gcm - + '/' - + scenario - + '/stats/' - + f'*{gcm}_{scenario}_{realizations[0]}_{calibration}_ba{bias_adj}_*_{gcm_startyear}_{gcm_endyear}_all.nc'.replace( - '__', '_' - ) - )[0] - else: - fn = glob.glob( - base_dir - + gcm - + '/stats/' - + f'*{gcm}_{calibration}_ba{bias_adj}_*_{gcm_startyear}_{gcm_endyear}_all.nc' - )[0] - nsets = fn.split('/')[-1].split('_')[-4] - - ds_glac = xr.open_dataset(fn) + # get first file + fp = glob.glob( + base_dir + + gcm + + '/' + + sim_climate_scenario + + '/stats/' + + f'*{gcm}_{sim_climate_scenario}_{realizations[0]}_{calibration}_ba{bias_adj}_*_{sim_startyear}_{sim_endyear}_all.nc'.replace( + '__', '_' + ) + )[0] + else: + fp = glob.glob( + base_dir + + gcm + + '/stats/' + + f'*{gcm}_{calibration}_ba{bias_adj}_*_{sim_startyear}_{sim_endyear}_all.nc' + )[0] + # get number of sets from file name + nsets = fp.split('/')[-1].split('_')[-4] + # open dataset and read time/year values + ds_glac = xr.open_dataset(fp) year_values = ds_glac.year.values time_values = ds_glac.time.values # check if desired vars are in ds @@ -155,8 +162,6 @@ def run(args): ### LEVEL I ### # loop through glacier batches of 1000 for nbatch, glacno_list in enumerate(glacno_list_batches): - print(f'Batch {nbatch}:') - # batch start timer loop_start = time.time() @@ -164,8 +169,6 @@ def run(args): batch_start = glacno_list[0].split('.')[1] batch_end = glacno_list[-1].split('.')[1] - print(nbatch, batch_start, batch_end) - # get all glacier info for glaciers in batch main_glac_rgi_batch = main_glac_rgi_all.loc[ main_glac_rgi_all.apply(lambda x: x.rgino_str in glacno_list, axis=1) @@ -191,19 +194,19 @@ def run(args): # for each batch, loop through GCM(s) and realization(s) for gcm in gcms: # get list of glacier simulation files - sim_dir = base_dir + gcm + '/' + scenario + '/stats/' + sim_dir = base_dir + gcm + '/' + sim_climate_scenario + '/stats/' ### LEVEL III ### for realization in realizations: print(f'GCM: {gcm} {realization}') fps = glob.glob( sim_dir - + f'*{gcm}_{scenario}_{realization}_{calibration}_ba{bias_adj}_{nsets}_{gcm_startyear}_{gcm_endyear}_all.nc'.replace( + + f'*{gcm}_{sim_climate_scenario}_{realization}_{calibration}_ba{bias_adj}_{nsets}_{sim_startyear}_{sim_endyear}_all.nc'.replace( '__', '_' - ) + ).replace('__', '_') ) - # during 0th batch, print the regional stats of glaciers and area successfully simulated for all regional glaciers for given gcm scenario + # during 0th batch, print the regional stats of glaciers and area successfully simulated for all regional glaciers for given gcm sim_climate_scenario if nbatch == 0: # Glaciers with successful runs to process glacno_ran = [x.split('/')[-1].split('_')[0] for x in fps] @@ -230,6 +233,7 @@ def run(args): reg_glac_gcm_melt_monthly = None reg_glac_gcm_refreeze_monthly = None reg_glac_gcm_frontalablation_monthly = None + reg_glac_gcm_snowline_monthly = None reg_glac_gcm_massbaltotal_monthly = None reg_glac_gcm_prec_monthly = None reg_glac_gcm_mass_monthly = None @@ -237,19 +241,21 @@ def run(args): # annual vars reg_glac_gcm_area_annual = None reg_glac_gcm_mass_annual = None + reg_glac_gcm_ELA_annual = None + reg_glac_gcm_mass_bsl_annual = None ### LEVEL IV ### # loop through each glacier in batch list for i, glacno in enumerate(glacno_list): # get glacier string and file name glacier_str = '{0:0.5f}'.format(float(glacno)) - glacno_fn = f'{sim_dir}/{glacier_str}_{gcm}_{scenario}_{realization}_{calibration}_ba{bias_adj}_{nsets}_{gcm_startyear}_{gcm_endyear}_all.nc'.replace( + glacno_fp = f'{sim_dir}/{glacier_str}_{gcm}_{sim_climate_scenario}_{realization}_{calibration}_ba{bias_adj}_{nsets}_{sim_startyear}_{sim_endyear}_all.nc'.replace( '__', '_' - ) + ).replace('__', '_') # try to load all glaciers in region try: # open netcdf file - ds_glac = xr.open_dataset(glacno_fn) + ds_glac = xr.open_dataset(glacno_fp) # get monthly vars glac_runoff_monthly = ds_glac.glac_runoff_monthly.values offglac_runoff_monthly = ds_glac.offglac_runoff_monthly.values @@ -261,6 +267,7 @@ def run(args): glac_frontalablation_monthly = ( ds_glac.glac_frontalablation_monthly.values ) + glac_snowline_monthly = ds_glac.glac_snowline_monthly.values glac_massbaltotal_monthly = ( ds_glac.glac_massbaltotal_monthly.values ) @@ -275,6 +282,9 @@ def run(args): glac_frontalablation_monthly = np.full( (1, len(time_values)), np.nan ) + glac_snowline_monthly = np.full( + (1, len(time_values)), np.nan + ) glac_massbaltotal_monthly = np.full( (1, len(time_values)), np.nan ) @@ -283,6 +293,8 @@ def run(args): # get annual vars glac_area_annual = ds_glac.glac_area_annual.values glac_mass_annual = ds_glac.glac_mass_annual.values + glac_ELA_annual = ds_glac.glac_ELA_annual.values + glac_mass_bsl_annual = ds_glac.glac_mass_bsl_annual.values # if glacier output DNE in sim output file, create empty nan arrays to keep record of missing glaciers except: @@ -295,6 +307,7 @@ def run(args): glac_frontalablation_monthly = np.full( (1, len(time_values)), np.nan ) + glac_snowline_monthly = np.full((1, len(time_values)), np.nan) glac_massbaltotal_monthly = np.full( (1, len(time_values)), np.nan ) @@ -303,6 +316,10 @@ def run(args): # annual vars glac_area_annual = np.full((1, year_values.shape[0]), np.nan) glac_mass_annual = np.full((1, year_values.shape[0]), np.nan) + glac_ELA_annual = np.full((1, year_values.shape[0]), np.nan) + glac_mass_bsl_annual = np.full( + (1, year_values.shape[0]), np.nan + ) # append each glacier output to master regional set of arrays if reg_glac_gcm_mass_annual is None: @@ -315,12 +332,15 @@ def run(args): reg_glac_gcm_frontalablation_monthly = ( glac_frontalablation_monthly ) + reg_glac_gcm_snowline_monthly = glac_snowline_monthly reg_glac_gcm_massbaltotal_monthly = glac_massbaltotal_monthly reg_glac_gcm_prec_monthly = glac_prec_monthly reg_glac_gcm_mass_monthly = glac_mass_monthly # annual vars reg_glac_gcm_area_annual = glac_area_annual reg_glac_gcm_mass_annual = glac_mass_annual + reg_glac_gcm_ELA_annual = glac_ELA_annual + reg_glac_gcm_mass_bsl_annual = glac_mass_bsl_annual # otherwise concatenate existing arrays else: # monthly vars @@ -348,6 +368,10 @@ def run(args): ), axis=0, ) + reg_glac_gcm_snowline_monthly = np.concatenate( + (reg_glac_gcm_snowline_monthly, glac_snowline_monthly), + axis=0, + ) reg_glac_gcm_massbaltotal_monthly = np.concatenate( ( reg_glac_gcm_massbaltotal_monthly, @@ -368,6 +392,12 @@ def run(args): reg_glac_gcm_mass_annual = np.concatenate( (reg_glac_gcm_mass_annual, glac_mass_annual), axis=0 ) + reg_glac_gcm_ELA_annual = np.concatenate( + (reg_glac_gcm_ELA_annual, glac_ELA_annual), axis=0 + ) + reg_glac_gcm_mass_bsl_annual = np.concatenate( + (reg_glac_gcm_mass_bsl_annual, glac_mass_bsl_annual), axis=0 + ) # aggregate gcms if reg_glac_allgcms_runoff_monthly is None: @@ -390,6 +420,9 @@ def run(args): reg_glac_allgcms_frontalablation_monthly = ( reg_glac_gcm_frontalablation_monthly[np.newaxis, :, :] ) + reg_glac_allgcms_snowline_monthly = reg_glac_gcm_snowline_monthly[ + np.newaxis, :, : + ] reg_glac_allgcms_massbaltotal_monthly = ( reg_glac_gcm_massbaltotal_monthly[np.newaxis, :, :] ) @@ -406,6 +439,12 @@ def run(args): reg_glac_allgcms_mass_annual = reg_glac_gcm_mass_annual[ np.newaxis, :, : ] + reg_glac_allgcms_ELA_annual = reg_glac_gcm_ELA_annual[ + np.newaxis, :, : + ] + reg_glac_allgcms_mass_bsl_annual = reg_glac_gcm_mass_bsl_annual[ + np.newaxis, :, : + ] else: # monthly vrs reg_glac_allgcms_runoff_monthly = np.concatenate( @@ -450,6 +489,13 @@ def run(args): ), axis=0, ) + reg_glac_allgcms_snowline_monthly = np.concatenate( + ( + reg_glac_allgcms_snowline_monthly, + reg_glac_gcm_snowline_monthly[np.newaxis, :, :], + ), + axis=0, + ) reg_glac_allgcms_massbaltotal_monthly = np.concatenate( ( reg_glac_allgcms_massbaltotal_monthly, @@ -486,6 +532,20 @@ def run(args): ), axis=0, ) + reg_glac_allgcms_ELA_annual = np.concatenate( + ( + reg_glac_allgcms_ELA_annual, + reg_glac_gcm_ELA_annual[np.newaxis, :, :], + ), + axis=0, + ) + reg_glac_allgcms_mass_bsl_annual = np.concatenate( + ( + reg_glac_allgcms_mass_bsl_annual, + reg_glac_gcm_mass_bsl_annual[np.newaxis, :, :], + ), + axis=0, + ) # ===== CREATE NETCDF FILES===== @@ -657,6 +717,29 @@ def run(args): ) ds.glac_frontalablation_monthly.attrs['grid_mapping'] = 'crs' + # glac_snowline_monthly + elif var == 'glac_snowline_monthly': + ds = xr.Dataset( + data_vars=dict( + glac_snowline_monthly=( + coord_order, + reg_glac_allgcms_snowline_monthly, + ), + crs=np.nan, + ), + coords=coords_dict, + attrs=attrs_dict, + ) + ds.glac_snowline_monthly.attrs['long_name'] = ( + 'transient snowline altitude above mean sea level' + ) + ds.glac_snowline_monthly.attrs['units'] = 'm' + ds.glac_snowline_monthly.attrs['temporal_resolution'] = 'monthly' + ds.glac_snowline_monthly.attrs['comment'] = ( + 'transient snowline is altitude separating snow from ice/firn' + ) + ds.glac_snowline_monthly.attrs['grid_mapping'] = 'crs' + # glac_massbaltotal_monthly elif var == 'glac_massbaltotal_monthly': ds = xr.Dataset( @@ -752,6 +835,49 @@ def run(args): ) ds.glac_mass_annual.attrs['grid_mapping'] = 'crs' + # glac_ELA_annual + elif var == 'glac_ELA_annual': + ds = xr.Dataset( + data_vars=dict( + glac_ELA_annual=(coord_order, reg_glac_allgcms_ELA_annual), + crs=np.nan, + ), + coords=coords_dict, + attrs=attrs_dict, + ) + ds.glac_ELA_annual.attrs['long_name'] = ( + 'annual equilibrium line altitude above mean sea level' + ) + ds.glac_ELA_annual.attrs['units'] = 'm' + ds.glac_ELA_annual.attrs['temporal_resolution'] = 'annual' + ds.glac_ELA_annual.attrs['comment'] = ( + 'equilibrium line altitude is the elevation where the climatic mass balance is zero' + ) + ds.glac_ELA_annual.attrs['grid_mapping'] = 'crs' + + # glac_mass_bsl_annual + elif var == 'glac_mass_bsl_annual': + ds = xr.Dataset( + data_vars=dict( + glac_mass_bsl_annual=( + coord_order, + reg_glac_allgcms_mass_bsl_annual, + ), + crs=np.nan, + ), + coords=coords_dict, + attrs=attrs_dict, + ) + ds.glac_mass_bsl_annual.attrs['long_name'] = ( + 'glacier mass below sea leve' + ) + ds.glac_mass_bsl_annual.attrs['units'] = 'kg' + ds.glac_mass_bsl_annual.attrs['temporal_resolution'] = 'annual' + ds.glac_mass_bsl_annual.attrs['comment'] = ( + 'mass of ice below sea level based on area and ice thickness at start of the year' + ) + ds.glac_mass_bsl_annual.attrs['grid_mapping'] = 'crs' + # crs attributes - same for all vars ds.crs.attrs['grid_mapping_name'] = 'latitude_longitude' ds.crs.attrs['longitude_of_prime_meridian'] = 0.0 @@ -800,66 +926,66 @@ def run(args): os.makedirs(vn_fp, exist_ok=True) if realizations[0]: - ds_fn = f'R{str(reg).zfill(2)}_{var}_{gcms[0]}_{scenario}_Batch-{str(batch_start)}-{str(batch_end)}_{calibration}_ba{bias_adj}_{nsets}_{gcm_startyear}_{gcm_endyear}_all.nc'.replace( + ds_fn = f'R{str(reg).zfill(2)}_{var}_{gcms[0]}_{sim_climate_scenario}_Batch-{str(batch_start)}-{str(batch_end)}_{calibration}_ba{bias_adj}_{nsets}_{sim_startyear}_{sim_endyear}_all.nc'.replace( '__', '_' ) else: - ds_fn = f'R{str(reg).zfill(2)}_{var}_{scenario}_Batch-{str(batch_start)}-{str(batch_end)}_{calibration}_ba{bias_adj}_{nsets}_{gcm_startyear}_{gcm_endyear}_all.nc'.replace( + ds_fn = f'R{str(reg).zfill(2)}_{var}_{sim_climate_scenario}_Batch-{str(batch_start)}-{str(batch_end)}_{calibration}_ba{bias_adj}_{nsets}_{sim_startyear}_{sim_endyear}_all.nc'.replace( '__', '_' - ) + ).replace('__', '_') ds.to_netcdf(vn_fp + ds_fn) loop_end = time.time() - print(f'Batch {nbatch} runtime:\t{np.round(loop_end - loop_start, 2)} seconds') + print( + f'Batch {nbatch} ([{str(reg).zfill(2)}.{batch_start}:{str(reg).zfill(2)}.{batch_end}]) runtime:\t{np.round(loop_end - loop_start, 2)} seconds' + ) ### MERGE BATCHES FOR ANNUAL VARS ### - vns = ['glac_mass_annual', 'glac_area_annual'] - - for vn in vns: - if vn in vars: - vn_fp = f'{comppath}glacier_stats/{vn}/{str(reg).zfill(2)}/' + for var in vars: + if 'annual' in var: + var_fp = f'{comppath}glacier_stats/{var}/{str(reg).zfill(2)}/' - fn_merge_list_start = [] + fp_merge_list_start = [] if realizations[0]: - fn_merge_list = glob.glob( - f'{vn_fp}/R{str(reg).zfill(2)}_{vn}_{gcms[0]}_{scenario}_Batch-*_{calibration}_ba{bias_adj}_{nsets}_{gcm_startyear}_{gcm_endyear}_all.nc'.replace( + fp_merge_list = glob.glob( + f'{var_fp}/R{str(reg).zfill(2)}_{var}_{gcms[0]}_{sim_climate_scenario}_Batch-*_{calibration}_ba{bias_adj}_{nsets}_{sim_startyear}_{sim_endyear}_all.nc'.replace( '__', '_' ) ) else: - fn_merge_list = glob.glob( - f'{vn_fp}/R{str(reg).zfill(2)}_{vn}_{scenario}_Batch-*_{calibration}_ba{bias_adj}_{nsets}_{gcm_startyear}_{gcm_endyear}_all.nc'.replace( + fp_merge_list = glob.glob( + f'{var_fp}/R{str(reg).zfill(2)}_{var}_{sim_climate_scenario}_Batch-*_{calibration}_ba{bias_adj}_{nsets}_{sim_startyear}_{sim_endyear}_all.nc'.replace( '__', '_' ) ) - fn_merge_list_start = [int(f.split('-')[-2]) for f in fn_merge_list] + fp_merge_list_start = [int(f.split('-')[-2]) for f in fp_merge_list] - if len(fn_merge_list) > 0: - fn_merge_list = [ - x for _, x in sorted(zip(fn_merge_list_start, fn_merge_list)) + if len(fp_merge_list) > 0: + fp_merge_list = [ + x for _, x in sorted(zip(fp_merge_list_start, fp_merge_list)) ] ds = None - for fn in fn_merge_list: - ds_batch = xr.open_dataset(fn) + for fp in fp_merge_list: + ds_batch = xr.open_dataset(fp) if ds is None: ds = ds_batch else: ds = xr.concat([ds, ds_batch], dim='glacier') # save - ds_fn = ( - fn.split('Batch')[0][:-1] - + f'_{calibration}_ba{bias_adj}_{nsets}_{gcm_startyear}_{gcm_endyear}_all.nc' + ds_fp = ( + fp.split('Batch')[0][:-1] + + f'_{calibration}_ba{bias_adj}_{nsets}_{sim_startyear}_{sim_endyear}_all.nc' ) - ds.to_netcdf(ds_fn) + ds.to_netcdf(ds_fp) ds_batch.close() - for fn in fn_merge_list: - os.remove(fn) + for fp in fp_merge_list: + os.remove(fp) def main(): @@ -867,7 +993,7 @@ def main(): # Set up CLI parser = argparse.ArgumentParser( - description="""description: program for compiling regional stats from the python glacier evolution model (PyGEM)\nnote, this script is not embarrassingly parallel\nit is currently set up to be parallelized by splitting into n jobs based on the number of regions and scenarios scecified\nfor example, the call below could be parallelized into 4 jobs (2 regions x 2 scenarios)\n\nexample call: $python compile_simulations -rgi_region 01 02 -scenario ssp345 ssp585 -gcm_startyear2000 -gcm_endyear 2100 -ncores 4 -vars glac_mass_annual glac_area_annual""", + description="""description: program for compiling regional stats from the python glacier evolution model (PyGEM)\nnote, this script is not embarrassingly parallel\nit is currently set up to be parallelized by splitting into n jobs based on the number of regions and sim_climate_scenarios scecified\nfor example, the call below could be parallelized into 4 jobs (2 regions x 2 sim_climate_scenarios)\n\nexample call: $python compile_simulations -rgi_region 01 02 -sim_climate_scenario ssp345 ssp585 -sim_startyear2000 -sim_endyear 2100 -ncores 4 -vars glac_mass_annual glac_area_annual""", formatter_class=argparse.RawTextHelpFormatter, ) requiredNamed = parser.add_argument_group('required named arguments') @@ -880,7 +1006,7 @@ def main(): help='Randoph Glacier Inventory region (can take multiple, e.g. "1 2 3")', ) requiredNamed.add_argument( - '-gcm_name', + '-sim_climate_name', type=str, default=None, required=True, @@ -888,12 +1014,12 @@ def main(): help='GCM name for which to compile simulations (can take multiple, ex. "ERA5" or "CESM2")', ) parser.add_argument( - '-scenario', + '-sim_climate_scenario', action='store', type=str, default=None, nargs='+', - help='rcp or ssp scenario used for model run (can take multiple, ex. "ssp245 ssp585")', + help='rcp or ssp sim_climate_scenario used for model run (can take multiple, ex. "ssp245 ssp585")', ) parser.add_argument( '-realization', @@ -904,17 +1030,17 @@ def main(): help='realization from large ensemble used for model run (cant take multiple, ex. "r1i1p1f1 r2i1p1f1 r3i1p1f1")', ) parser.add_argument( - '-gcm_startyear', + '-sim_startyear', action='store', type=int, - default=pygem_prms['climate']['gcm_startyear'], + default=pygem_prms['climate']['sim_startyear'], help='start year for the model run', ) parser.add_argument( - '-gcm_endyear', + '-sim_endyear', action='store', type=int, - default=pygem_prms['climate']['gcm_endyear'], + default=pygem_prms['climate']['sim_endyear'], help='start year for the model run', ) parser.add_argument( @@ -948,11 +1074,14 @@ def main(): 'glac_melt_monthly', 'glac_refreeze_monthly', 'glac_frontalablation_monthly', + 'glac_snowline_monthly', 'glac_massbaltotal_monthly', 'glac_prec_monthly', 'glac_mass_monthly', 'glac_mass_annual', 'glac_area_annual', + 'glac_ELA_annual', + 'glac_mass_bsl_annual', ], nargs='+', ) @@ -967,13 +1096,13 @@ def main(): args = parser.parse_args() simpath = args.sim_path region = args.rgi_region01 - gcms = args.gcm_name - scenarios = args.scenario + gcms = args.sim_climate_name + sim_climate_scenarios = args.sim_climate_scenario realizations = args.realization calib = args.option_calibration bias_adj = args.option_bias_adjustment - gcm_startyear = args.gcm_startyear - gcm_endyear = args.gcm_endyear + sim_startyear = args.sim_startyear + sim_endyear = args.sim_endyear vars = args.vars if not simpath: @@ -988,19 +1117,21 @@ def main(): if not isinstance(gcms, list): gcms = [gcms] - if scenarios: - if not isinstance(scenarios, list): - scenarios = [scenarios] + if sim_climate_scenarios: + if not isinstance(sim_climate_scenarios, list): + sim_climate_scenarios = [sim_climate_scenarios] if set(['ERA5', 'ERA-Interim', 'COAWST']) & set(gcms): raise ValueError( - f'Cannot compile present-day and future data simulataneously. A scenario was specified, which does not exist for one of the specified GCMs.\nGCMs: {gcms}\nScenarios: {scenarios}' + f'Cannot compile present-day and future data simulataneously. A sim_climate_scenario was specified, which does not exist for one of the specified GCMs.\nGCMs: {gcms}\nScenarios: {sim_climate_scenarios}' ) else: - scenarios = [''] + sim_climate_scenarios = [''] if set(gcms) - set(['ERA5', 'ERA-Interim', 'COAWST']): raise ValueError( - f'Must specify a scenario for future GCM runs\nGCMs: {gcms}\nscenarios: {scenarios}' + f'Must specify a sim_climate_scenario for future GCM runs\nGCMs: {gcms}\nsim_climate_scenarios: {sim_climate_scenarios}' ) + # no bias adjustment if not a future GCM - ensure bias_adj == 0 + bias_adj = 0 if realizations is None: realizations = [''] @@ -1020,15 +1151,18 @@ def main(): 'glac_melt_monthly', 'glac_refreeze_monthly', 'glac_frontalablation_monthly', + 'glac_snowline_monthly', 'glac_massbaltotal_monthly', 'glac_prec_monthly', 'glac_mass_monthly', 'glac_mass_annual', 'glac_area_annual', + 'glac_ELA_annual', + 'glac_mass_bsl_annual', ] # get number of jobs and split into desired number of cores - njobs = int(len(region) * len(scenarios)) + njobs = int(len(region) * len(sim_climate_scenarios)) # number of cores for parallel processing if args.ncores > 1: num_cores = int(np.min([njobs, args.ncores])) @@ -1042,16 +1176,16 @@ def main(): 'simpath', 'gcms', 'realizations', - 'scenario', + 'sim_climate_scenario', 'calib', 'bias_adj', - 'gcm_startyear', - 'gcm_endyear', + 'sim_startyear', + 'sim_endyear', 'vars', ] i = 0 - # if realizations specified, aggregate all realizations for each gcm and scenario by region - for sce in scenarios: + # if realizations specified, aggregate all realizations for each gcm and sim_climate_scenario by region + for sce in sim_climate_scenarios: for reg in region: list_packed_vars.append( [ @@ -1062,8 +1196,8 @@ def main(): sce, calib, bias_adj, - gcm_startyear, - gcm_endyear, + sim_startyear, + sim_endyear, vars, ] ) diff --git a/pygem/bin/preproc/preproc_wgms_estimate_kp.py b/pygem/bin/preproc/preproc_wgms_estimate_kp.py index de07d8d0..841cdbf3 100644 --- a/pygem/bin/preproc/preproc_wgms_estimate_kp.py +++ b/pygem/bin/preproc/preproc_wgms_estimate_kp.py @@ -213,7 +213,7 @@ def est_kp( # ===== LOAD CLIMATE DATA ===== # Climate class - gcm = class_climate.GCM(name=pygem_prms['climate']['ref_gcm_name']) + gcm = class_climate.GCM(name=pygem_prms['climate']['ref_climate_name']) # Air temperature [degC] gcm_temp, gcm_dates = gcm.importGCMvarnearestneighbor_xarray( @@ -233,7 +233,7 @@ def est_kp( ) # ===== PROCESS THE OBSERVATIONS ====== - prec_cn = pygem_prms['climate']['ref_gcm_name'] + '_prec' + prec_cn = pygem_prms['climate']['ref_climate_name'] + '_prec' wgms_df[prec_cn] = np.nan wgms_df['kp'] = np.nan wgms_df['ndays'] = np.nan diff --git a/pygem/bin/run/run_calibration.py b/pygem/bin/run/run_calibration.py index 0d715d51..fbcbd9c1 100755 --- a/pygem/bin/run/run_calibration.py +++ b/pygem/bin/run/run_calibration.py @@ -59,7 +59,7 @@ def getparser(): Parameters ---------- - ref_gcm_name (optional) : str + ref_climate_name (optional) : str reference gcm name ncores (optional) : int number of cores to use in parallels @@ -95,10 +95,10 @@ def getparser(): help='Randoph Glacier Inventory subregion (either `all` or multiple spaced integers, e.g. `-run_region02 1 2 3`)', ) parser.add_argument( - '-ref_gcm_name', + '-ref_climate_name', action='store', type=str, - default=pygem_prms['climate']['ref_gcm_name'], + default=pygem_prms['climate']['ref_climate_name'], help='reference gcm name', ) parser.add_argument( @@ -122,7 +122,7 @@ def getparser(): type=str, default=None, help='filepath containing list of rgi_glac_number, helpful for running batches on spc', - ), + ) ) parser.add_argument( '-rgi_glac_number', @@ -541,7 +541,7 @@ def run(list_packed_vars): # Unpack variables glac_no = list_packed_vars[1] - gcm_name = list_packed_vars[2] + ref_climate_name = list_packed_vars[2] parser = getparser() args = parser.parse_args() @@ -554,21 +554,20 @@ def run(list_packed_vars): dates_table = modelsetup.datesmodelrun( startyear=args.ref_startyear, endyear=args.ref_endyear, - spinupyears=pygem_prms['climate']['ref_spinupyears'], option_wateryear=pygem_prms['climate']['ref_wateryear'], ) # ===== LOAD CLIMATE DATA ===== # Climate class - assert gcm_name in ['ERA5', 'ERA-Interim'], ( - 'Error: Calibration not set up for ' + gcm_name + assert ref_climate_name in ['ERA5', 'ERA-Interim'], ( + 'Error: Calibration not set up for ' + ref_climate_name ) - gcm = class_climate.GCM(name=gcm_name) + gcm = class_climate.GCM(name=ref_climate_name) # Air temperature [degC] gcm_temp, gcm_dates = gcm.importGCMvarnearestneighbor_xarray( gcm.temp_fn, gcm.temp_vn, main_glac_rgi, dates_table ) - if pygem_prms['mb']['option_ablation'] == 2 and gcm_name in ['ERA5']: + if pygem_prms['mb']['option_ablation'] == 2 and ref_climate_name in ['ERA5']: gcm_tempstd, gcm_dates = gcm.importGCMvarnearestneighbor_xarray( gcm.tempstd_fn, gcm.tempstd_vn, main_glac_rgi, dates_table ) @@ -591,7 +590,7 @@ def run(list_packed_vars): for glac in range(main_glac_rgi.shape[0]): if debug or glac == 0 or glac == main_glac_rgi.shape[0]: print( - gcm_name, + ref_climate_name, ':', main_glac_rgi.loc[main_glac_rgi.index.values[glac], 'RGIId'], ) @@ -3002,13 +3001,13 @@ def main(): ) # Read GCM names from argument parser - gcm_name = args.ref_gcm_name - print('Processing:', gcm_name) + ref_climate_name = args.ref_climate_name + print('Processing:', ref_climate_name) # Pack variables for multiprocessing list_packed_vars = [] for count, glac_no_lst in enumerate(glac_no_lsts): - list_packed_vars.append([count, glac_no_lst, gcm_name]) + list_packed_vars.append([count, glac_no_lst, ref_climate_name]) # Parallel processing if num_cores > 1: print('Processing in parallel with ' + str(num_cores) + ' cores...') diff --git a/pygem/bin/run/run_calibration_frontalablation.py b/pygem/bin/run/run_calibration_frontalablation.py index 7862f794..ccdca552 100644 --- a/pygem/bin/run/run_calibration_frontalablation.py +++ b/pygem/bin/run/run_calibration_frontalablation.py @@ -126,21 +126,19 @@ def reg_calving_flux( dates_table = modelsetup.datesmodelrun( startyear=args.ref_startyear, endyear=args.ref_endyear, - spinupyears=pygem_prms['climate']['ref_spinupyears'], - option_wateryear=pygem_prms['climate']['ref_wateryear'], ) # ===== LOAD CLIMATE DATA ===== # Climate class - assert args.ref_gcm_name in ['ERA5', 'ERA-Interim'], ( - 'Error: Calibration not set up for ' + args.ref_gcm_name + assert args.ref_climate_name in ['ERA5', 'ERA-Interim'], ( + 'Error: Calibration not set up for ' + args.ref_climate_name ) - gcm = class_climate.GCM(name=args.ref_gcm_name) + gcm = class_climate.GCM(name=args.ref_climate_name) # Air temperature [degC] gcm_temp, gcm_dates = gcm.importGCMvarnearestneighbor_xarray( gcm.temp_fn, gcm.temp_vn, main_glac_rgi, dates_table ) - if pygem_prms['mb']['option_ablation'] == 2 and args.ref_gcm_name in ['ERA5']: + if pygem_prms['mb']['option_ablation'] == 2 and args.ref_climate_name in ['ERA5']: gcm_tempstd, gcm_dates = gcm.importGCMvarnearestneighbor_xarray( gcm.tempstd_fn, gcm.tempstd_vn, main_glac_rgi, dates_table ) @@ -3059,10 +3057,10 @@ def main(): nargs='+', ) parser.add_argument( - '-ref_gcm_name', + '-ref_climate_name', action='store', type=str, - default=pygem_prms['climate']['ref_gcm_name'], + default=pygem_prms['climate']['ref_climate_name'], help='reference gcm name', ) parser.add_argument( diff --git a/pygem/bin/run/run_calibration_reg_glena.py b/pygem/bin/run/run_calibration_reg_glena.py index c674562b..503ed0ae 100644 --- a/pygem/bin/run/run_calibration_reg_glena.py +++ b/pygem/bin/run/run_calibration_reg_glena.py @@ -47,7 +47,7 @@ def getparser(): Parameters ---------- - ref_gcm_name (optional) : str + ref_climate_name (optional) : str reference gcm name rgi_glac_number_fn : str filename of .pkl file containing a list of glacier numbers which is used to run batches on the supercomputer @@ -74,10 +74,10 @@ def getparser(): nargs='+', ) parser.add_argument( - '-ref_gcm_name', + '-ref_climate_name', action='store', type=str, - default=pygem_prms['climate']['ref_gcm_name'], + default=pygem_prms['climate']['ref_climate_name'], help='reference gcm name', ) parser.add_argument( @@ -334,16 +334,16 @@ def main(): # ===== LOAD CLIMATE DATA ===== # Climate class - gcm_name = args.ref_gcm_name - assert gcm_name in ['ERA5', 'ERA-Interim'], ( - 'Error: Calibration not set up for ' + gcm_name + sim_climate_name = args.ref_climate_name + assert sim_climate_name in ['ERA5', 'ERA-Interim'], ( + 'Error: Calibration not set up for ' + sim_climate_name ) - gcm = class_climate.GCM(name=gcm_name) + gcm = class_climate.GCM(name=sim_climate_name) # Air temperature [degC] gcm_temp, gcm_dates = gcm.importGCMvarnearestneighbor_xarray( gcm.temp_fn, gcm.temp_vn, main_glac_rgi_subset, dates_table ) - if pygem_prms['mbmod']['option_ablation'] == 2 and gcm_name in ['ERA5']: + if pygem_prms['mbmod']['option_ablation'] == 2 and sim_climate_name in ['ERA5']: gcm_tempstd, gcm_dates = gcm.importGCMvarnearestneighbor_xarray( gcm.tempstd_fn, gcm.tempstd_vn, main_glac_rgi_subset, dates_table ) diff --git a/pygem/bin/run/run_simulation.py b/pygem/bin/run/run_simulation.py index 8f9d8fa8..4a4c45e4 100755 --- a/pygem/bin/run/run_simulation.py +++ b/pygem/bin/run/run_simulation.py @@ -77,9 +77,9 @@ def getparser(): ---------- gcm_list_fn (optional) : str text file that contains the climate data to be used in the model simulation - gcm_name (optional) : str + sim_climate_name (optional) : str gcm name - scenario (optional) : str + sim_climate_scenario (optional) : str representative concentration pathway or shared socioeconomic pathway (ex. 'rcp26', 'ssp585') realization (optional) : str single realization from large ensemble (ex. '1011.001', '1301.020') @@ -128,10 +128,10 @@ def getparser(): help='Randoph Glacier Inventory glacier number (can take multiple)', ) parser.add_argument( - '-ref_gcm_name', + '-ref_climate_name', action='store', type=str, - default=pygem_prms['climate']['ref_gcm_name'], + default=pygem_prms['climate']['ref_climate_name'], help='reference gcm name', ) parser.add_argument( @@ -159,22 +159,22 @@ def getparser(): '-gcm_list_fn', action='store', type=str, - default=pygem_prms['climate']['ref_gcm_name'], + default=pygem_prms['climate']['ref_climate_name'], help='text file full of commands to run (ex. CanESM2 or CESM2)', ) parser.add_argument( - '-gcm_name', + '-sim_climate_name', action='store', type=str, - default=pygem_prms['climate']['gcm_name'], + default=pygem_prms['climate']['sim_climate_name'], help='GCM name used for model run', ) parser.add_argument( - '-scenario', + '-sim_climate_scenario', action='store', type=none_or_value, - default=pygem_prms['climate']['scenario'], - help='rcp or ssp scenario used for model run (ex. rcp26 or ssp585)', + default=pygem_prms['climate']['sim_climate_scenario'], + help='rcp or ssp sim_climate_scenario used for model run (ex. rcp26 or ssp585)', ) parser.add_argument( '-realization', @@ -191,17 +191,17 @@ def getparser(): help='text file full of realizations to run', ) parser.add_argument( - '-gcm_startyear', + '-sim_startyear', action='store', type=int, - default=pygem_prms['climate']['gcm_startyear'], + default=pygem_prms['climate']['sim_startyear'], help='start year for the model run', ) parser.add_argument( - '-gcm_endyear', + '-sim_endyear', action='store', type=int, - default=pygem_prms['climate']['gcm_endyear'], + default=pygem_prms['climate']['sim_endyear'], help='start year for the model run', ) parser.add_argument( @@ -356,15 +356,17 @@ def run(list_packed_vars): args = parser.parse_args() count = list_packed_vars[0] glac_no = list_packed_vars[1] - gcm_name = list_packed_vars[2] + sim_climate_name = list_packed_vars[2] realization = list_packed_vars[3] - if (gcm_name != args.ref_gcm_name) and (args.scenario is None): - scenario = os.path.basename(args.gcm_list_fn).split('_')[1] + if (sim_climate_name != args.ref_climate_name) and ( + args.sim_climate_scenario is None + ): + sim_climate_scenario = os.path.basename(args.gcm_list_fn).split('_')[1] else: - scenario = args.scenario + sim_climate_scenario = args.sim_climate_scenario debug = args.debug if debug: - print(f'scenario:{scenario}') + print(f'sim_climate_scenario:{sim_climate_scenario}') # ===== LOAD GLACIERS ===== main_glac_rgi = modelsetup.selectglaciersrgitable(glac_no=glac_no) @@ -372,50 +374,51 @@ def run(list_packed_vars): # ===== TIME PERIOD ===== # Reference Calibration Period # adjust end year in event that gcm_end year precedes ref_startyear - ref_endyear = min([args.ref_endyear, args.gcm_endyear]) + ref_endyear = min([args.ref_endyear, args.sim_endyear]) dates_table_ref = modelsetup.datesmodelrun( startyear=args.ref_startyear, endyear=ref_endyear, - spinupyears=pygem_prms['climate']['ref_spinupyears'], option_wateryear=pygem_prms['climate']['ref_wateryear'], ) # GCM Full Period (includes reference and simulation periods) dates_table_full = modelsetup.datesmodelrun( - startyear=min([args.ref_startyear, args.gcm_startyear]), - endyear=args.gcm_endyear, - spinupyears=pygem_prms['climate']['gcm_spinupyears'], - option_wateryear=pygem_prms['climate']['gcm_wateryear'], + startyear=min([args.ref_startyear, args.sim_startyear]), + endyear=args.sim_endyear, + option_wateryear=pygem_prms['climate']['sim_wateryear'], ) # GCM Simulation Period dates_table = modelsetup.datesmodelrun( - startyear=args.gcm_startyear, - endyear=args.gcm_endyear, - spinupyears=pygem_prms['climate']['gcm_spinupyears'], - option_wateryear=pygem_prms['climate']['gcm_wateryear'], + startyear=args.sim_startyear, + endyear=args.sim_endyear, + option_wateryear=pygem_prms['climate']['sim_wateryear'], ) if debug: print('ref years:', args.ref_startyear, ref_endyear) - print('sim years:', args.gcm_startyear, args.gcm_endyear) + print('sim years:', args.sim_startyear, args.sim_endyear) # ===== LOAD CLIMATE DATA ===== # Climate class - if gcm_name in ['ERA5', 'ERA-Interim', 'COAWST']: - gcm = class_climate.GCM(name=gcm_name) + if sim_climate_name in ['ERA5', 'ERA-Interim', 'COAWST']: + gcm = class_climate.GCM(name=sim_climate_name) ref_gcm = gcm dates_table_ref = dates_table_full else: # GCM object if realization is None: - gcm = class_climate.GCM(name=gcm_name, scenario=scenario) + gcm = class_climate.GCM( + name=sim_climate_name, sim_climate_scenario=sim_climate_scenario + ) else: gcm = class_climate.GCM( - name=gcm_name, scenario=scenario, realization=realization + name=sim_climate_name, + sim_climate_scenario=sim_climate_scenario, + realization=realization, ) # Reference GCM - ref_gcm = class_climate.GCM(name=args.ref_gcm_name) + ref_gcm = class_climate.GCM(name=args.ref_climate_name) # ----- Select Temperature and Precipitation Data ----- # Air temperature [degC] @@ -445,12 +448,12 @@ def run(list_packed_vars): # ----- Temperature and Precipitation Bias Adjustments ----- # No adjustments - if args.option_bias_adjustment == 0 or gcm_name == args.ref_gcm_name: - if pygem_prms['climate']['gcm_wateryear'] == 'hydro': + if args.option_bias_adjustment == 0 or sim_climate_name == args.ref_climate_name: + if pygem_prms['climate']['sim_wateryear'] == 'hydro': dates_cn = 'wateryear' else: dates_cn = 'year' - sim_idx_start = dates_table_full[dates_cn].to_list().index(args.gcm_startyear) + sim_idx_start = dates_table_full[dates_cn].to_list().index(args.sim_startyear) gcm_elev_adj = gcm_elev gcm_temp_adj = gcm_temp[:, sim_idx_start:] gcm_prec_adj = gcm_prec[:, sim_idx_start:] @@ -465,10 +468,8 @@ def run(list_packed_vars): gcm_temp, dates_table_ref, dates_table_full, - args.gcm_startyear, + args.sim_startyear, args.ref_startyear, - ref_spinupyears=pygem_prms['climate']['ref_spinupyears'], - gcm_spinupyears=pygem_prms['climate']['gcm_spinupyears'], ) # Precipitation bias correction gcm_prec_adj, gcm_elev_adj = gcmbiasadj.prec_biasadj_opt1( @@ -477,10 +478,8 @@ def run(list_packed_vars): gcm_prec, dates_table_ref, dates_table_full, - args.gcm_startyear, + args.sim_startyear, args.ref_startyear, - ref_spinupyears=pygem_prms['climate']['ref_spinupyears'], - gcm_spinupyears=pygem_prms['climate']['gcm_spinupyears'], ) # OPTION 2: Adjust temp and prec using Huss and Hock (2015) elif args.option_bias_adjustment == 2: @@ -491,10 +490,8 @@ def run(list_packed_vars): gcm_temp, dates_table_ref, dates_table_full, - args.gcm_startyear, + args.sim_startyear, args.ref_startyear, - ref_spinupyears=pygem_prms['climate']['ref_spinupyears'], - gcm_spinupyears=pygem_prms['climate']['gcm_spinupyears'], ) # Precipitation bias correction gcm_prec_adj, gcm_elev_adj = gcmbiasadj.prec_biasadj_HH2015( @@ -503,8 +500,6 @@ def run(list_packed_vars): gcm_prec, dates_table_ref, dates_table_full, - ref_spinupyears=pygem_prms['climate']['ref_spinupyears'], - gcm_spinupyears=pygem_prms['climate']['gcm_spinupyears'], ) # OPTION 3: Adjust temp and prec using quantile delta mapping, Cannon et al. (2015) elif args.option_bias_adjustment == 3: @@ -515,10 +510,8 @@ def run(list_packed_vars): gcm_temp, dates_table_ref, dates_table_full, - args.gcm_startyear, + args.sim_startyear, args.ref_startyear, - ref_spinupyears=pygem_prms['climate']['ref_spinupyears'], - gcm_spinupyears=pygem_prms['climate']['gcm_spinupyears'], ) # Precipitation bias correction @@ -528,10 +521,8 @@ def run(list_packed_vars): gcm_prec, dates_table_ref, dates_table_full, - args.gcm_startyear, + args.sim_startyear, args.ref_startyear, - ref_spinupyears=pygem_prms['climate']['ref_spinupyears'], - gcm_spinupyears=pygem_prms['climate']['gcm_spinupyears'], ) # assert that the gcm_elev_adj is not None @@ -542,12 +533,12 @@ def run(list_packed_vars): if pygem_prms['mb']['option_ablation'] != 2: gcm_tempstd = np.zeros((main_glac_rgi.shape[0], dates_table.shape[0])) ref_tempstd = np.zeros((main_glac_rgi.shape[0], dates_table_ref.shape[0])) - elif pygem_prms['mb']['option_ablation'] == 2 and gcm_name in ['ERA5']: + elif pygem_prms['mb']['option_ablation'] == 2 and sim_climate_name in ['ERA5']: gcm_tempstd, gcm_dates = gcm.importGCMvarnearestneighbor_xarray( gcm.tempstd_fn, gcm.tempstd_vn, main_glac_rgi, dates_table ) ref_tempstd = gcm_tempstd - elif pygem_prms['mb']['option_ablation'] == 2 and args.ref_gcm_name in ['ERA5']: + elif pygem_prms['mb']['option_ablation'] == 2 and args.ref_climate_name in ['ERA5']: # Compute temp std based on reference climate data ref_tempstd, ref_dates = ref_gcm.importGCMvarnearestneighbor_xarray( ref_gcm.tempstd_fn, ref_gcm.tempstd_vn, main_glac_rgi, dates_table_ref @@ -561,7 +552,7 @@ def run(list_packed_vars): ref_tempstd = np.zeros((main_glac_rgi.shape[0], dates_table_ref.shape[0])) # Lapse rate - if gcm_name in ['ERA-Interim', 'ERA5']: + if sim_climate_name in ['ERA-Interim', 'ERA5']: gcm_lr, gcm_dates = gcm.importGCMvarnearestneighbor_xarray( gcm.lr_fn, gcm.lr_vn, main_glac_rgi, dates_table ) @@ -576,7 +567,7 @@ def run(list_packed_vars): ref_lr, dates_table_ref, dates_table_full, - args.gcm_startyear, + args.sim_startyear, args.ref_startyear, ) @@ -597,7 +588,7 @@ def run(list_packed_vars): for glac in range(main_glac_rgi.shape[0]): if glac == 0: print( - gcm_name, + sim_climate_name, ':', main_glac_rgi.loc[main_glac_rgi.index.values[glac], 'RGIId'], ) @@ -642,13 +633,6 @@ def run(list_packed_vars): } gdir_ref.dates_table = dates_table_ref - # Add climate data to glacier directory - if pygem_prms['climate']['hindcast'] == True: - gcm_temp_adj = gcm_temp_adj[::-1] - gcm_tempstd = gcm_tempstd[::-1] - gcm_prec_adj = gcm_prec_adj[::-1] - gcm_lr = gcm_lr[::-1] - gdir.historical_climate = { 'elev': gcm_elev_adj[glac], 'temp': gcm_temp_adj[glac, :], @@ -864,7 +848,7 @@ def run(list_packed_vars): ] # Time attributes and values - if pygem_prms['climate']['gcm_wateryear'] == 'hydro': + if pygem_prms['climate']['sim_wateryear'] == 'hydro': annual_columns = np.unique(dates_table['wateryear'].values)[ 0 : int(dates_table.shape[0] / 12) ] @@ -873,9 +857,7 @@ def run(list_packed_vars): 0 : int(dates_table.shape[0] / 12) ] # append additional year to year_values to account for mass and area at end of period - year_values = annual_columns[ - pygem_prms['climate']['gcm_spinupyears'] : annual_columns.shape[0] - ] + year_values = annual_columns year_values = np.concatenate( (year_values, np.array([annual_columns[-1] + 1])) ) @@ -1167,11 +1149,15 @@ def run(list_packed_vars): + '/Output/simulations/fail-exceed_domain/' + reg_str + '/' - + gcm_name + + sim_climate_name + '/' ) - if gcm_name not in ['ERA-Interim', 'ERA5', 'COAWST']: - fail_domain_fp += scenario + '/' + if sim_climate_name not in [ + 'ERA-Interim', + 'ERA5', + 'COAWST', + ]: + fail_domain_fp += sim_climate_scenario + '/' if not os.path.exists(fail_domain_fp): os.makedirs(fail_domain_fp, exist_ok=True) txt_fn_fail = glacier_str + '-sim_failed.txt' @@ -1198,9 +1184,6 @@ def run(list_packed_vars): fs=fs, is_tidewater=gdir.is_tidewater, water_level=water_level, - spinupyears=pygem_prms['climate'][ - 'ref_spinupyears' - ], ) _, diag = ev_model.run_until_and_store(nyears) ev_model.mb_model.glac_wide_volume_annual = ( @@ -1374,11 +1357,15 @@ def run(list_packed_vars): + '/Output/simulations/fail-exceed_domain/' + reg_str + '/' - + gcm_name + + sim_climate_name + '/' ) - if gcm_name not in ['ERA-Interim', 'ERA5', 'COAWST']: - fail_domain_fp += scenario + '/' + if sim_climate_name not in [ + 'ERA-Interim', + 'ERA5', + 'COAWST', + ]: + fail_domain_fp += sim_climate_scenario + '/' if not os.path.exists(fail_domain_fp): os.makedirs(fail_domain_fp, exist_ok=True) txt_fn_fail = glacier_str + '-sim_failed.txt' @@ -1406,7 +1393,7 @@ def run(list_packed_vars): option_areaconstant=True, ) # ----- MODEL RUN WITH CONSTANT GLACIER AREA ----- - years = np.arange(args.gcm_startyear, args.gcm_endyear + 1) + years = np.arange(args.sim_startyear, args.sim_endyear + 1) mb_all = [] for year in years - years[0]: mb_annual = mbmod.get_annual_mb( @@ -1768,14 +1755,14 @@ def run(list_packed_vars): glacier_rgi_table=glacier_rgi_table, dates_table=dates_table, nsims=1, - gcm_name=gcm_name, - scenario=scenario, + sim_climate_name=sim_climate_name, + sim_climate_scenario=sim_climate_scenario, realization=realization, modelprms=modelprms, ref_startyear=args.ref_startyear, ref_endyear=ref_endyear, - gcm_startyear=args.gcm_startyear, - gcm_endyear=args.gcm_endyear, + sim_startyear=args.sim_startyear, + sim_endyear=args.sim_endyear, option_calibration=args.option_calibration, option_bias_adjustment=args.option_bias_adjustment, ) @@ -1867,14 +1854,14 @@ def run(list_packed_vars): glacier_rgi_table=glacier_rgi_table, dates_table=dates_table, nsims=nsims, - gcm_name=gcm_name, - scenario=scenario, + sim_climate_name=sim_climate_name, + sim_climate_scenario=sim_climate_scenario, realization=realization, modelprms=modelprms, ref_startyear=args.ref_startyear, ref_endyear=ref_endyear, - gcm_startyear=args.gcm_startyear, - gcm_endyear=args.gcm_endyear, + sim_startyear=args.sim_startyear, + sim_endyear=args.sim_endyear, option_calibration=args.option_calibration, option_bias_adjustment=args.option_bias_adjustment, ) @@ -2092,14 +2079,14 @@ def run(list_packed_vars): nsims=1, nbins=surface_h_initial.shape[0], binned_components=args.export_binned_components, - gcm_name=gcm_name, - scenario=scenario, + sim_climate_name=sim_climate_name, + sim_climate_scenario=sim_climate_scenario, realization=realization, modelprms=modelprms, ref_startyear=args.ref_startyear, ref_endyear=ref_endyear, - gcm_startyear=args.gcm_startyear, - gcm_endyear=args.gcm_endyear, + sim_startyear=args.sim_startyear, + sim_endyear=args.sim_endyear, option_calibration=args.option_calibration, option_bias_adjustment=args.option_bias_adjustment, ) @@ -2170,14 +2157,14 @@ def run(list_packed_vars): nsims=nsims, nbins=surface_h_initial.shape[0], binned_components=args.export_binned_components, - gcm_name=gcm_name, - scenario=scenario, + sim_climate_name=sim_climate_name, + sim_climate_scenario=sim_climate_scenario, realization=realization, modelprms=modelprms, ref_startyear=args.ref_startyear, ref_endyear=ref_endyear, - gcm_startyear=args.gcm_startyear, - gcm_endyear=args.gcm_endyear, + sim_startyear=args.sim_startyear, + sim_endyear=args.sim_endyear, option_calibration=args.option_calibration, option_bias_adjustment=args.option_bias_adjustment, ) @@ -2259,11 +2246,11 @@ def run(list_packed_vars): + '/Output/simulations/failed/' + reg_str + '/' - + gcm_name + + sim_climate_name + '/' ) - if gcm_name not in ['ERA-Interim', 'ERA5', 'COAWST']: - fail_fp += scenario + '/' + if sim_climate_name not in ['ERA-Interim', 'ERA5', 'COAWST']: + fail_fp += sim_climate_scenario + '/' if not os.path.exists(fail_fp): os.makedirs(fail_fp, exist_ok=True) txt_fn_fail = glacier_str + '-sim_failed.txt' @@ -2286,8 +2273,8 @@ def main(): assert args.ref_startyear < args.ref_endyear, ( f'ref_startyear [{args.ref_startyear}] must be less than ref_endyear [{args.ref_endyear}]' ) - assert args.gcm_startyear < args.gcm_endyear, ( - f'gcm_startyear [{args.gcm_startyear}] must be less than gcm_endyear [{args.gcm_endyear}]' + assert args.sim_startyear < args.sim_endyear, ( + f'sim_startyear [{args.sim_startyear}] must be less than sim_endyear [{args.sim_endyear}]' ) except AssertionError as err: print('error: ', err) @@ -2324,17 +2311,17 @@ def main(): ) # Read GCM names from argument parser - gcm_name = args.gcm_list_fn - if args.gcm_name is not None: - gcm_list = [args.gcm_name] - scenario = args.scenario - elif args.gcm_list_fn == args.ref_gcm_name: - gcm_list = [args.ref_gcm_name] - scenario = args.scenario + sim_climate_name = args.gcm_list_fn + if args.sim_climate_name is not None: + gcm_list = [args.sim_climate_name] + sim_climate_scenario = args.sim_climate_scenario + elif args.gcm_list_fn == args.ref_climate_name: + gcm_list = [args.ref_climate_name] + sim_climate_scenario = args.sim_climate_scenario else: with open(args.gcm_list_fn, 'r') as gcm_fn: gcm_list = gcm_fn.read().splitlines() - scenario = os.path.basename(args.gcm_list_fn).split('_')[1] + sim_climate_scenario = os.path.basename(args.gcm_list_fn).split('_')[1] print('Found %d gcms to process' % (len(gcm_list))) # Read realizations from argument parser @@ -2352,20 +2339,24 @@ def main(): # If passing this through the list_packed_vars, then don't go back and get from arg parser again! # Loop through all GCMs - for gcm_name in gcm_list: - if args.scenario is None: - print('Processing:', gcm_name) - elif args.scenario is not None: - print('Processing:', gcm_name, scenario) + for sim_climate_name in gcm_list: + if args.sim_climate_scenario is None: + print('Processing:', sim_climate_name) + elif args.sim_climate_scenario is not None: + print('Processing:', sim_climate_name, sim_climate_scenario) # Pack variables for multiprocessing list_packed_vars = [] if realizations is not None: for realization in realizations: for count, glac_no_lst in enumerate(glac_no_lsts): - list_packed_vars.append([count, glac_no_lst, gcm_name, realization]) + list_packed_vars.append( + [count, glac_no_lst, sim_climate_name, realization] + ) else: for count, glac_no_lst in enumerate(glac_no_lsts): - list_packed_vars.append([count, glac_no_lst, gcm_name, realizations]) + list_packed_vars.append( + [count, glac_no_lst, sim_climate_name, realizations] + ) print('Processing with ' + str(num_cores) + ' cores...') # Parallel processing diff --git a/pygem/class_climate.py b/pygem/class_climate.py index d0193e09..e39e19dd 100755 --- a/pygem/class_climate.py +++ b/pygem/class_climate.py @@ -32,13 +32,13 @@ class GCM: ---------- name : str name of climate dataset. - scenario : str - rcp or ssp scenario (example: 'rcp26' or 'ssp585') + sim_climate_scenario : str + rcp or ssp sim_climate_scenario (example: 'rcp26' or 'ssp585') realization : str realization from large ensemble (example: '1011.001' or '1301.020') """ - def __init__(self, name=str(), scenario=str(), realization=None): + def __init__(self, name=str(), sim_climate_scenario=str(), realization=None): """ Add variable name and specific properties associated with each gcm. """ @@ -51,10 +51,10 @@ def __init__(self, name=str(), scenario=str(), realization=None): # Source of climate data self.name = name - # If multiple realizations from each model+scenario are being used, + # If multiple realizations from each model+sim_climate_scenario are being used, # then self.realization = realization. # Otherwise, the realization attribute is not considered for single - # realization model+scenario simulations. + # realization model+sim_climate_scenario simulations. if realization is not None: self.realization = realization @@ -72,7 +72,7 @@ def __init__(self, name=str(), scenario=str(), realization=None): self.temp_fn = ( self.temp_vn + '_mon_' - + scenario + + sim_climate_scenario + '_' + name + '-' @@ -82,7 +82,7 @@ def __init__(self, name=str(), scenario=str(), realization=None): self.prec_fn = ( self.prec_vn + '_mon_' - + scenario + + sim_climate_scenario + '_' + name + '-' @@ -90,26 +90,31 @@ def __init__(self, name=str(), scenario=str(), realization=None): + '.cam.h0.1980-2100.nc' ) self.elev_fn = ( - self.elev_vn + '_fx_' + scenario + '_' + name + '.cam.h0.nc' + self.elev_vn + + '_fx_' + + sim_climate_scenario + + '_' + + name + + '.cam.h0.nc' ) # Variable filepaths self.var_fp = ( pygem_prms['root'] + pygem_prms['climate']['paths']['cesm2_relpath'] - + scenario + + sim_climate_scenario + pygem_prms['climate']['paths']['cesm2_fp_var_ending'] ) self.fx_fp = ( pygem_prms['root'] + pygem_prms['climate']['paths']['cesm2_relpath'] - + scenario + + sim_climate_scenario + pygem_prms['climate']['paths']['cesm2_fp_fx_ending'] ) # Extra information self.timestep = pygem_prms['time']['timestep'] self.rgi_lat_colname = pygem_prms['rgi']['rgi_lat_colname'] self.rgi_lon_colname = pygem_prms['rgi']['rgi_lon_colname'] - self.scenario = scenario + self.sim_climate_scenario = sim_climate_scenario # Set parameters for GFDL SPEAR Large Ensemble elif self.name == 'GFDL-SPEAR-MED': @@ -125,7 +130,7 @@ def __init__(self, name=str(), scenario=str(), realization=None): self.temp_fn = ( self.temp_vn + '_mon_' - + scenario + + sim_climate_scenario + '_' + name + '-' @@ -135,32 +140,34 @@ def __init__(self, name=str(), scenario=str(), realization=None): self.prec_fn = ( self.prec_vn + '_mon_' - + scenario + + sim_climate_scenario + '_' + name + '-' + realization + 'i1p1f1_gr3_1980-2100.nc' ) - self.elev_fn = self.elev_vn + '_fx_' + scenario + '_' + name + '.nc' + self.elev_fn = ( + self.elev_vn + '_fx_' + sim_climate_scenario + '_' + name + '.nc' + ) # Variable filepaths self.var_fp = ( pygem_prms['root'] + pygem_prms['climate']['paths']['gfdl_relpath'] - + scenario + + sim_climate_scenario + pygem_prms['climate']['paths']['gfdl_fp_var_ending'] ) self.fx_fp = ( pygem_prms['root'] + pygem_prms['climate']['paths']['gfdl_relpath'] - + scenario + + sim_climate_scenario + pygem_prms['climate']['paths']['gfdl_fp_fx_ending'] ) # Extra information self.timestep = pygem_prms['time']['timestep'] self.rgi_lat_colname = pygem_prms['rgi']['rgi_lat_colname'] self.rgi_lon_colname = pygem_prms['rgi']['rgi_lon_colname'] - self.scenario = scenario + self.sim_climate_scenario = sim_climate_scenario else: self.realization = [] @@ -223,7 +230,7 @@ def __init__(self, name=str(), scenario=str(), realization=None): self.rgi_lon_colname = pygem_prms['rgi']['rgi_lon_colname'] # Standardized CMIP5 format (GCM/RCP) - elif 'rcp' in scenario: + elif 'rcp' in sim_climate_scenario: # Variable names self.temp_vn = 'tas' self.prec_vn = 'pr' @@ -233,25 +240,40 @@ def __init__(self, name=str(), scenario=str(), realization=None): self.time_vn = 'time' # Variable filenames self.temp_fn = ( - self.temp_vn + '_mon_' + name + '_' + scenario + '_r1i1p1_native.nc' + self.temp_vn + + '_mon_' + + name + + '_' + + sim_climate_scenario + + '_r1i1p1_native.nc' ) self.prec_fn = ( - self.prec_vn + '_mon_' + name + '_' + scenario + '_r1i1p1_native.nc' + self.prec_vn + + '_mon_' + + name + + '_' + + sim_climate_scenario + + '_r1i1p1_native.nc' ) self.elev_fn = ( - self.elev_vn + '_fx_' + name + '_' + scenario + '_r0i0p0.nc' + self.elev_vn + + '_fx_' + + name + + '_' + + sim_climate_scenario + + '_r0i0p0.nc' ) # Variable filepaths self.var_fp = ( pygem_prms['root'] + pygem_prms['climate']['paths']['cmip5_relpath'] - + scenario + + sim_climate_scenario + pygem_prms['climate']['paths']['cmip5_fp_var_ending'] ) self.fx_fp = ( pygem_prms['root'] + pygem_prms['climate']['paths']['cmip5_relpath'] - + scenario + + sim_climate_scenario + pygem_prms['climate']['paths']['cmip5_fp_fx_ending'] ) if not os.path.exists(self.var_fp) and os.path.exists( @@ -276,10 +298,10 @@ def __init__(self, name=str(), scenario=str(), realization=None): self.timestep = pygem_prms['time']['timestep'] self.rgi_lat_colname = pygem_prms['rgi']['rgi_lat_colname'] self.rgi_lon_colname = pygem_prms['rgi']['rgi_lon_colname'] - self.scenario = scenario + self.sim_climate_scenario = sim_climate_scenario # Standardized CMIP6 format (GCM/SSP) - elif 'ssp' in scenario: + elif 'ssp' in sim_climate_scenario: # Variable names self.temp_vn = 'tas' self.prec_vn = 'pr' @@ -289,10 +311,20 @@ def __init__(self, name=str(), scenario=str(), realization=None): self.time_vn = 'time' # Variable filenames self.temp_fn = ( - name + '_' + scenario + '_r1i1p1f1_' + self.temp_vn + '.nc' + name + + '_' + + sim_climate_scenario + + '_r1i1p1f1_' + + self.temp_vn + + '.nc' ) self.prec_fn = ( - name + '_' + scenario + '_r1i1p1f1_' + self.prec_vn + '.nc' + name + + '_' + + sim_climate_scenario + + '_r1i1p1f1_' + + self.prec_vn + + '.nc' ) self.elev_fn = name + '_' + self.elev_vn + '.nc' # Variable filepaths @@ -312,7 +344,7 @@ def __init__(self, name=str(), scenario=str(), realization=None): self.timestep = pygem_prms['time']['timestep'] self.rgi_lat_colname = pygem_prms['rgi']['rgi_lat_colname'] self.rgi_lon_colname = pygem_prms['rgi']['rgi_lon_colname'] - self.scenario = scenario + self.sim_climate_scenario = sim_climate_scenario def importGCMfxnearestneighbor_xarray(self, filename, vn, main_glac_rgi): """ diff --git a/pygem/gcmbiasadj.py b/pygem/gcmbiasadj.py index eb7e2987..c894ea2e 100755 --- a/pygem/gcmbiasadj.py +++ b/pygem/gcmbiasadj.py @@ -73,7 +73,7 @@ def temp_biasadj_HH2015( gcm_temp, dates_table_ref, dates_table, - gcm_startyear, + sim_startyear, ref_startyear, ref_spinupyears=0, gcm_spinupyears=0, @@ -146,14 +146,14 @@ def temp_biasadj_HH2015( # if/else statement for whether or not the full GCM period is the same as the simulation period # create GCM subset for applying bias-correction (e.g., 2000-2100), # that does not include the earlier reference years (e.g., 1981-2000) - if gcm_startyear == ref_startyear: + if sim_startyear == ref_startyear: bc_temp = gcm_temp else: - if pygem_prms['climate']['gcm_wateryear'] == 'hydro': + if pygem_prms['climate']['sim_wateryear'] == 'hydro': dates_cn = 'wateryear' else: dates_cn = 'year' - sim_idx_start = dates_table[dates_cn].to_list().index(gcm_startyear) + sim_idx_start = dates_table[dates_cn].to_list().index(sim_startyear) bc_temp = gcm_temp[:, sim_idx_start:] # Monthly temperature bias adjusted according to monthly average @@ -184,7 +184,7 @@ def temp_biasadj_HH2015( :, gcm_subset_idx_start : gcm_subset_idx_end + 1 ][:, ref_spinupyears * 12 :] - if gcm_startyear == ref_startyear: + if sim_startyear == ref_startyear: if debug: print( ( @@ -221,7 +221,7 @@ def prec_biasadj_HH2015( gcm_prec, dates_table_ref, dates_table, - gcm_startyear, + sim_startyear, ref_startyear, ref_spinupyears=0, gcm_spinupyears=0, @@ -276,14 +276,14 @@ def prec_biasadj_HH2015( # if/else statement for whether or not the full GCM period is the same as the simulation period # create GCM subset for applying bias-correction (e.g., 2000-2100), # that does not include the earlier reference years (e.g., 1985-2000) - if gcm_startyear == ref_startyear: + if sim_startyear == ref_startyear: bc_prec = gcm_prec else: - if pygem_prms['climate']['gcm_wateryear'] == 'hydro': + if pygem_prms['climate']['sim_wateryear'] == 'hydro': dates_cn = 'wateryear' else: dates_cn = 'year' - sim_idx_start = dates_table[dates_cn].to_list().index(gcm_startyear) + sim_idx_start = dates_table[dates_cn].to_list().index(sim_startyear) bc_prec = gcm_prec[:, sim_idx_start:] # Bias adjusted precipitation accounting for differences in monthly mean @@ -320,7 +320,7 @@ def prec_biasadj_opt1( gcm_prec, dates_table_ref, dates_table, - gcm_startyear, + sim_startyear, ref_startyear, ref_spinupyears=0, gcm_spinupyears=0, @@ -375,14 +375,14 @@ def prec_biasadj_opt1( # if/else statement for whether or not the full GCM period is the same as the simulation period # create GCM subset for applying bias-correction (e.g., 2000-2100), # that does not include the earlier reference years (e.g., 1985-2000) - if gcm_startyear == ref_startyear: + if sim_startyear == ref_startyear: bc_prec = gcm_prec else: - if pygem_prms['climate']['gcm_wateryear'] == 'hydro': + if pygem_prms['climate']['sim_wateryear'] == 'hydro': dates_cn = 'wateryear' else: dates_cn = 'year' - sim_idx_start = dates_table[dates_cn].to_list().index(gcm_startyear) + sim_idx_start = dates_table[dates_cn].to_list().index(sim_startyear) bc_prec = gcm_prec[:, sim_idx_start:] # Bias adjusted precipitation accounting for differences in monthly mean @@ -476,7 +476,7 @@ def temp_biasadj_QDM( gcm_temp, dates_table_ref, dates_table, - gcm_startyear, + sim_startyear, ref_startyear, ref_spinupyears=0, gcm_spinupyears=0, @@ -531,14 +531,14 @@ def temp_biasadj_QDM( # if/else statement for whether or not the full GCM period is the same as the simulation period # create GCM subset for applying bias-correction (e.g., 2000-2100), # that does not include the earlier reference years (e.g., 1981-2000) - if gcm_startyear == ref_startyear: + if sim_startyear == ref_startyear: bc_temp = gcm_temp else: - if pygem_prms['climate']['gcm_wateryear'] == 'hydro': + if pygem_prms['climate']['sim_wateryear'] == 'hydro': dates_cn = 'wateryear' else: dates_cn = 'year' - sim_idx_start = dates_table[dates_cn].to_list().index(gcm_startyear) + sim_idx_start = dates_table[dates_cn].to_list().index(sim_startyear) bc_temp = gcm_temp[:, sim_idx_start:] # create an empty array for the bias-corrected GCM data @@ -600,7 +600,7 @@ def prec_biasadj_QDM( gcm_prec, dates_table_ref, dates_table, - gcm_startyear, + sim_startyear, ref_startyear, ref_spinupyears=0, gcm_spinupyears=0, @@ -656,14 +656,14 @@ def prec_biasadj_QDM( # if/else statement for whether or not the full GCM period is the same as the simulation period # create GCM subset for applying bias-correction (e.g., 2000-2100), # that does not include the earlier reference years (e.g., 1981-2000) - if gcm_startyear == ref_startyear: + if sim_startyear == ref_startyear: bc_prec = gcm_prec else: - if pygem_prms['climate']['gcm_wateryear'] == 'hydro': + if pygem_prms['climate']['sim_wateryear'] == 'hydro': dates_cn = 'wateryear' else: dates_cn = 'year' - sim_idx_start = dates_table[dates_cn].to_list().index(gcm_startyear) + sim_idx_start = dates_table[dates_cn].to_list().index(sim_startyear) bc_prec = gcm_prec[:, sim_idx_start:] # create an empty array for the bias-corrected GCM data @@ -712,7 +712,7 @@ def prec_biasadj_QDM( def monthly_avg_array_rolled( - ref_array, dates_table_ref, dates_table, gcm_startyear, ref_startyear + ref_array, dates_table_ref, dates_table, sim_startyear, ref_startyear ): """Monthly average array from reference data rolled to ensure proper months @@ -743,12 +743,12 @@ def monthly_avg_array_rolled( # if/else statement for whether or not the full GCM period is the same as the simulation period # create GCM subset for applying bias-correction (e.g., 2000-2100), # that does not include the earlier reference years (e.g., 1981-2000) - if gcm_startyear != ref_startyear: - if pygem_prms['climate']['gcm_wateryear'] == 'hydro': + if sim_startyear != ref_startyear: + if pygem_prms['climate']['sim_wateryear'] == 'hydro': dates_cn = 'wateryear' else: dates_cn = 'year' - sim_idx_start = dates_table[dates_cn].to_list().index(gcm_startyear) + sim_idx_start = dates_table[dates_cn].to_list().index(sim_startyear) gcm_array = gcm_array[:, sim_idx_start:] return gcm_array diff --git a/pygem/glacierdynamics.py b/pygem/glacierdynamics.py index a66a5edf..5ba4b8ad 100755 --- a/pygem/glacierdynamics.py +++ b/pygem/glacierdynamics.py @@ -622,7 +622,6 @@ def _massredistributionHuss( glac_idx_initial, heights, debug=False, - hindcast=0, sec_in_year=365 * 24 * 3600, ): """ @@ -665,10 +664,6 @@ def _massredistributionHuss( glac_bin_massbalclim_annual * sec_in_year * glacier_area_t0 ).sum() - # For hindcast simulations, volume change is the opposite - if hindcast == 1: - glacier_volumechange = -1 * glacier_volumechange - if debug: print('\nDebugging Mass Redistribution Huss function\n') print('glacier volume change:', glacier_volumechange) diff --git a/pygem/massbalance.py b/pygem/massbalance.py index 9d2c439c..85475f8f 100644 --- a/pygem/massbalance.py +++ b/pygem/massbalance.py @@ -36,7 +36,6 @@ def __init__( modelprms, glacier_rgi_table, option_areaconstant=False, - hindcast=pygem_prms['climate']['hindcast'], frontalablation_k=None, debug=pygem_prms['debug']['mb'], debug_refreeze=pygem_prms['debug']['refreeze'], @@ -63,8 +62,6 @@ def __init__( option to turn on print statements for development or debugging of code debug_refreeze : Boolean option to turn on print statements for development/debugging of refreezing code - hindcast : Boolean - switch to run the model in reverse or not (may be irrelevant after converting to OGGM's setup) """ if debug: print('\n\nDEBUGGING MASS BALANCE FUNCTION\n\n') @@ -106,12 +103,6 @@ def __init__( self.glacier_gcm_lrgcm = gdir.historical_climate['lr'] self.glacier_gcm_lrglac = gdir.historical_climate['lr'] - if pygem_prms['climate']['hindcast'] == True: - self.glacier_gcm_prec = self.glacier_gcm_prec[::-1] - self.glacier_gcm_temp = self.glacier_gcm_temp[::-1] - self.glacier_gcm_lrgcm = self.glacier_gcm_lrgcm[::-1] - self.glacier_gcm_lrglac = self.glacier_gcm_lrglac[::-1] - self.repeat_period = repeat_period # Variables to store (consider storing in xarray) @@ -242,8 +233,8 @@ def get_annual_mb( year = int(year) if self.repeat_period: year = year % ( - pygem_prms['climate']['gcm_endyear'] - - pygem_prms['climate']['gcm_startyear'] + pygem_prms['climate']['sim_endyear'] + - pygem_prms['climate']['sim_startyear'] ) fl = fls[fl_id] @@ -309,46 +300,44 @@ def get_annual_mb( # if (pygem_prms['time']['timestep'] == 'monthly') and (glac_idx_t0.shape[0] != 0): # AIR TEMPERATURE: Downscale the gcm temperature [deg C] to each bin - if pygem_prms['mb']['option_temp2bins'] == 1: - # Downscale using gcm and glacier lapse rates - # T_bin = T_gcm + lr_gcm * (z_ref - z_gcm) + lr_glac * (z_bin - z_ref) + tempchange - self.bin_temp[:, 12 * year : 12 * (year + 1)] = ( - self.glacier_gcm_temp[12 * year : 12 * (year + 1)] - + self.glacier_gcm_lrgcm[12 * year : 12 * (year + 1)] - * ( - self.glacier_rgi_table.loc[ - pygem_prms['mb']['option_elev_ref_downscale'] - ] - - self.glacier_gcm_elev - ) - + self.glacier_gcm_lrglac[12 * year : 12 * (year + 1)] + # Downscale using gcm and glacier lapse rates + # T_bin = T_gcm + lr_gcm * (z_ref - z_gcm) + lr_glac * (z_bin - z_ref) + tempchange + self.bin_temp[:, 12 * year : 12 * (year + 1)] = ( + self.glacier_gcm_temp[12 * year : 12 * (year + 1)] + + self.glacier_gcm_lrgcm[12 * year : 12 * (year + 1)] + * ( + self.glacier_rgi_table.loc[ + pygem_prms['mb']['option_elev_ref_downscale'] + ] + - self.glacier_gcm_elev + ) + + self.glacier_gcm_lrglac[12 * year : 12 * (year + 1)] + * ( + heights + - self.glacier_rgi_table.loc[ + pygem_prms['mb']['option_elev_ref_downscale'] + ] + )[:, np.newaxis] + + self.modelprms['tbias'] + ) + + # PRECIPITATION/ACCUMULATION: Downscale the precipitation (liquid and solid) to each bin + # Precipitation using precipitation factor and precipitation gradient + # P_bin = P_gcm * prec_factor * (1 + prec_grad * (z_bin - z_ref)) + bin_precsnow[:, 12 * year : 12 * (year + 1)] = ( + self.glacier_gcm_prec[12 * year : 12 * (year + 1)] + * self.modelprms['kp'] + * ( + 1 + + self.modelprms['precgrad'] * ( heights - self.glacier_rgi_table.loc[ pygem_prms['mb']['option_elev_ref_downscale'] ] - )[:, np.newaxis] - + self.modelprms['tbias'] - ) - - # PRECIPITATION/ACCUMULATION: Downscale the precipitation (liquid and solid) to each bin - if pygem_prms['mb']['option_prec2bins'] == 1: - # Precipitation using precipitation factor and precipitation gradient - # P_bin = P_gcm * prec_factor * (1 + prec_grad * (z_bin - z_ref)) - bin_precsnow[:, 12 * year : 12 * (year + 1)] = ( - self.glacier_gcm_prec[12 * year : 12 * (year + 1)] - * self.modelprms['kp'] - * ( - 1 - + self.modelprms['precgrad'] - * ( - heights - - self.glacier_rgi_table.loc[ - pygem_prms['mb']['option_elev_ref_downscale'] - ] - ) - )[:, np.newaxis] - ) + ) + )[:, np.newaxis] + ) # Option to adjust prec of uppermost 25% of glacier for wind erosion and reduced moisture content if pygem_prms['mb']['option_preclimit'] == 1: # Elevation range based on all flowlines diff --git a/pygem/output.py b/pygem/output.py index 0ebc658a..eb151808 100644 --- a/pygem/output.py +++ b/pygem/output.py @@ -53,10 +53,10 @@ class single_glacier: DataFrame containing metadata and characteristics of the glacier from the Randolph Glacier Inventory. dates_table : pd.DataFrame DataFrame containing the time series of dates associated with the model output. - gcm_name : str + sim_climate_name : str Name of the General Circulation Model (GCM) used for climate forcing. - scenario : str - Emission or climate scenario under which the simulation is run. + sim_climate_scenario : str + Emission or climate sim_climate_scenario under which the simulation is run. realization : str Specific realization or ensemble member of the GCM simulation. nsims : int @@ -67,9 +67,9 @@ class single_glacier: Start year of the reference period for model calibration or comparison. ref_endyear : int End year of the reference period for model calibration or comparison. - gcm_startyear : int + sim_startyear : int Start year of the GCM forcing data used in the simulation. - gcm_endyear : int + sim_endyear : int End year of the GCM forcing data used in the simulation. option_calibration : str Model calibration method. @@ -79,15 +79,15 @@ class single_glacier: glacier_rgi_table: pd.DataFrame dates_table: pd.DataFrame - gcm_name: str - scenario: str + sim_climate_name: str + sim_climate_scenario: str realization: str nsims: int modelprms: dict ref_startyear: int ref_endyear: int - gcm_startyear: int - gcm_endyear: int + sim_startyear: int + sim_endyear: int option_calibration: str option_bias_adjustment: str @@ -123,23 +123,23 @@ def set_fn(self, outfn=None): if outfn: self.outfn = outfn else: - self.outfn = self.glacier_str + '_' + self.gcm_name + '_' - if self.scenario: - self.outfn += f'{self.scenario}_' + self.outfn = self.glacier_str + '_' + self.sim_climate_name + '_' + if self.sim_climate_scenario: + self.outfn += f'{self.sim_climate_scenario}_' if self.realization: self.outfn += f'{self.realization}_' if self.option_calibration: self.outfn += f'{self.option_calibration}_' else: self.outfn += f'kp{self.modelprms["kp"]}_ddfsnow{self.modelprms["ddfsnow"]}_tbias{self.modelprms["tbias"]}_' - if self.gcm_name not in ['ERA-Interim', 'ERA5', 'COAWST']: + if self.sim_climate_name not in ['ERA-Interim', 'ERA5', 'COAWST']: self.outfn += f'ba{self.option_bias_adjustment}_' else: self.outfn += 'ba0_' if self.option_calibration: self.outfn += 'SETS_' - self.outfn += f'{self.gcm_startyear}_' - self.outfn += f'{self.gcm_endyear}_' + self.outfn += f'{self.sim_startyear}_' + self.outfn += f'{self.sim_endyear}_' def get_fn(self): """Return the output dataset filename.""" @@ -163,30 +163,24 @@ def set_modelprms(self, modelprms): def _set_time_vals(self): """Set output dataset time and year values from dates_table.""" - if pygem_prms['climate']['gcm_wateryear'] == 'hydro': + if pygem_prms['climate']['sim_wateryear'] == 'hydro': self.year_type = 'water year' self.annual_columns = np.unique(self.dates_table['wateryear'].values)[ 0 : int(self.dates_table.shape[0] / 12) ] - elif pygem_prms['climate']['gcm_wateryear'] == 'calendar': + elif pygem_prms['climate']['sim_wateryear'] == 'calendar': self.year_type = 'calendar year' self.annual_columns = np.unique(self.dates_table['year'].values)[ 0 : int(self.dates_table.shape[0] / 12) ] - elif pygem_prms['climate']['gcm_wateryear'] == 'custom': + elif pygem_prms['climate']['sim_wateryear'] == 'custom': self.year_type = 'custom year' - self.time_values = self.dates_table.loc[ - pygem_prms['climate']['gcm_spinupyears'] * 12 : self.dates_table.shape[0] - + 1, - 'date', - ].tolist() + self.time_values = self.dates_table['date'].values.tolist() self.time_values = [ cftime.DatetimeNoLeap(x.year, x.month, x.day) for x in self.time_values ] # append additional year to self.year_values to account for mass and area at end of period - self.year_values = self.annual_columns[ - pygem_prms['climate']['gcm_spinupyears'] : self.annual_columns.shape[0] - ] + self.year_values = self.annual_columns self.year_values = np.concatenate( (self.year_values, np.array([self.annual_columns[-1] + 1])) ) @@ -198,11 +192,11 @@ def _model_params_record(self): # overwrite variables that are possibly different from pygem_input self.mdl_params_dict['ref_startyear'] = self.ref_startyear self.mdl_params_dict['ref_endyear'] = self.ref_endyear - self.mdl_params_dict['gcm_startyear'] = self.gcm_startyear - self.mdl_params_dict['gcm_endyear'] = self.gcm_endyear - self.mdl_params_dict['gcm_name'] = self.gcm_name + self.mdl_params_dict['sim_startyear'] = self.sim_startyear + self.mdl_params_dict['sim_endyear'] = self.sim_endyear + self.mdl_params_dict['sim_climate_name'] = self.sim_climate_name self.mdl_params_dict['realization'] = self.realization - self.mdl_params_dict['scenario'] = self.scenario + self.mdl_params_dict['sim_climate_scenario'] = self.sim_climate_scenario self.mdl_params_dict['option_calibration'] = self.option_calibration self.mdl_params_dict['option_bias_adjustment'] = self.option_bias_adjustment # record manually defined modelprms if calibration option is None @@ -370,9 +364,9 @@ def __post_init__(self): def _set_outdir(self): """Set the output directory path. Create if it does not already exist.""" - self.outdir += self.reg_str + '/' + self.gcm_name + '/' - if self.gcm_name not in ['ERA-Interim', 'ERA5', 'COAWST']: - self.outdir += self.scenario + '/' + self.outdir += self.reg_str + '/' + self.sim_climate_name + '/' + if self.sim_climate_name not in ['ERA-Interim', 'ERA5', 'COAWST']: + self.outdir += self.sim_climate_scenario + '/' self.outdir += 'stats/' # Create filepath if it does not exist os.makedirs(self.outdir, exist_ok=True) @@ -818,9 +812,9 @@ def __post_init__(self): def _set_outdir(self): """Set the output directory path. Create if it does not already exist.""" - self.outdir += self.reg_str + '/' + self.gcm_name + '/' - if self.gcm_name not in ['ERA-Interim', 'ERA5', 'COAWST']: - self.outdir += self.scenario + '/' + self.outdir += self.reg_str + '/' + self.sim_climate_name + '/' + if self.sim_climate_name not in ['ERA-Interim', 'ERA5', 'COAWST']: + self.outdir += self.sim_climate_scenario + '/' self.outdir += 'binned/' # Create filepath if it does not exist os.makedirs(self.outdir, exist_ok=True) diff --git a/pygem/pygem_modelsetup.py b/pygem/pygem_modelsetup.py index 8de10b07..9f25eed6 100755 --- a/pygem/pygem_modelsetup.py +++ b/pygem/pygem_modelsetup.py @@ -28,7 +28,7 @@ def datesmodelrun( startyear=pygem_prms['climate']['ref_startyear'], endyear=pygem_prms['climate']['ref_endyear'], - spinupyears=pygem_prms['climate']['ref_spinupyears'], + spinupyears=0, option_wateryear=pygem_prms['climate']['ref_wateryear'], ): """ diff --git a/pygem/setup/config.py b/pygem/setup/config.py index 11262932..33cdc1d0 100644 --- a/pygem/setup/config.py +++ b/pygem/setup/config.py @@ -146,19 +146,16 @@ def _validate_config(self, config): 'oggm.overwrite_gdirs': bool, 'oggm.has_internet': bool, 'climate': dict, - 'climate.ref_gcm_name': str, + 'climate.ref_climate_name': str, 'climate.ref_startyear': int, 'climate.ref_endyear': int, 'climate.ref_wateryear': str, - 'climate.ref_spinupyears': int, - 'climate.gcm_name': str, - 'climate.scenario': (str, type(None)), - 'climate.gcm_startyear': int, - 'climate.gcm_endyear': int, - 'climate.gcm_wateryear': str, + 'climate.sim_climate_name': str, + 'climate.sim_climate_scenario': (str, type(None)), + 'climate.sim_startyear': int, + 'climate.sim_endyear': int, + 'climate.sim_wateryear': str, 'climate.constantarea_years': int, - 'climate.gcm_spinupyears': int, - 'climate.hindcast': bool, 'climate.paths': dict, 'climate.paths.era5_relpath': str, 'climate.paths.era5_temp_fn': str, @@ -263,7 +260,7 @@ def _validate_config(self, config): 'calib.data.frontalablation.frontalablation_relpath': str, 'calib.data.frontalablation.frontalablation_cal_fn': str, 'calib.data.icethickness': dict, - 'calib.data.icethickness.h_consensus_relpath': str, + 'calib.data.icethickness.h_ref_relpath': str, 'calib.icethickness_cal_frac_byarea': float, 'sim': dict, 'sim.option_dynamics': (str, type(None)), @@ -301,9 +298,7 @@ def _validate_config(self, config): 'mb.include_debris': bool, 'mb.debris_relpath': str, 'mb.option_elev_ref_downscale': str, - 'mb.option_temp2bins': int, 'mb.option_adjusttemp_surfelev': int, - 'mb.option_prec2bins': int, 'mb.option_preclimit': int, 'mb.option_accumulation': int, 'mb.option_ablation': int, @@ -336,23 +331,16 @@ def _validate_config(self, config): 'time.wateryear_month_start': int, 'time.winter_month_start': int, 'time.summer_month_start': int, - 'time.option_dates': int, 'time.timestep': str, 'constants': dict, 'constants.density_ice': int, 'constants.density_water': int, - 'constants.area_ocean': float, 'constants.k_ice': float, 'constants.k_air': float, 'constants.ch_ice': int, 'constants.ch_air': int, 'constants.Lh_rf': int, 'constants.tolerance': float, - 'constants.gravity': float, - 'constants.pressure_std': int, - 'constants.temp_std': float, - 'constants.R_gas': float, - 'constants.molarmass_air': float, 'debug': dict, 'debug.refreeze': bool, 'debug.mb': bool, diff --git a/pygem/setup/config.yaml b/pygem/setup/config.yaml index b86044c6..57b6ae27 100644 --- a/pygem/setup/config.yaml +++ b/pygem/setup/config.yaml @@ -52,24 +52,20 @@ oggm: climate: # Reference period runs (reference period refers to the calibration period) # This will typically vary between 1980-present - ref_gcm_name: ERA5 # reference climate dataset + ref_climate_name: ERA5 # reference climate dataset ref_startyear: 2000 # first year of model run (reference dataset) ref_endyear: 2019 # last year of model run (reference dataset) ref_wateryear: calendar # options for years: 'calendar', 'hydro', 'custom' - ref_spinupyears: 0 # spin up years # GCM period used for simulation run - gcm_name: ERA5 # simulation climate dataset - scenario: null # simulation scenario - gcm_startyear: 2000 # first year of model run (simulation dataset) - gcm_endyear: 2019 # last year of model run (simulation dataset) - gcm_wateryear: calendar # options for years: 'calendar', 'hydro', 'custom' + sim_climate_name: ERA5 # simulation climate dataset + sim_climate_scenario: null # simulation scenario + sim_startyear: 2000 # first year of model run (simulation dataset) + sim_endyear: 2019 # last year of model run (simulation dataset) + sim_wateryear: calendar # options for years: 'calendar', 'hydro', 'custom' constantarea_years: 0 # number of years to not let the area or volume change - gcm_spinupyears: 0 # spin up years for simulation (output not set up for spinup years at present) - # Hindcast option (flips array so 1960-2000 would run 2000-1960 ensuring that glacier area at 2000 is correct) - hindcast: false # true: run hindcast simulation, false: do not # ===== CLIMATE DATA FILEPATHS AND FILENAMES ===== paths: - # ERA4 (default reference climate data) + # ERA5 (default reference climate data) era5_relpath: /climate_data/ERA5/ era5_temp_fn: ERA5_temp_monthly.nc era5_tempstd_fn: ERA5_tempstd_monthly.nc @@ -200,7 +196,7 @@ calib: frontalablation_cal_fn: all-frontalablation_cal_ind.csv # merged frontalablation calibration data (produced by run_calibration_frontalablation.py) # ice thickness icethickness: - h_consensus_relpath: /IceThickness_Farinotti/composite_thickness_RGI60-all_regions/ + h_ref_relpath: /IceThickness_Farinotti/composite_thickness_RGI60-all_regions/ icethickness_cal_frac_byarea: 0.9 # Regional glacier area fraction that is used to calibrate the ice thickness # e.g., 0.9 means only the largest 90% of glaciers by area will be used to calibrate @@ -228,11 +224,11 @@ sim: export_binned_area_threshold: 0 # Area threshold for exporting binned ice thickness # ===== OGGM DYNAMICS ===== oggm_dynamics: - cfl_number: 0.02 - cfl_number_calving: 0.01 - glena_reg_relpath: /Output/calibration/glena_region.csv + cfl_number: 0.02 # Time step threshold (seconds) + cfl_number_calving: 0.01 # Time step threshold for marine-terimating glaciers (seconds) use_reg_glena: true - # glens_a multiplier if not using regionally calibrated glens_a + glena_reg_relpath: /Output/calibration/glena_region.csv + # glen_a multiplier if not using regionally calibrated glens_a fs: 0 glen_a_multiplier: 1 # Mass redistribution / Glacier geometry change options @@ -269,10 +265,8 @@ mb: # Reference elevation options for downscaling climate variables option_elev_ref_downscale: Zmed # 'Zmed', 'Zmax', or 'Zmin' for median, maximum or minimum glacier elevations # Downscale temperature to bins options - option_temp2bins: 1 # 1: lr_gcm and lr_glac to adjust temp from gcm to the glacier bins option_adjusttemp_surfelev: 1 # 1: adjust temps based on surface elev changes; 0: no adjustment # Downscale precipitation to bins options - option_prec2bins: 1 # 1: prec_factor and prec_grad to adjust precip from gcm to the glacier bins option_preclimit: 0 # 1: limit the uppermost 25% using an expontial fxn # Accumulation model options @@ -330,26 +324,18 @@ time: wateryear_month_start: 10 # water year starting month winter_month_start: 10 # first month of winter (for HMA winter is October 1 - April 30) summer_month_start: 5 # first month of summer (for HMA summer is May 1 - Sept 30) - option_dates: 1 # 1: use dates from date table (first of each month), 2: dates from climate data timestep: monthly # time step ('monthly' only option at present) # ===== MODEL CONSTANTS ===== constants: density_ice: 900 # Density of ice [kg m-3] (or Gt / 1000 km3) density_water: 1000 # Density of water [kg m-3] - area_ocean: 362.5e+12 # Area of ocean [m2] (Cogley, 2012 from Marzeion et al. 2020) k_ice: 2.33 # Thermal conductivity of ice [J s-1 K-1 m-1] recall (W = J s-1) k_air: 0.023 # Thermal conductivity of air [J s-1 K-1 m-1] (Mellor, 1997) - #k_air: 0.001 # Thermal conductivity of air [J s-1 K-1 m-1] ch_ice: 1890000 # Volumetric heat capacity of ice [J K-1 m-3] (density=900, heat_capacity=2100 J K-1 kg-1) ch_air: 1297 # Volumetric Heat capacity of air [J K-1 m-3] (density=1.29, heat_capacity=1005 J K-1 kg-1) Lh_rf: 333550 # Latent heat of fusion [J kg-1] tolerance: 1.0e-12 # Model tolerance (used to remove low values caused by rounding errors) - gravity: 9.81 # Gravity [m s-2] - pressure_std: 101325 # Standard pressure [Pa] - temp_std: 288.15 # Standard temperature [K] - R_gas: 8.3144598 # Universal gas constant [J mol-1 K-1] - molarmass_air: 0.0289644 # Molar mass of Earth's air [kg mol-1] # ===== DEBUGGING OPTIONS ===== debug: diff --git a/pygem/shop/icethickness.py b/pygem/shop/icethickness.py index 9e1b9b7e..256bf526 100755 --- a/pygem/shop/icethickness.py +++ b/pygem/shop/icethickness.py @@ -43,7 +43,7 @@ @entity_task(log, writes=['consensus_mass']) def consensus_gridded( gdir, - h_consensus_fp=f'{pygem_prms["root"]}/{pygem_prms["calib"]["data"]["icethickness"]["h_consensus_relpath"]}', + h_consensus_fp=f'{pygem_prms["root"]}/{pygem_prms["calib"]["data"]["icethickness"]["h_ref_relpath"]}', add_mass=True, add_to_gridded=True, ): diff --git a/pygem/tests/test_notebooks.py b/pygem/tests/test_notebooks.py index 78c81cb8..d57d5b19 100644 --- a/pygem/tests/test_notebooks.py +++ b/pygem/tests/test_notebooks.py @@ -7,17 +7,23 @@ nb_dir = os.environ.get('PYGEM_NOTEBOOKS_DIRPATH') or os.path.join( os.path.expanduser('~'), 'PyGEM-notebooks' ) -notebooks = [f for f in os.listdir(nb_dir) if f.endswith('.ipynb')] +# TODO #54: Test all notebooks +# notebooks = [f for f in os.listdir(nb_dir) if f.endswith('.ipynb')] + +# list of notebooks to test, in the desired order +notebooks = [ + 'simple_test.ipynb', + 'advanced_test.ipynb', + 'advanced_test_tw.ipynb', +] @pytest.mark.parametrize('notebook', notebooks) def test_notebook(notebook): - # TODO #54: Test all notebooks + """ + Run pytest with nbmake on the specified notebook. - if notebook not in ( - 'simple_test.ipynb', - 'advanced_test.ipynb', - 'advanced_test_tw.ipynb', - ): - pytest.skip() + This test is parameterized to run each notebook individually, + preserving the order defined in the `notebooks` list. + """ subprocess.check_call(['pytest', '--nbmake', os.path.join(nb_dir, notebook)]) diff --git a/pygem/utils/_funcs.py b/pygem/utils/_funcs.py index a17dc127..1eb74738 100755 --- a/pygem/utils/_funcs.py +++ b/pygem/utils/_funcs.py @@ -68,6 +68,40 @@ def annualweightedmean_array(var, dates_table): return var_annual +def haversine_dist(grid_lons, grid_lats, target_lons, target_lats): + """ + Compute haversine distances between each (lon_target, lat_target) + and all (grid_lons, grid_lats) positions. + + Parameters: + - grid_lons: (ncol,) array of longitudes from data + - grid_lats: (ncol,) array of latitudes from data + - target_lons: (n_targets,) array of target longitudes + - target_lats: (n_targets,) array of target latitudes + + Returns: + - distances: (n_targets, ncol) array of distances in km to each grid location for each target + """ + R = 6371.0 # Earth radius in kilometers + + # Convert degrees to radians + grid_lons = np.radians(grid_lons)[np.newaxis, :] # (1, ncol) + grid_lats = np.radians(grid_lats)[np.newaxis, :] + target_lons = np.radians(target_lons)[:, np.newaxis] # (n_targets, 1) + target_lats = np.radians(target_lats)[:, np.newaxis] + + dlon = grid_lons - target_lons + dlat = grid_lats - target_lats + + a = ( + np.sin(dlat / 2.0) ** 2 + + np.cos(target_lats) * np.cos(grid_lats) * np.sin(dlon / 2.0) ** 2 + ) + c = 2 * np.arcsin(np.sqrt(a)) + + return R * c # (n_targets, ncol) + + def append_json(file_path, new_key, new_value): """ Opens a JSON file, reads its content, adds a new key-value pair,