Lab 8-2: Modeling a snowpack with a gridded, numerical model - openAMUNDSEN#
From the documents:
openAMUNDSEN is a fully distributed model, designed primarily for resolving the mass and energy balance of snow and ice covered surfaces in mountain regions. Typically, it is applied in areas ranging from the point scale to the regional scale (i.e., up to some hundreds to thousands of square kilometers), using a spatial resolution of 10–100 m and a temporal resolution of 1–3 h, however its potential applications are very versatile.
You can learn more about the model here http://doc.openamundsen.org.
We run the model the entire upper East River Basin. We provide as inputs precipitation and meteorological (air temp, \(SW_{in}\), etc) conditions measured at Kettle Ponds as well a 50-meter resolution digital elevation model (DEM) of the upper East River Basin.
pip install openamundsen
Collecting openamundsen
Downloading openamundsen-1.1.6-py3-none-any.whl.metadata (15 kB)
Collecting cerberus (from openamundsen)
Downloading cerberus-1.3.8-py3-none-any.whl.metadata (5.5 kB)
Collecting loguru (from openamundsen)
Downloading loguru-0.7.3-py3-none-any.whl.metadata (22 kB)
Collecting munch (from openamundsen)
Downloading munch-4.0.0-py2.py3-none-any.whl.metadata (5.9 kB)
Requirement already satisfied: netCDF4 in /opt/hostedtoolcache/Python/3.11.14/x64/lib/python3.11/site-packages (from openamundsen) (1.7.4)
Collecting numba>=0.50.1 (from openamundsen)
Downloading numba-0.64.0-cp311-cp311-manylinux2014_x86_64.manylinux_2_17_x86_64.whl.metadata (2.9 kB)
Requirement already satisfied: numpy in /opt/hostedtoolcache/Python/3.11.14/x64/lib/python3.11/site-packages (from openamundsen) (2.4.2)
Requirement already satisfied: pandas>=1.1.0 in /opt/hostedtoolcache/Python/3.11.14/x64/lib/python3.11/site-packages (from openamundsen) (2.3.3)
Collecting pwlf (from openamundsen)
Downloading pwlf-2.5.2-py3-none-any.whl.metadata (6.4 kB)
Requirement already satisfied: pyproj in /opt/hostedtoolcache/Python/3.11.14/x64/lib/python3.11/site-packages (from openamundsen) (3.7.2)
Collecting setuptools_scm (from openamundsen)
Downloading setuptools_scm-9.2.2-py3-none-any.whl.metadata (7.7 kB)
Requirement already satisfied: setuptools>=61 in /opt/hostedtoolcache/Python/3.11.14/x64/lib/python3.11/site-packages (from openamundsen) (79.0.1)
Requirement already satisfied: scipy in /opt/hostedtoolcache/Python/3.11.14/x64/lib/python3.11/site-packages (from openamundsen) (1.17.1)
Collecting ruamel.yaml>=0.15.0 (from openamundsen)
Downloading ruamel_yaml-0.19.1-py3-none-any.whl.metadata (16 kB)
Requirement already satisfied: rasterio>=1.1.0 in /opt/hostedtoolcache/Python/3.11.14/x64/lib/python3.11/site-packages (from openamundsen) (1.4.4)
Requirement already satisfied: xarray>=0.14.0 in /opt/hostedtoolcache/Python/3.11.14/x64/lib/python3.11/site-packages (from openamundsen) (2026.2.0)
Collecting llvmlite<0.47,>=0.46.0dev0 (from numba>=0.50.1->openamundsen)
Downloading llvmlite-0.46.0-cp311-cp311-manylinux2014_x86_64.manylinux_2_17_x86_64.whl.metadata (5.0 kB)
Requirement already satisfied: python-dateutil>=2.8.2 in /opt/hostedtoolcache/Python/3.11.14/x64/lib/python3.11/site-packages (from pandas>=1.1.0->openamundsen) (2.9.0.post0)
Requirement already satisfied: pytz>=2020.1 in /opt/hostedtoolcache/Python/3.11.14/x64/lib/python3.11/site-packages (from pandas>=1.1.0->openamundsen) (2026.1.post1)
Requirement already satisfied: tzdata>=2022.7 in /opt/hostedtoolcache/Python/3.11.14/x64/lib/python3.11/site-packages (from pandas>=1.1.0->openamundsen) (2025.3)
Requirement already satisfied: six>=1.5 in /opt/hostedtoolcache/Python/3.11.14/x64/lib/python3.11/site-packages (from python-dateutil>=2.8.2->pandas>=1.1.0->openamundsen) (1.17.0)
Requirement already satisfied: affine in /opt/hostedtoolcache/Python/3.11.14/x64/lib/python3.11/site-packages (from rasterio>=1.1.0->openamundsen) (2.4.0)
Requirement already satisfied: attrs in /opt/hostedtoolcache/Python/3.11.14/x64/lib/python3.11/site-packages (from rasterio>=1.1.0->openamundsen) (25.4.0)
Requirement already satisfied: certifi in /opt/hostedtoolcache/Python/3.11.14/x64/lib/python3.11/site-packages (from rasterio>=1.1.0->openamundsen) (2026.2.25)
Requirement already satisfied: click!=8.2.*,>=4.0 in /opt/hostedtoolcache/Python/3.11.14/x64/lib/python3.11/site-packages (from rasterio>=1.1.0->openamundsen) (8.3.1)
Requirement already satisfied: cligj>=0.5 in /opt/hostedtoolcache/Python/3.11.14/x64/lib/python3.11/site-packages (from rasterio>=1.1.0->openamundsen) (0.7.2)
Requirement already satisfied: click-plugins in /opt/hostedtoolcache/Python/3.11.14/x64/lib/python3.11/site-packages (from rasterio>=1.1.0->openamundsen) (1.1.1.2)
Requirement already satisfied: pyparsing in /opt/hostedtoolcache/Python/3.11.14/x64/lib/python3.11/site-packages (from rasterio>=1.1.0->openamundsen) (3.3.2)
Requirement already satisfied: packaging>=24.1 in /opt/hostedtoolcache/Python/3.11.14/x64/lib/python3.11/site-packages (from xarray>=0.14.0->openamundsen) (26.0)
Requirement already satisfied: cftime in /opt/hostedtoolcache/Python/3.11.14/x64/lib/python3.11/site-packages (from netCDF4->openamundsen) (1.6.5)
Downloading openamundsen-1.1.6-py3-none-any.whl (133 kB)
Downloading numba-0.64.0-cp311-cp311-manylinux2014_x86_64.manylinux_2_17_x86_64.whl (3.7 MB)
?25l ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 0.0/3.7 MB ? eta -:--:--
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 3.7/3.7 MB 146.2 MB/s 0:00:00
?25hDownloading llvmlite-0.46.0-cp311-cp311-manylinux2014_x86_64.manylinux_2_17_x86_64.whl (56.3 MB)
?25l ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 0.0/56.3 MB ? eta -:--:--
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╸ 55.8/56.3 MB 279.1 MB/s eta 0:00:01
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 56.3/56.3 MB 208.1 MB/s 0:00:00
?25hDownloading ruamel_yaml-0.19.1-py3-none-any.whl (118 kB)
Downloading cerberus-1.3.8-py3-none-any.whl (30 kB)
Downloading loguru-0.7.3-py3-none-any.whl (61 kB)
Downloading munch-4.0.0-py2.py3-none-any.whl (9.9 kB)
Downloading pwlf-2.5.2-py3-none-any.whl (17 kB)
Downloading setuptools_scm-9.2.2-py3-none-any.whl (62 kB)
Installing collected packages: setuptools_scm, ruamel.yaml, munch, loguru, llvmlite, cerberus, pwlf, numba, openamundsen
?25l
━━━━╺━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 1/9 [ruamel.yaml]
━━━━━━━━━━━━━━━━━╸━━━━━━━━━━━━━━━━━━━━━━ 4/9 [llvmlite]
━━━━━━━━━━━━━━━━━╸━━━━━━━━━━━━━━━━━━━━━━ 4/9 [llvmlite]
━━━━━━━━━━━━━━━━━╸━━━━━━━━━━━━━━━━━━━━━━ 4/9 [llvmlite]
━━━━━━━━━━━━━━━━━╸━━━━━━━━━━━━━━━━━━━━━━ 4/9 [llvmlite]
━━━━━━━━━━━━━━━━━━━━━━╺━━━━━━━━━━━━━━━━━ 5/9 [cerberus]
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╺━━━━━━━━ 7/9 [numba]
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╺━━━━━━━━ 7/9 [numba]
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╺━━━━━━━━ 7/9 [numba]
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╺━━━━━━━━ 7/9 [numba]
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╺━━━━━━━━ 7/9 [numba]
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╺━━━━━━━━ 7/9 [numba]
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╺━━━━━━━━ 7/9 [numba]
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╺━━━━━━━━ 7/9 [numba]
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╺━━━━━━━━ 7/9 [numba]
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╺━━━━━━━━ 7/9 [numba]
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╺━━━━━━━━ 7/9 [numba]
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╺━━━━━━━━ 7/9 [numba]
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 9/9 [openamundsen]
?25h
Successfully installed cerberus-1.3.8 llvmlite-0.46.0 loguru-0.7.3 munch-4.0.0 numba-0.64.0 openamundsen-1.1.6 pwlf-2.5.2 ruamel.yaml-0.19.1 setuptools_scm-9.2.2
Note: you may need to restart the kernel to use updated packages.
import openamundsen as oa
import xarray as xr
import rioxarray as rix
import pandas as pd
from shapely.geometry import Point
import geopandas as gpd
from metpy.units import units
from metpy import constants
import altair as alt
import matplotlib.pyplot as plt
Run Model#
The model requires the following files to run:
a csv file of inputs (
1.csv),a csv of the locations of the input measurements (
stations.csv),a DEM of the model domain (
dem_uppereastriver_50.asc, also files with the same name ending with extensions.prjand.asc.aux.xml)a gridded mask indicating where to run the model (
roi_uppereastriver_50.asc|.prj) All of these files are within a compressed directory, which you can download form the Module 8 web page.
To get the model ready to run, download the compressed directory of model inputs by clicking the link openAMUNDSEN Inputs on the Module 8 webpage.
Once downloaded, unzip this directory and place it in a directory ../data/ such that there exists a directory ../data/openamundsen/.
Then, download the `openAMUNDSEN Configuration File`` (also linked on the Module 8 web page), and place it in the same directory as this notebook file.
Once you follow the steps above, the following command should show you a bunch of model input files:
ls ../data/openamundsen
1.csv dem_uppereastriver_50.prj stations.csv
dem_uppereastriver_50.asc roi_uppereastriver_50.asc
dem_uppereastriver_50.asc.aux.xml roi_uppereastriver_50.prj
Here is how we run the model.
read the config and set up a model object.
config = oa.read_config('open_amundsen_config.yml') # read in configuration file
model = oa.OpenAmundsen(config) # create OpenAmundsen object and populate unspecified parameters with default values
Initialize the model.
If there are any problems with your inputs, errors will be thrown here.
model.initialize() # read in input data files, initialize state variables etc.
2026-03-04 21:15:32 | INFO | Initializing model grid
2026-03-04 21:15:32 | INFO | Grid has dimensions 337x280
2026-03-04 21:15:32 | INFO | Reading DEM (../data/openamundsen/dem_uppereastriver_50.asc)
2026-03-04 21:15:32 | INFO | Calculating sky view factor
2026-03-04 21:15:32 | INFO | Calculating sky view factor for azimuth=0° (1/36)
2026-03-04 21:15:38 | INFO | Calculating sky view factor for azimuth=10° (2/36)
2026-03-04 21:15:38 | INFO | Calculating sky view factor for azimuth=20° (3/36)
2026-03-04 21:15:39 | INFO | Calculating sky view factor for azimuth=30° (4/36)
2026-03-04 21:15:39 | INFO | Calculating sky view factor for azimuth=40° (5/36)
2026-03-04 21:15:39 | INFO | Calculating sky view factor for azimuth=50° (6/36)
2026-03-04 21:15:39 | INFO | Calculating sky view factor for azimuth=60° (7/36)
2026-03-04 21:15:39 | INFO | Calculating sky view factor for azimuth=70° (8/36)
2026-03-04 21:15:39 | INFO | Calculating sky view factor for azimuth=80° (9/36)
2026-03-04 21:15:39 | INFO | Calculating sky view factor for azimuth=90° (10/36)
2026-03-04 21:15:40 | INFO | Calculating sky view factor for azimuth=100° (11/36)
2026-03-04 21:15:40 | INFO | Calculating sky view factor for azimuth=110° (12/36)
2026-03-04 21:15:40 | INFO | Calculating sky view factor for azimuth=120° (13/36)
2026-03-04 21:15:40 | INFO | Calculating sky view factor for azimuth=130° (14/36)
2026-03-04 21:15:40 | INFO | Calculating sky view factor for azimuth=140° (15/36)
2026-03-04 21:15:40 | INFO | Calculating sky view factor for azimuth=150° (16/36)
2026-03-04 21:15:40 | INFO | Calculating sky view factor for azimuth=160° (17/36)
2026-03-04 21:15:40 | INFO | Calculating sky view factor for azimuth=170° (18/36)
2026-03-04 21:15:41 | INFO | Calculating sky view factor for azimuth=180° (19/36)
2026-03-04 21:15:41 | INFO | Calculating sky view factor for azimuth=190° (20/36)
2026-03-04 21:15:41 | INFO | Calculating sky view factor for azimuth=200° (21/36)
2026-03-04 21:15:41 | INFO | Calculating sky view factor for azimuth=210° (22/36)
2026-03-04 21:15:41 | INFO | Calculating sky view factor for azimuth=220° (23/36)
2026-03-04 21:15:41 | INFO | Calculating sky view factor for azimuth=230° (24/36)
2026-03-04 21:15:41 | INFO | Calculating sky view factor for azimuth=240° (25/36)
2026-03-04 21:15:42 | INFO | Calculating sky view factor for azimuth=250° (26/36)
2026-03-04 21:15:42 | INFO | Calculating sky view factor for azimuth=260° (27/36)
2026-03-04 21:15:42 | INFO | Calculating sky view factor for azimuth=270° (28/36)
2026-03-04 21:15:42 | INFO | Calculating sky view factor for azimuth=280° (29/36)
2026-03-04 21:15:42 | INFO | Calculating sky view factor for azimuth=290° (30/36)
2026-03-04 21:15:42 | INFO | Calculating sky view factor for azimuth=300° (31/36)
2026-03-04 21:15:42 | INFO | Calculating sky view factor for azimuth=310° (32/36)
2026-03-04 21:15:42 | INFO | Calculating sky view factor for azimuth=320° (33/36)
2026-03-04 21:15:43 | INFO | Calculating sky view factor for azimuth=330° (34/36)
2026-03-04 21:15:43 | INFO | Calculating sky view factor for azimuth=340° (35/36)
2026-03-04 21:15:43 | INFO | Calculating sky view factor for azimuth=350° (36/36)
2026-03-04 21:15:43 | INFO | Reading ROI (../data/openamundsen/roi_uppereastriver_50.asc)
2026-03-04 21:15:43 | INFO | Reading meteo file: ../data/openamundsen/1.csv
2026-03-04 21:15:43 | INFO | Correcting station precipitation with method: wmo
2026-03-04 21:15:43 | INFO | Calculating terrain parameters
Run the model. This will take about two minutes with our current configuration.
model.run() # run the model
2026-03-04 21:15:43 | INFO | Starting model run
2026-03-04 21:15:43 | INFO | Processing time step 2022-11-16 00:00
2026-03-04 21:15:56 | INFO | Processing time step 2022-11-16 03:00
2026-03-04 21:15:56 | INFO | Processing time step 2022-11-16 06:00
2026-03-04 21:15:56 | INFO | Processing time step 2022-11-16 09:00
2026-03-04 21:15:56 | INFO | Processing time step 2022-11-16 12:00
2026-03-04 21:15:57 | INFO | Processing time step 2022-11-16 15:00
2026-03-04 21:15:57 | INFO | Processing time step 2022-11-16 18:00
2026-03-04 21:15:57 | INFO | Processing time step 2022-11-16 21:00
2026-03-04 21:15:57 | INFO | Processing time step 2022-11-17 00:00
2026-03-04 21:15:57 | INFO | Processing time step 2022-11-17 03:00
2026-03-04 21:15:57 | INFO | Processing time step 2022-11-17 06:00
2026-03-04 21:15:57 | INFO | Processing time step 2022-11-17 09:00
2026-03-04 21:15:57 | INFO | Processing time step 2022-11-17 12:00
2026-03-04 21:15:58 | INFO | Processing time step 2022-11-17 15:00
2026-03-04 21:15:58 | INFO | Processing time step 2022-11-17 18:00
2026-03-04 21:15:58 | INFO | Processing time step 2022-11-17 21:00
2026-03-04 21:15:58 | INFO | Processing time step 2022-11-18 00:00
2026-03-04 21:15:58 | INFO | Processing time step 2022-11-18 03:00
2026-03-04 21:15:58 | INFO | Processing time step 2022-11-18 06:00
2026-03-04 21:15:58 | INFO | Processing time step 2022-11-18 09:00
2026-03-04 21:15:58 | INFO | Processing time step 2022-11-18 12:00
2026-03-04 21:15:59 | INFO | Processing time step 2022-11-18 15:00
2026-03-04 21:15:59 | INFO | Processing time step 2022-11-18 18:00
2026-03-04 21:15:59 | INFO | Processing time step 2022-11-18 21:00
2026-03-04 21:15:59 | INFO | Processing time step 2022-11-19 00:00
2026-03-04 21:15:59 | INFO | Processing time step 2022-11-19 03:00
2026-03-04 21:15:59 | INFO | Processing time step 2022-11-19 06:00
2026-03-04 21:15:59 | INFO | Processing time step 2022-11-19 09:00
2026-03-04 21:15:59 | INFO | Processing time step 2022-11-19 12:00
2026-03-04 21:16:00 | INFO | Processing time step 2022-11-19 15:00
2026-03-04 21:16:00 | INFO | Processing time step 2022-11-19 18:00
2026-03-04 21:16:00 | INFO | Processing time step 2022-11-19 21:00
2026-03-04 21:16:00 | INFO | Processing time step 2022-11-20 00:00
2026-03-04 21:16:00 | INFO | Processing time step 2022-11-20 03:00
2026-03-04 21:16:00 | INFO | Processing time step 2022-11-20 06:00
2026-03-04 21:16:00 | INFO | Processing time step 2022-11-20 09:00
2026-03-04 21:16:00 | INFO | Processing time step 2022-11-20 12:00
2026-03-04 21:16:01 | INFO | Processing time step 2022-11-20 15:00
2026-03-04 21:16:01 | INFO | Processing time step 2022-11-20 18:00
2026-03-04 21:16:01 | INFO | Processing time step 2022-11-20 21:00
2026-03-04 21:16:01 | INFO | Processing time step 2022-11-21 00:00
2026-03-04 21:16:01 | INFO | Processing time step 2022-11-21 03:00
2026-03-04 21:16:01 | INFO | Processing time step 2022-11-21 06:00
2026-03-04 21:16:01 | INFO | Processing time step 2022-11-21 09:00
2026-03-04 21:16:01 | INFO | Processing time step 2022-11-21 12:00
2026-03-04 21:16:01 | INFO | Processing time step 2022-11-21 15:00
2026-03-04 21:16:02 | INFO | Processing time step 2022-11-21 18:00
2026-03-04 21:16:02 | INFO | Processing time step 2022-11-21 21:00
2026-03-04 21:16:02 | INFO | Processing time step 2022-11-22 00:00
2026-03-04 21:16:02 | INFO | Processing time step 2022-11-22 03:00
2026-03-04 21:16:02 | INFO | Processing time step 2022-11-22 06:00
2026-03-04 21:16:02 | INFO | Processing time step 2022-11-22 09:00
2026-03-04 21:16:02 | INFO | Processing time step 2022-11-22 12:00
2026-03-04 21:16:02 | INFO | Processing time step 2022-11-22 15:00
2026-03-04 21:16:03 | INFO | Processing time step 2022-11-22 18:00
2026-03-04 21:16:03 | INFO | Processing time step 2022-11-22 21:00
2026-03-04 21:16:03 | INFO | Processing time step 2022-11-23 00:00
2026-03-04 21:16:03 | INFO | Processing time step 2022-11-23 03:00
2026-03-04 21:16:03 | INFO | Processing time step 2022-11-23 06:00
2026-03-04 21:16:03 | INFO | Processing time step 2022-11-23 09:00
2026-03-04 21:16:03 | INFO | Processing time step 2022-11-23 12:00
2026-03-04 21:16:03 | INFO | Processing time step 2022-11-23 15:00
2026-03-04 21:16:04 | INFO | Processing time step 2022-11-23 18:00
2026-03-04 21:16:04 | INFO | Processing time step 2022-11-23 21:00
2026-03-04 21:16:04 | INFO | Processing time step 2022-11-24 00:00
2026-03-04 21:16:04 | INFO | Processing time step 2022-11-24 03:00
2026-03-04 21:16:04 | INFO | Processing time step 2022-11-24 06:00
2026-03-04 21:16:04 | INFO | Processing time step 2022-11-24 09:00
2026-03-04 21:16:04 | INFO | Processing time step 2022-11-24 12:00
2026-03-04 21:16:04 | INFO | Processing time step 2022-11-24 15:00
2026-03-04 21:16:05 | INFO | Processing time step 2022-11-24 18:00
2026-03-04 21:16:05 | INFO | Processing time step 2022-11-24 21:00
2026-03-04 21:16:05 | INFO | Processing time step 2022-11-25 00:00
2026-03-04 21:16:05 | INFO | Processing time step 2022-11-25 03:00
2026-03-04 21:16:05 | INFO | Processing time step 2022-11-25 06:00
2026-03-04 21:16:05 | INFO | Processing time step 2022-11-25 09:00
2026-03-04 21:16:05 | INFO | Processing time step 2022-11-25 12:00
2026-03-04 21:16:05 | INFO | Processing time step 2022-11-25 15:00
2026-03-04 21:16:06 | INFO | Processing time step 2022-11-25 18:00
2026-03-04 21:16:06 | INFO | Processing time step 2022-11-25 21:00
2026-03-04 21:16:06 | INFO | Processing time step 2022-11-26 00:00
2026-03-04 21:16:06 | INFO | Processing time step 2022-11-26 03:00
2026-03-04 21:16:06 | INFO | Processing time step 2022-11-26 06:00
2026-03-04 21:16:06 | INFO | Processing time step 2022-11-26 09:00
2026-03-04 21:16:06 | INFO | Processing time step 2022-11-26 12:00
2026-03-04 21:16:06 | INFO | Processing time step 2022-11-26 15:00
2026-03-04 21:16:07 | INFO | Processing time step 2022-11-26 18:00
2026-03-04 21:16:07 | INFO | Processing time step 2022-11-26 21:00
2026-03-04 21:16:07 | INFO | Processing time step 2022-11-27 00:00
2026-03-04 21:16:07 | INFO | Processing time step 2022-11-27 03:00
2026-03-04 21:16:07 | INFO | Processing time step 2022-11-27 06:00
2026-03-04 21:16:07 | INFO | Processing time step 2022-11-27 09:00
2026-03-04 21:16:07 | INFO | Processing time step 2022-11-27 12:00
2026-03-04 21:16:07 | INFO | Processing time step 2022-11-27 15:00
2026-03-04 21:16:08 | INFO | Processing time step 2022-11-27 18:00
2026-03-04 21:16:08 | INFO | Processing time step 2022-11-27 21:00
2026-03-04 21:16:08 | INFO | Processing time step 2022-11-28 00:00
2026-03-04 21:16:08 | INFO | Processing time step 2022-11-28 03:00
2026-03-04 21:16:08 | INFO | Processing time step 2022-11-28 06:00
2026-03-04 21:16:08 | INFO | Processing time step 2022-11-28 09:00
2026-03-04 21:16:08 | INFO | Processing time step 2022-11-28 12:00
2026-03-04 21:16:08 | INFO | Processing time step 2022-11-28 15:00
2026-03-04 21:16:09 | INFO | Processing time step 2022-11-28 18:00
2026-03-04 21:16:09 | INFO | Processing time step 2022-11-28 21:00
2026-03-04 21:16:09 | INFO | Processing time step 2022-11-29 00:00
2026-03-04 21:16:09 | INFO | Processing time step 2022-11-29 03:00
2026-03-04 21:16:09 | INFO | Processing time step 2022-11-29 06:00
2026-03-04 21:16:09 | INFO | Processing time step 2022-11-29 09:00
2026-03-04 21:16:09 | INFO | Processing time step 2022-11-29 12:00
2026-03-04 21:16:09 | INFO | Processing time step 2022-11-29 15:00
2026-03-04 21:16:10 | INFO | Processing time step 2022-11-29 18:00
2026-03-04 21:16:10 | INFO | Processing time step 2022-11-29 21:00
2026-03-04 21:16:10 | INFO | Processing time step 2022-11-30 00:00
2026-03-04 21:16:10 | INFO | Processing time step 2022-11-30 03:00
2026-03-04 21:16:10 | INFO | Processing time step 2022-11-30 06:00
2026-03-04 21:16:10 | INFO | Processing time step 2022-11-30 09:00
2026-03-04 21:16:10 | INFO | Processing time step 2022-11-30 12:00
2026-03-04 21:16:11 | INFO | Processing time step 2022-11-30 15:00
2026-03-04 21:16:11 | INFO | Processing time step 2022-11-30 18:00
2026-03-04 21:16:11 | INFO | Processing time step 2022-11-30 21:00
2026-03-04 21:16:11 | INFO | Processing time step 2022-12-01 00:00
2026-03-04 21:16:11 | INFO | Processing time step 2022-12-01 03:00
2026-03-04 21:16:11 | INFO | Processing time step 2022-12-01 06:00
2026-03-04 21:16:11 | INFO | Processing time step 2022-12-01 09:00
2026-03-04 21:16:11 | INFO | Processing time step 2022-12-01 12:00
2026-03-04 21:16:12 | INFO | Processing time step 2022-12-01 15:00
2026-03-04 21:16:12 | INFO | Processing time step 2022-12-01 18:00
2026-03-04 21:16:12 | INFO | Processing time step 2022-12-01 21:00
2026-03-04 21:16:12 | INFO | Processing time step 2022-12-02 00:00
2026-03-04 21:16:12 | INFO | Processing time step 2022-12-02 03:00
2026-03-04 21:16:12 | INFO | Processing time step 2022-12-02 06:00
2026-03-04 21:16:12 | INFO | Processing time step 2022-12-02 09:00
2026-03-04 21:16:13 | INFO | Processing time step 2022-12-02 12:00
2026-03-04 21:16:13 | INFO | Processing time step 2022-12-02 15:00
2026-03-04 21:16:13 | INFO | Processing time step 2022-12-02 18:00
2026-03-04 21:16:13 | INFO | Processing time step 2022-12-02 21:00
2026-03-04 21:16:13 | INFO | Processing time step 2022-12-03 00:00
2026-03-04 21:16:13 | INFO | Processing time step 2022-12-03 03:00
---------------------------------------------------------------------------
KeyboardInterrupt Traceback (most recent call last)
Cell In[6], line 1
----> 1 model.run() # run the model
File /opt/hostedtoolcache/Python/3.11.14/x64/lib/python3.11/site-packages/openamundsen/model.py:216, in OpenAmundsen.run(self)
213 start_time = time.time()
215 for _ in range(len(self.dates)):
--> 216 self.run_single()
218 time_diff = pd.Timedelta(seconds=(time.time() - start_time))
219 logger.success("Model run finished. Runtime: " + str(time_diff))
File /opt/hostedtoolcache/Python/3.11.14/x64/lib/python3.11/site-packages/openamundsen/model.py:252, in OpenAmundsen.run_single(self)
249 self._meteo_callback(self)
251 self._process_meteo_data()
--> 252 self._model_interface()
253 self.point_output.update()
254 self.gridded_output.update()
File /opt/hostedtoolcache/Python/3.11.14/x64/lib/python3.11/site-packages/openamundsen/model.py:316, in OpenAmundsen._model_interface(self)
311 def _model_interface(self):
312 """
313 Interface for calling the different submodules. This method is called
314 in every time step after the meteorological fields have been prepared.
315 """
--> 316 self._irradiance()
318 if self._require_land_cover:
319 self.land_cover.lai()
File /opt/hostedtoolcache/Python/3.11.14/x64/lib/python3.11/site-packages/openamundsen/model.py:823, in OpenAmundsen._irradiance(self)
820 roi = self.grid.roi
821 m.sw_in[roi] = m.sw_in_clearsky[roi] * m.cloud_factor[roi]
--> 823 modules.radiation.longwave_irradiance(self)
File /opt/hostedtoolcache/Python/3.11.14/x64/lib/python3.11/site-packages/openamundsen/modules/radiation/irradiance.py:298, in longwave_irradiance(model)
291 rock_emission_factor = 0.01 # (K W-1 m2) temperature of emitting rocks during daytime is assumed to be higher than the air temperature by this factor multiplied by the incoming shortwave radiation (Greuell et al., 1997) # noqa: E501
293 # Incoming longwave radiation from the clear sky
294 lw_in_clearsky = (
295 clear_sky_emissivity
296 * constants.STEFAN_BOLTZMANN
297 * m.temp[roi] ** 4
--> 298 * model.state.base.svf[roi]
299 * (1 - m.cloud_fraction[roi] ** 2)
300 )
302 # Incoming longwave radiation from clouds
303 lw_in_clouds = (
304 cloud_emissivity
305 * constants.STEFAN_BOLTZMANN
(...) 308 * m.cloud_fraction[roi] ** 2
309 )
File /opt/hostedtoolcache/Python/3.11.14/x64/lib/python3.11/site-packages/munch/__init__.py:91, in Munch.__getattr__(self, k)
88 self.update(*args, **kwargs)
90 # only called if k not found in normal places
---> 91 def __getattr__(self, k):
92 """ Gets key if it exists, otherwise throws AttributeError.
93
94 nb. __getattr__ is only called if key is not found in normal places.
(...) 112 True
113 """
114 try:
115 # Throws exception if not in prototype chain
KeyboardInterrupt:
Examine results#
Now, let’s open up the model results and check them out! There are two types of model result output files.
Point data results are results at Kettle Ponds. We configured the model (in the yml file) to save results at the Kettle Ponds location. You can add additional points at which to save model results. Point results have the full suite of model outputs, for each 3-hour timestep. By default, all outputs are saved for point data results.
Gridded data results are results for the entire watershed, i.e. modeled values at each grid cell, where grid cells are 4. By default, limited gridded model results are saved to file because the size of the outputs are large. We configured the model (in the yml file) to only save the snow melt variable. If you are interested in checking out other gridded results, look in the configuration file, under the section grids: variables:, and add variables there, and then rerun the model.
We will look at model results and compare them to measurements from the Kettle Ponds site. So let’s open up the measured SOS data.
measured = xr.open_dataset('../data/sos_full_dataset_30min.nc')
Point Results#
# We grab the point `kps` right off the bat b/c it's the only point with data inside this dataset.
modeled = xr.open_dataset("openamundsen_results/output_timeseries.nc").sel(point='kps')
Resample the 30-minute measurements to match the 3-hour model results. Resampling using a pandas dataframe is way faster than using an xarray dataset, so we convert to dataframe, perform the resampling, then convert back to an xarray dataset. We also trim the measurements to match the time period of the model.
measured = measured.to_dataframe().resample('3h').mean().to_xarray()
measured = measured.sel(time = slice(modeled.time.min(), modeled.time.max()))
Compare measured and modeled SWE#
measured['SWE_p3_c'].sel(time=slice('20221110', '20230630')).plot(label='Measured tow', color='orange')
modeled.swe.sel(time=slice('20221110', '20230630')).plot(label='Modeled')
plt.legend()
plt.show()
We see that modeled SWE divergences from meausred SWE at the beginning of April. Let’s examine the measured and modeled energy balances during the first half of April.
Calculate energy balance terms#
We will focus on the following energy balance terms: \(R_{net}\), \(H_L\), \(H_s\), and \(P\), where \(P\) is the precipitation energy flux. We found in previous labs that the ground heat flux \(G\) and internal energy \(dU/dt\) terms are generally small, so we omit them here for simplicity. \(P\) is generally small because the snowfall is usually pretty close in temperature to the snowpack, but for the homework assignment, we will model a rain-on-snow event, where \(P\) can be large (i.e. rain is warmer than the snowpack, so rainfall results in warming of the snowpack).
Modeled#
# Net Radiation
modeled['Rnet'] = modeled['sw_in'] - modeled['sw_out'] + modeled['lw_in'] - modeled['lw_out']
# Precipitation energy flux
temp_difference = (
modeled.temp - modeled.surface_temp
) * units("kelvin")
precip_per_3hours = modeled.precip * units("kg/m^2")
seconds_in_3hours = 10800 * units('seconds')
precip_rate = precip_per_3hours / seconds_in_3hours
modeled['P'] = (precip_rate * constants.water_specific_heat * temp_difference)
# Note that HL and Hs, in the modeled modeled, are already in units of W/m^2
print(modeled['sens_heat_flux'].attrs['units'])
print(modeled['lat_heat_flux'].attrs['units'])
Measured#
measured['Rnet'] = measured['Rsw_in_9m_d'] - measured['Rsw_out_9m_d'] + measured['Rlw_in_9m_d'] - measured['Rlw_out_9m_d']
# Precipitation energy flux, use modeled precip b/c that's an input, and we don't have the precip dataset loaded, so this is convenient
temp_difference = (
measured['T_3m_c'] - measured['Tsurf_c']
) * units("kelvin")
precip_per_3hours = modeled.precip * units("kg/m^2")
seconds_in_3hours = 10800 * units('seconds')
precip_rate = precip_per_3hours / seconds_in_3hours
measured['P'] = (precip_rate * constants.water_specific_heat * temp_difference)
# convert measured sensible and latent heat fluxes to W/m^2. See Lab5-1.
# Note that while the modeled sensible and latent heat fluxes are **into** the snowpack,
# our measurements of heat flux are **upwards/away** from the snowpack. So we negate
# the measured values to match the modeled values.
latent_heat_sublimation = 2590 #J/g
measured['lat_heat_flux'] = - 1 * measured['w_h2o__3m_c'] * latent_heat_sublimation
specific_heat_capacity_air = 1.005 # J/K/g
air_density = 1000 # g/m^3
measured['sens_heat_flux'] = - 1 * measured['w_tc__3m_c'] * specific_heat_capacity_air * air_density
Compare#
fig, axes = plt.subplots(1,2, figsize=((12.8, 4.8)), sharex=True, sharey=True)
time_slice = slice('20230401', '20230415')
modeled['Rnet'].sel(time=time_slice).plot(label = 'Rnet', ax=axes[0])
modeled['sens_heat_flux'].sel(time=time_slice).plot(label = 'Hs', ax=axes[0])
modeled['lat_heat_flux'].sel(time=time_slice).plot(label = 'HL', ax=axes[0])
modeled['P'].sel(time=time_slice).plot(label = 'P', ax=axes[0])
measured['Rnet'].sel(time=time_slice).plot(label = 'Rnet', ax=axes[1])
measured['sens_heat_flux'].sel(time=time_slice).plot(label = 'Hs', ax=axes[1])
measured['lat_heat_flux'].sel(time=time_slice).plot(label = 'HL', ax=axes[1])
measured['P'].sel(time=time_slice).plot(label = 'P', ax=axes[1])
axes[0].set_title('Modeled'); axes[0].legend()
axes[1].set_title('Measured'); axes[1].legend()
plt.show()
Ok, so we see some interesting things here. Most obviously, modeled net radiation is smaller than measured net radiation, particularly starting on April 8th. What do we recall from previous labs about what happened around this date that the model doesn’t pick up on?
Second, it appears that the latent and sensible heat fluxes are larger in the model-world than in the real world.
Gridded results#
modeled_gridded = xr.open_dataset("openamundsen_results/output_grids.nc")
Let’s check out melt rates at the four dates throughout May and June. The plot of modeled SWE (above) shows that modeled melt at Kettle Ponds occurs all throughout May, but higher altitude locations probably melted out later, in June.
fig, axes = plt.subplots(2, 2, figsize=(10, 8), sharex=True, sharey=True)
axes = axes.flatten()
# Plot for May 1
modeled_gridded['melt'].sel(time='2023-05-01T00:00:00.000000000').plot(ax=axes[0], cmap='viridis', vmin=0, vmax=1.5)
axes[0].set_title('May 1 Melt Rate')
# Plot for May 15
modeled_gridded['melt'].sel(time='2023-05-15T00:00:00.000000000').plot(ax=axes[1], cmap='viridis', vmin=0, vmax=1.5)
axes[1].set_title('May 15 Melt Rate')
# Plot for June 2
modeled_gridded['melt'].sel(time='2023-06-02T00:00:00.000000000').plot(ax=axes[2], cmap='viridis', vmin=0, vmax=1.5)
axes[2].set_title('June 2 Melt Rate')
# Plot for June 15
modeled_gridded['melt'].sel(time='2023-06-14T00:00:00.000000000').plot(ax=axes[3], cmap='viridis', vmin=0, vmax=1.5)
axes[3].set_title('June 15 Melt Rate')
kps_loc = gpd.GeoSeries([Point(-106.972983, 38.941817, 0)],
index = ['Kettle Ponds'],
crs = 'epsg:4326'
).to_crs(modeled_gridded.rio.crs)
for ax in axes.flatten():
# Set all the plots to have the same x and y-axis limits, and set the aspect ratio to be equal
ax.set_xlim(322500, 337000)
ax.set_ylim(4305000, 4322500)
ax.set_aspect('equal')
# plot the KPS location with a red dot
kps_loc.plot(ax=ax, color='red')
plt.tight_layout()
plt.show()
Cool! So we see on May 1 that melt was beginning at the lowest elevation parts of the basin (see the DEM for the upper East River Valley below).
We also see that by May 15, melt was creeping up in elevation, mostly occuring in the low elevation, valley center.
By June 2, the lowest parts of the valley had melted out completely and so no more melt could occur. Surrounding the lowest elevation areas, melt was most intense.
By June 15, even more low parts of the valley had melted out, and a higher elevation band was melting out. These plots also show that for higher elevation areas in the basin, not much melt had ocurred by June 15th. Unfortunately, our Kettle Ponds measurements did not go past June 19, and we only ran the model till June 15th. To model the full snowpack at higher elevations, we would need to run the model up until a later date.
dem = rix.open_rasterio("../data/openamundsen/dem_uppereastriver_50.asc")
dem.plot(cbar_kwargs={'label': "Elevation (m)"})
plt.show()