Choose Your Own Basin¶
This tutorial provides an example of how to integrate RAPID2 into a workflow using Python to handle model setup and examination of model output. The code below allows for the selection of several model parameters at the start, such as basin, simulation date, and input runoff model, and then performs model setup automatically. Several python packages are used in this tutorial, but all are included in the RAPID2 pypi package except for Matplotlib, which is a commonly used Python plotting package.
This tutorial is provided as an interactive python notebook (ipynb), and can be opend and run in your browser by clicking the Binder badge:
Before starting, esnure that RAPID2 has been installed in your vistual environment (see the Quickstart Guide). If running this tutorial in Binder, this setup is done automatically. We can check that rapid2 is installed and verify the version:
First, ensure that the Python package for RAPID2 is installed.
%%bash
rapid2 --version && echo "[OK] - rapid2 is installed."
Next, import Python packages needed for the tutorial.
import os
import zipfile
from pathlib import Path
import urllib.request
import yaml
import netCDF4
import numpy as np
import matplotlib.pyplot as plt
from datetime import datetime
import getpass
import earthaccess
1. Set Parameters¶
Here is where we set the parameters for the run:
- BASIN: Set the basin ID based on global Pfafstetter Level 2 hydrologic regions (e.g. pfaf_74 for the Mississippi Basin)
- PHASE: Phase of GLDAS forcing to use (e.g. 2.1)
- MODEL: The hydrologic model to use for input (options are: VIC, CLSM, NOAH)
- TIME: The year and month for the simulation (e.g.'2010-01')
Note that all input paramaters must be provided as strings.
BASIN = "pfaf_31"
PHASE = '2.1'
MODEL = 'VIC'
TIME = '2010-01'
Next we set up the working directory, and input and output directories.
WORK_DIR = Path("rapid2_tutorial")
INPUT_DIR = WORK_DIR / "input"
OUTPUT_DIR = WORK_DIR / "output"
INPUT_DIR.mkdir(parents=True, exist_ok=True)
OUTPUT_DIR.mkdir(parents=True, exist_ok=True)
2. Donwload Static River Parameters¶
Next, we need to get six static parameter files for the Mississippi Basin. We'll be using prepared files from a Zenodo repository:
RAPID2 static files for global implementation with MERIT Basins v1.0 and GLDAS v2.0
This repository includes 61 zipped folders corresponding to the global Pfafstetter Level 2 hydrologic regions. The code below uses the pfaf ID set in the BASIN variable above to download and unpack the corresponding zipped folder.
ZENODO_RECORD = "20672740"
ZENODO_BASE = f"https://zenodo.org/records/{ZENODO_RECORD}/files"
print(f"Downloading {BASIN} ZIP file from Zenodo record {ZENODO_RECORD}...")
url = f"{ZENODO_BASE}/MERIT_Basins_v1.0_{BASIN}_GLDAS_v2.0_RAPID_v2.0.zip?download=1"
dest = INPUT_DIR / BASIN
if not dest.exists():
urllib.request.urlretrieve(url, dest)
else:
print(f'\tFiles already downloaded')
with zipfile.ZipFile(dest, "r") as z:
z.extractall(INPUT_DIR)
print(f"\nStatic files for {BASIN}:")
for f in sorted(INPUT_DIR.glob(f"*{BASIN}*")):
print(f" {f.name}")
3. Downlod NLDAS Runoff Data¶
The RAPID2 routing model requires runoff inputs from a land surface model. We'll use runoff data downloaded from the Global Land Data Assimilation System (GLDAS) using the dgldas2 command that is part of the rapid2 package.
A NASA Earthdata account is required to download the GLDAS data. If you have not yet done so, you may also need to accept an end-user license agreement from the NASA Goddard Earth Sciences (GES) Data and Information Services Center (DISC).
⚠️ If running this tutorial in Binder, you first need to authenticate with NASA Earthdata using an Earthdata token. Generate one at urs.earthdata.nasa.gov under "Generate Token", then run the cell below and paste it at the prompt. If running everything in your terminal, the
dgldas2command will handle your login credentials and the cell below can be skipped.
⚠️ Because Binder runs whatever code is in the repository it points at, launch this tutorial only from the badge on the official rapid-hub repository.
# Skip this cell if running in terminal
# Held in an environment variable, which the %%bash cells below inherit so that
# dgldas2 can authenticate.
os.environ["EARTHDATA_TOKEN"] = getpass.getpass("Earthdata token: ").strip()
if not os.environ["EARTHDATA_TOKEN"]:
raise RuntimeError(
"[ERROR] No token entered. Generate one at https://urs.earthdata.nasa.gov "
"under 'Generate Token', then re-run this cell."
)
auth = earthaccess.login(strategy="environment")
if not auth.authenticated:
raise RuntimeError("[ERROR] Earthdata login failed - re-run this cell.")
print("[OK] - Earthdata token registered for this session.")
GLDAS_FILE = INPUT_DIR / f'GLDAS_2.1_VIC_{TIME}.nc4'
!dgldas2 \
--phase {PHASE} \
--model {MODEL} \
--time {TIME} \
--land_surface_model \
{GLDAS_FILE}
The new file is placed in input/Tutorial. Other files are also automatically
downloaded in that directory and removed after use.
4. Coupling the Data with cpllsm¶
Once you have the raw land surface model data and the static river basin parameters, the two must be linked by mapping the runoff to the river network.. The cpllsm tool transforms the gridded runoff data into RAPID2's reach-based external inflow format. By providing the land surface model data, connectivity, coordinates, and coupling files, this tool generates the external inflow NetCDF file required for routing.
QEXT_FILE = INPUT_DIR / f"Qext_{BASIN}_GLDAS_2.1_VIC_2010-01.nc4"
!cpllsm \
--land_surface_model \
{GLDAS_FILE} \
--connectivity \
{INPUT_DIR}/con_{BASIN}.parquet \
--coordinates \
{INPUT_DIR}/crd_{BASIN}.parquet \
--coupling \
{INPUT_DIR}/cpl_{BASIN}_GLDAS.parquet \
--external_inflow \
{QEXT_FILE}
print(f"Wrote {QEXT_FILE}")
5. Create Cold Start File with zeroqinit¶
At the start of our simulation, we need to provide RAPID2 with initial values of discharge. We'll use the next command to create a cold start file with all zeroes.
QINIT_FILE = INPUT_DIR / f"Qinit_{BASIN}_GLDAS_2.1_VIC_2010-01.nc4"
!zeroqinit \
--external_inflow \
{QEXT_FILE} \
--initial_outflow \
{QINIT_FILE}
print(f"Created cold start file: {QINIT_FILE}")
6. Run the Routing Model with rapid2¶
RAPID2 uses a YAML configuration file (namelist) to specify input files, parameters, and output locations. IS_dtR is the routing time step (in seconds).
Before running the model, we need to generate the namelist file for our tutorial run. Here, we do this using the yaml python package to write the necessary input and output parameters to the namelist.
NL = {
'Qex_ncf': f'./{QEXT_FILE}',
'Q00_ncf': f'./{QINIT_FILE}',
'con_pqt': f'./{INPUT_DIR}/con_{BASIN}.parquet' ,
'kpr_pqt': f'./{INPUT_DIR}/kpr_{BASIN}_nrm.parquet',
'xpr_pqt': f'./{INPUT_DIR}/xpr_{BASIN}_nrm.parquet',
'bas_pqt': f'./{INPUT_DIR}/bas_{BASIN}_topo.parquet',
'IS_dtR': 900,
'Qou_ncf': f'./{OUTPUT_DIR}/Qout_{BASIN}_GLDAS_{PHASE}_{MODEL}_{TIME}_tst.nc4',
'Qfi_ncf': f'./{OUTPUT_DIR}/Qfinal_{BASIN}_GLDAS_{PHASE}_{MODEL}_{TIME}_tst.nc4'
}
with open(f'{INPUT_DIR}/name_list.yml', 'w') as yaml_file:
yaml.dump(NL,yaml_file,default_flow_style=False)
print(f"Wrote namelist file to: {INPUT_DIR}/name_list.yml")
Run the model using the rapid2 command.
!rapid2 --namelist {INPUT_DIR}/name_list.yml
7. Plot Output¶
Now we can examine the output dishcarge. In this tutorial, we use the netCDF4 python package to read our modeled output discharge file. This file contains time series of discharge for our selected month with time step set by the IS_dtR parameter in the namelist file. Here, we compute the mean discharge of every river reach in the basin, and then plot the discharge of reach with the highest mean.
with netCDF4.Dataset(f'./{OUTPUT_DIR}/Qout_{BASIN}_GLDAS_{PHASE}_{MODEL}_{TIME}_tst.nc4','r') as ds:
qout = ds.variables["Qout"][:]
rivid = ds.variables["rivid"][:]
time = ds.variables['time'][:]
date = []
#print(ds.variables)
for t in time:
date.append(datetime.fromtimestamp(t))
mean_per_reach = np.mean(qout, axis=0)
main = int(np.argmax(mean_per_reach))
peak_q=float(np.max(qout[:, main]))
fig, ax = plt.subplots(figsize=(15,5))
ax.set_ylim(0,peak_q*1.05)
ax.set_xlim(date[0],date[-1])
ax.set_ylabel(r'Discharge $[m^3/s]$')
ax.set_title(f'Outflow for {BASIN} Time:{TIME}')
plt.plot(date,qout[:,main])
#Run these lines to plot every reach on the same plot. This can be slow depending on the number of reaches in the basin.
#for i in range(len(rivid)):
# plt.plot(date,qout[:,i],alpha=0.2)