{
"cells": [
{
"cell_type": "markdown",
"id": "1348f884",
"metadata": {},
"source": [
"# Jupyter Notebook for Validation of Results\n",
"\n",
"## Overview\n",
"\n",
"This notebook compares the results obtained from the two implementations, `shbundle` and `pyshbundle`, for a specific task.
\n",
"We will focus on evaluating the following:\n",
"1. Root Mean Squared Error (RMSE)\n",
"2. Normalized Root Mean Squared Error (NRMSE)\n",
"3. Difference in global area weighted water budget closure.\n",
"4. The difference in time series for a given basin.\n",
"\n",
"## Evaluation Metrics\n",
"\n",
"### 1. Root Mean Squared Error (RMSE)\n",
"\n",
"The RMSE is a measure of the average magnitude of the errors between predicted and observed values.
It is calculated as follows:\n",
"\n",
"$$ RMSE = \\sqrt{\\frac{1}{n}\\Sigma_{i=1}^{n}{\\Big({y_i - \\bar{y}}\\Big)^2}}$$\n",
"\n",
"where:\n",
"- $n$ is the number of observations\n",
"- ${y_i}$ is the $i_{th}$ observation\n",
"- $\\bar{y}$ is the mean of observations\n",
"\n",
"### 2. Normalized Root Mean Squared Error (NRMSE)\n",
"\n",
"NRMSE is a normalized version of RMSE, which provides a relative measure of the error compared to the range of the observed values. It is calculated as:\n",
"\n",
"$$ NRMSE = \\frac{RMSE}{\\max(y) - \\min(y)}$$\n",
"\n",
"### 3. Difference in global area weighted water budget closure.\n",
"\n",
"Since the total mass of Earth is conserved (at the accuracy we can currently observe), the mass change observed over any month should sum up to zero. To verify this, we area weight the solutions from `shbundle` & `pyshbundle` using the area_weighting() function by creating a global area grid. Then we multiply this with the global solutions at each time step and take the average, followed by plotting the results.\n",
"\n",
"### 4. Difference in basin-average Time Series\n",
"\n",
"Finally we create a basin-average TWSA time-series from both the implementations."
]
},
{
"cell_type": "markdown",
"id": "61ad9a48",
"metadata": {},
"source": [
"## Data Preparation\n",
"\n",
"Before proceeding with the analysis, ensure that the data from both `shbundle` and `pyshbundle` are loaded into the notebook."
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "b244d043-7498-4e8a-a638-b9d0703a5755",
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import pandas as pd\n",
"import xarray as xr\n",
"import matplotlib.pyplot as plt\n",
"import os"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "a07fc51a-4640-4f03-a2e2-d65fc8da751b",
"metadata": {
"tags": []
},
"outputs": [
{
"data": {
"text/plain": [
"array([[[ -2.5984123, -2.5491967, -2.4989052, ..., -2.7395167,\n",
" -2.693581 , -2.6465435],\n",
" [ -8.085175 , -7.9737053, -7.859015 , ..., -8.400428 ,\n",
" -8.298518 , -8.1934395],\n",
" [ -12.36135 , -12.230325 , -12.094314 , ..., -12.725532 ,\n",
" -12.608853 , -12.487493 ],\n",
" ...,\n",
" [ -63.64362 , -64.29353 , -64.944595 , ..., -61.70406 ,\n",
" -62.34856 , -62.995186 ],\n",
" [ -59.735207 , -60.10004 , -60.464703 , ..., -58.641277 ,\n",
" -59.00567 , -59.37036 ],\n",
" [ -56.28362 , -56.395493 , -56.50705 , ..., -55.946453 ,\n",
" -56.05906 , -56.17146 ]],\n",
"\n",
" [[ -1.6946957, -1.6432129, -1.5906656, ..., -1.8426621,\n",
" -1.7944294, -1.7451043],\n",
" [ -7.195656 , -7.077713 , -6.9565587, ..., -7.530373 ,\n",
" -7.42197 , -7.310403 ],\n",
" [ -11.6058655, -11.462069 , -11.313266 , ..., -12.008124 ,\n",
" -11.878803 , -11.744747 ],\n",
" ...,\n",
" [ -60.940014 , -61.585247 , -62.23357 , ..., -59.026295 ,\n",
" -59.660206 , -60.29822 ],\n",
" [ -58.615997 , -58.973927 , -59.332268 , ..., -57.546318 ,\n",
" -57.902035 , -58.258648 ],\n",
" [ -56.071903 , -56.180115 , -56.288097 , ..., -55.7463 ,\n",
" -55.854958 , -55.963505 ]],\n",
"\n",
" [[ 28.409975 , 28.435665 , 28.462093 , ..., 28.337358 ,\n",
" 28.360819 , 28.385025 ],\n",
" [ 24.354195 , 24.405733 , 24.459566 , ..., 24.213146 ,\n",
" 24.257923 , 24.304932 ],\n",
" [ 20.799301 , 20.849913 , 20.904062 , ..., 20.667936 ,\n",
" 20.708385 , 20.75215 ],\n",
" ...,\n",
" [ -53.207798 , -53.74962 , -54.29029 , ..., -51.579227 ,\n",
" -52.12223 , -52.66521 ],\n",
" [ -48.65204 , -48.940662 , -49.228447 , ..., -47.782604 ,\n",
" -48.07286 , -48.36272 ],\n",
" [ -45.042896 , -45.128735 , -45.214226 , ..., -44.783646 ,\n",
" -44.870323 , -44.956753 ]],\n",
"\n",
" ...,\n",
"\n",
" [[ -91.37819 , -91.51497 , -91.65443 , ..., -90.984245 ,\n",
" -91.1128 , -91.244125 ],\n",
" [ -78.2787 , -78.60071 , -78.93178 , ..., -77.366776 ,\n",
" -77.66176 , -77.96573 ],\n",
" [ -67.629616 , -68.0414 , -68.46865 , ..., -66.48567 ,\n",
" -66.85189 , -67.23316 ],\n",
" ...,\n",
" [ 169.38713 , 170.59637 , 171.79698 , ..., 165.71684 ,\n",
" 166.94644 , 168.17017 ],\n",
" [ 155.91563 , 156.5555 , 157.19012 , ..., 153.96843 ,\n",
" 154.62172 , 155.27092 ],\n",
" [ 143.16539 , 143.34567 , 143.52426 , ..., 142.61514 ,\n",
" 142.80006 , 142.9835 ]],\n",
"\n",
" [[-125.82837 , -125.960236 , -126.09478 , ..., -125.44908 ,\n",
" -125.57277 , -125.6992 ],\n",
" [-112.4685 , -112.77872 , -113.098 , ..., -111.592 ,\n",
" -111.87517 , -112.16732 ],\n",
" [-101.242386 , -101.638016 , -102.049255 , ..., -100.147675 ,\n",
" -100.49736 , -100.86221 ],\n",
" ...,\n",
" [ 149.24057 , 150.4772 , 151.70389 , ..., 145.48059 ,\n",
" 146.7413 , 147.99495 ],\n",
" [ 134.43509 , 135.08853 , 135.73586 , ..., 132.44183 ,\n",
" 133.11136 , 133.7759 ],\n",
" [ 120.18574 , 120.36971 , 120.551735 , ..., 119.62282 ,\n",
" 119.81222 , 119.999886 ]],\n",
"\n",
" [[-101.70098 , -101.829636 , -101.96099 , ..., -101.33151 ,\n",
" -101.4519 , -101.575066 ],\n",
" [ -88.18868 , -88.48786 , -88.79633 , ..., -87.346596 ,\n",
" -87.61808 , -87.898766 ],\n",
" [ -76.74753 , -77.12447 , -77.51744 , ..., -75.71127 ,\n",
" -76.041084 , -76.38645 ],\n",
" ...,\n",
" [ 150.29921 , 151.56467 , 152.82253 , ..., 146.46725 ,\n",
" 147.74953 , 149.02718 ],\n",
" [ 137.3161 , 137.98648 , 138.65187 , ..., 135.27888 ,\n",
" 135.96191 , 136.64113 ],\n",
" [ 124.77108 , 124.96145 , 125.15015 , ..., 124.190674 ,\n",
" 124.38561 , 124.57911 ]]], dtype=float32)"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"''' Load tws processed data from shbundle as .mat file & from pyshbundle as xarray dataset '''\n",
"import scipy.io\n",
"# Load the .mat file\n",
"data = scipy.io.loadmat('../pyshbundle/data/validation_data/tws_sh.mat')\n",
"\n",
"# Access the variables in the .mat file\n",
"var1 = data['tws_m']\n",
"\n",
"temp=xr.open_dataset('../examples/tws_pysh.nc', engine=\"netcdf4\")\n",
"\n",
"var2=temp.tws.values"
]
},
{
"cell_type": "markdown",
"id": "25c2914c",
"metadata": {},
"source": [
"#### Lets convert both datasets to a netcdf format for easier calculations."
]
},
{
"cell_type": "markdown",
"id": "d18571e3",
"metadata": {},
"source": [
"* Converting `pyshbundle` processed data into netcdf format using xarray, to `ds_pysh`"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "8f7bbce8",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
<xarray.Dataset> Size: 51MB\n",
"Dimensions: (time: 196, lat: 180, lon: 360)\n",
"Coordinates:\n",
" * time (time) datetime64[ns] 2kB 2002-04-30 2002-05-31 ... 2022-06-30\n",
" * lat (lat) int64 1kB 89 88 87 86 85 84 83 ... -85 -86 -87 -88 -89 -90\n",
" * lon (lon) int64 3kB -180 -179 -178 -177 -176 ... 175 176 177 178 179\n",
"Data variables:\n",
" tws (time, lat, lon) float32 51MB -2.598 -2.549 -2.499 ... 124.4 124.6\n",
"Attributes:\n",
" units: mm<xarray.Dataset> Size: 102MB\n",
"Dimensions: (time: 196, lat: 180, lon: 360)\n",
"Coordinates:\n",
" * time (time) datetime64[ns] 2kB 2002-04-30 2002-05-31 ... 2022-06-30\n",
" * lat (lat) int64 1kB 89 88 87 86 85 84 83 ... -85 -86 -87 -88 -89 -90\n",
" * lon (lon) int64 3kB -180 -179 -178 -177 -176 ... 175 176 177 178 179\n",
"Data variables:\n",
" tws (time, lat, lon) float64 102MB -2.598 -2.549 -2.499 ... 124.4 124.6\n",
"Attributes:\n",
" units: mm<xarray.DataArray 'tws' (time: 196)> Size: 2kB\n",
"array([ 5.87514242e-09, -4.26223627e-09, -1.39690889e-08, -7.80649822e-09,\n",
" -6.42470429e-09, -8.19055354e-09, 6.63720314e-09, 9.31254960e-09,\n",
" -5.73921703e-09, 1.89799448e-09, 8.24561639e-09, 3.79555209e-09,\n",
" 1.93735599e-09, -1.83972199e-09, 3.81898356e-09, -2.62595984e-09,\n",
" -2.38502441e-09, -1.87216996e-10, -1.11288437e-09, 9.30366366e-10,\n",
" 2.59444324e-09, 2.16750338e-10, 2.78564623e-09, -7.79717517e-09,\n",
" 1.23423911e-09, -6.35305366e-09, 1.16946797e-10, 1.88048715e-09,\n",
" -4.28015162e-09, -2.50727195e-09, 9.53246503e-09, 6.83303813e-09,\n",
" 2.47475427e-09, 8.28576528e-10, -1.35008986e-09, -1.28247128e-09,\n",
" -5.63241678e-09, 2.24805231e-09, -3.87196555e-10, -6.36321565e-09,\n",
" -6.67633858e-09, 1.87588679e-09, 2.03895521e-09, 2.05828367e-09,\n",
" 1.07927121e-11, -1.64327905e-09, 2.35859809e-09, 3.71154031e-10,\n",
" -2.67425211e-09, -6.57371404e-09, 5.83918020e-10, -2.58793085e-09,\n",
" -1.10072471e-09, 3.78740363e-09, 1.29334811e-09, -3.38107645e-09,\n",
" 7.70597853e-11, 3.39068548e-11, -9.98683181e-09, -6.64013796e-09,\n",
" -7.95609364e-12, -1.23654871e-09, -1.55860070e-09, -2.62608386e-09,\n",
" 5.01679307e-09, -3.52707908e-09, -6.91500705e-09, -1.66992629e-09,\n",
" 1.49260014e-10, 4.41058990e-09, -2.50259149e-09, -1.05710635e-08,\n",
" -5.70400550e-09, 3.75981063e-09, 1.72891964e-09, 4.58359305e-09,\n",
" -6.86300811e-09, 3.69887100e-09, 9.94802864e-09, -6.34540876e-09,\n",
"...\n",
" -1.11768902e-09, -3.87988484e-09, 1.34124587e-08, -1.02895904e-08,\n",
" 1.17854943e-08, 4.70972124e-10, 2.21671436e-09, 2.73956553e-09,\n",
" -1.18020195e-09, 2.38018681e-10, 6.36812381e-09, -1.67224391e-09,\n",
" -6.25125249e-09, -2.01698382e-09, 8.31228049e-09, -5.87187740e-09,\n",
" 1.25964770e-09, -1.22692640e-09, 6.42588051e-09, -2.44977039e-09,\n",
" -3.71793112e-09, 3.94672049e-09, 2.38468517e-09, -7.19730792e-10,\n",
" 9.95087958e-09, 1.84939524e-09, 4.66824203e-09, 1.02633004e-08,\n",
" -6.80598554e-09, 1.77905552e-08, 8.67085266e-09, 1.34698485e-08,\n",
" 1.69740135e-10, -1.13187689e-08, 6.77686765e-09, -2.11015178e-10,\n",
" 8.87092181e-09, 6.98568608e-10, -1.24434164e-08, 1.12154009e-08,\n",
" 6.78988250e-09, -1.17979583e-09, -1.28954413e-09, 3.04643606e-10,\n",
" 4.69958849e-09, 4.11638867e-09, 1.30736892e-08, -1.37415186e-08,\n",
" 1.06418072e-08, 9.08527564e-09, 4.90118193e-09, 6.11483914e-09,\n",
" -6.40441328e-09, 1.84915147e-08, -6.10342800e-09, 6.43397894e-09,\n",
" 1.89346129e-08, 4.24687298e-09, 3.45359202e-09, 2.99639133e-10,\n",
" -7.83841821e-10, 4.93612486e-09, -2.34304746e-09, 4.80692647e-09,\n",
" -9.03242491e-09, -2.54978475e-09, 6.18848239e-09, 1.40658390e-08,\n",
" -4.17258348e-09, 6.13656311e-09, 4.88619985e-09, 1.56514421e-08,\n",
" 3.92759882e-09, -5.24303300e-09, 1.21710124e-08, 8.05693334e-09,\n",
" 9.01684360e-10, 8.92289255e-09, 1.87684341e-09, -1.58011433e-08])\n",
"Coordinates:\n",
" * time (time) datetime64[ns] 2kB 2002-04-30 2002-05-31 ... 2022-06-30\n",
" lat int64 8B 89\n",
" lon int64 8B -180<xarray.DataArray 'tws' (time: 243)> Size: 2kB\n",
"array([-9.69094529e-07, 2.36812042e-07, nan, nan,\n",
" 8.72978152e-08, -1.22263451e-07, 4.77928275e-08, 6.80909089e-08,\n",
" -1.37740276e-07, 1.01072999e-08, -7.46400161e-07, 9.89647361e-07,\n",
" 5.14377902e-07, -2.70129902e-07, nan, -6.01012488e-07,\n",
" -3.09933981e-08, 6.28074375e-07, -5.55916685e-07, 3.76842447e-07,\n",
" -5.70394079e-08, 2.17199208e-07, 3.39967215e-07, 8.08937116e-07,\n",
" 1.29339418e-06, 5.10734992e-07, -9.80370345e-08, nan,\n",
" nan, nan, nan, 1.06307496e-07,\n",
" 8.68036309e-09, -7.27863338e-07, 2.53009333e-07, -1.59850345e-07,\n",
" -2.06998760e-06, 8.39926741e-07, -2.76245459e-07, 9.98629304e-08,\n",
" 7.68413599e-07, 5.36044709e-07, -6.23643089e-07, -3.18224664e-07,\n",
" -5.31668150e-08, -1.18423031e-08, -7.10703212e-08, -5.07394873e-08,\n",
" 3.89335582e-07, 1.04053072e-07, 4.57618640e-07, -4.59484362e-09,\n",
" 2.88781294e-07, -1.25456124e-06, 7.81874121e-09, 1.07865930e-07,\n",
" -3.85457177e-07, 2.34572131e-08, -2.31180174e-07, -3.60701080e-07,\n",
" 2.07377070e-07, -1.33825765e-07, -5.05677363e-07, 4.24036536e-07,\n",
" 1.33390401e-07, 1.30523659e-07, 9.86716245e-07, -2.10345803e-08,\n",
" 7.35452730e-08, -8.95017749e-09, -9.91731675e-08, 1.29662244e-07,\n",
" -1.78431570e-07, -5.90120649e-07, 1.17968489e-06, -3.21184501e-08,\n",
" -1.75152451e-07, 6.74823468e-07, 6.38817767e-07, -3.48368417e-07,\n",
"...\n",
" -1.13959693e-07, 7.35506944e-08, nan, -5.70835113e-08,\n",
" nan, -1.35550025e-06, -2.58927741e-08, 1.28360909e-07,\n",
" -1.60221802e-07, nan, nan, -7.22614857e-10,\n",
" 1.00813597e-07, 8.52889684e-07, nan, 2.22721461e-07,\n",
" 9.06335572e-07, -2.43086902e-07, nan, nan,\n",
" nan, nan, nan, nan,\n",
" nan, nan, nan, nan,\n",
" nan, nan, 8.60203485e-07, 4.81101566e-07,\n",
" nan, nan, 4.96292998e-08, -3.21344373e-08,\n",
" -1.32220993e-07, -1.40416368e-07, nan, 3.48717322e-08,\n",
" -5.65086509e-07, 1.99339013e-07, 1.42135383e-06, -1.75112262e-07,\n",
" 1.64046334e-07, 4.76149893e-07, -6.60354260e-08, -1.97107710e-07,\n",
" 3.18719103e-07, -2.36057595e-07, -2.61367182e-08, 1.02312743e-07,\n",
" 1.44168624e-08, -2.41245061e-07, 1.09676563e-07, 1.71292669e-08,\n",
" -4.57230357e-07, 3.95559994e-07, -2.18310092e-06, -4.80689124e-07,\n",
" -8.21612119e-07, -6.15576880e-07, -1.20758642e-07, -3.08953387e-08,\n",
" 1.51928063e-07, 4.51007708e-07, 4.07144896e-09, 3.63877753e-07,\n",
" -2.96939973e-07, 6.93324552e-07, 8.75828334e-08, -6.22073145e-07,\n",
" 6.41594710e-07, -5.43634883e-07, 6.60334599e-07, -3.24580469e-08,\n",
" -2.22809993e-09, 1.79623873e-07, -1.27318494e-08])\n",
"Coordinates:\n",
" * time (time) datetime64[ns] 2kB 2002-04-30 2002-05-31 ... 2022-06-30