Jupyter Notebook for Validation of Results¶

Overview¶

This notebook compares the results obtained from the two implementations, shbundle and pyshbundle, for a specific task.
We will focus on evaluating the following:

Root Mean Squared Error (RMSE)
Normalized Root Mean Squared Error (NRMSE)
Difference in global area weighted water budget closure.
The difference in time series for a given basin.

Evaluation Metrics¶

1. Root Mean Squared Error (RMSE)¶

The RMSE is a measure of the average magnitude of the errors between predicted and observed values.
It is calculated as follows:

$$ RMSE = \sqrt{\frac{1}{n}\Sigma_{i=1}^{n}{\Big({y_i - \bar{y}}\Big)^2}}$$

where:

$n$ is the number of observations
${y_i}$ is the $i_{th}$ observation
$\bar{y}$ is the mean of observations

2. Normalized Root Mean Squared Error (NRMSE)¶

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:

$$ NRMSE = \frac{RMSE}{\max(y) - \min(y)}$$

3. Difference in global area weighted water budget closure.¶

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.

4. Difference in basin-average Time Series¶

Finally we create a basin-average TWSA time-series from both the implementations.

Data Preparation¶

Before proceeding with the analysis, ensure that the data from both shbundle and pyshbundle are loaded into the notebook.

In [1]:

import numpy as np
import pandas as pd
import xarray as xr
import matplotlib.pyplot as plt
import os

In [3]:

''' Load tws processed data from shbundle as .mat file & from pyshbundle as xarray dataset '''
import scipy.io
# Load the .mat file
data = scipy.io.loadmat('../pyshbundle/data/validation_data/tws_sh.mat')

# Access the variables in the .mat file
var1 = data['tws_m']

temp=xr.open_dataset('../examples/tws_pysh.nc', engine="netcdf4")

var2=temp.tws.values

Out[3]:

array([[[  -2.5984123,   -2.5491967,   -2.4989052, ...,   -2.7395167,
           -2.693581 ,   -2.6465435],
        [  -8.085175 ,   -7.9737053,   -7.859015 , ...,   -8.400428 ,
           -8.298518 ,   -8.1934395],
        [ -12.36135  ,  -12.230325 ,  -12.094314 , ...,  -12.725532 ,
          -12.608853 ,  -12.487493 ],
        ...,
        [ -63.64362  ,  -64.29353  ,  -64.944595 , ...,  -61.70406  ,
          -62.34856  ,  -62.995186 ],
        [ -59.735207 ,  -60.10004  ,  -60.464703 , ...,  -58.641277 ,
          -59.00567  ,  -59.37036  ],
        [ -56.28362  ,  -56.395493 ,  -56.50705  , ...,  -55.946453 ,
          -56.05906  ,  -56.17146  ]],

       [[  -1.6946957,   -1.6432129,   -1.5906656, ...,   -1.8426621,
           -1.7944294,   -1.7451043],
        [  -7.195656 ,   -7.077713 ,   -6.9565587, ...,   -7.530373 ,
           -7.42197  ,   -7.310403 ],
        [ -11.6058655,  -11.462069 ,  -11.313266 , ...,  -12.008124 ,
          -11.878803 ,  -11.744747 ],
        ...,
        [ -60.940014 ,  -61.585247 ,  -62.23357  , ...,  -59.026295 ,
          -59.660206 ,  -60.29822  ],
        [ -58.615997 ,  -58.973927 ,  -59.332268 , ...,  -57.546318 ,
          -57.902035 ,  -58.258648 ],
        [ -56.071903 ,  -56.180115 ,  -56.288097 , ...,  -55.7463   ,
          -55.854958 ,  -55.963505 ]],

       [[  28.409975 ,   28.435665 ,   28.462093 , ...,   28.337358 ,
           28.360819 ,   28.385025 ],
        [  24.354195 ,   24.405733 ,   24.459566 , ...,   24.213146 ,
           24.257923 ,   24.304932 ],
        [  20.799301 ,   20.849913 ,   20.904062 , ...,   20.667936 ,
           20.708385 ,   20.75215  ],
        ...,
        [ -53.207798 ,  -53.74962  ,  -54.29029  , ...,  -51.579227 ,
          -52.12223  ,  -52.66521  ],
        [ -48.65204  ,  -48.940662 ,  -49.228447 , ...,  -47.782604 ,
          -48.07286  ,  -48.36272  ],
        [ -45.042896 ,  -45.128735 ,  -45.214226 , ...,  -44.783646 ,
          -44.870323 ,  -44.956753 ]],

       ...,

       [[ -91.37819  ,  -91.51497  ,  -91.65443  , ...,  -90.984245 ,
          -91.1128   ,  -91.244125 ],
        [ -78.2787   ,  -78.60071  ,  -78.93178  , ...,  -77.366776 ,
          -77.66176  ,  -77.96573  ],
        [ -67.629616 ,  -68.0414   ,  -68.46865  , ...,  -66.48567  ,
          -66.85189  ,  -67.23316  ],
        ...,
        [ 169.38713  ,  170.59637  ,  171.79698  , ...,  165.71684  ,
          166.94644  ,  168.17017  ],
        [ 155.91563  ,  156.5555   ,  157.19012  , ...,  153.96843  ,
          154.62172  ,  155.27092  ],
        [ 143.16539  ,  143.34567  ,  143.52426  , ...,  142.61514  ,
          142.80006  ,  142.9835   ]],

       [[-125.82837  , -125.960236 , -126.09478  , ..., -125.44908  ,
         -125.57277  , -125.6992   ],
        [-112.4685   , -112.77872  , -113.098    , ..., -111.592    ,
         -111.87517  , -112.16732  ],
        [-101.242386 , -101.638016 , -102.049255 , ..., -100.147675 ,
         -100.49736  , -100.86221  ],
        ...,
        [ 149.24057  ,  150.4772   ,  151.70389  , ...,  145.48059  ,
          146.7413   ,  147.99495  ],
        [ 134.43509  ,  135.08853  ,  135.73586  , ...,  132.44183  ,
          133.11136  ,  133.7759   ],
        [ 120.18574  ,  120.36971  ,  120.551735 , ...,  119.62282  ,
          119.81222  ,  119.999886 ]],

       [[-101.70098  , -101.829636 , -101.96099  , ..., -101.33151  ,
         -101.4519   , -101.575066 ],
        [ -88.18868  ,  -88.48786  ,  -88.79633  , ...,  -87.346596 ,
          -87.61808  ,  -87.898766 ],
        [ -76.74753  ,  -77.12447  ,  -77.51744  , ...,  -75.71127  ,
          -76.041084 ,  -76.38645  ],
        ...,
        [ 150.29921  ,  151.56467  ,  152.82253  , ...,  146.46725  ,
          147.74953  ,  149.02718  ],
        [ 137.3161   ,  137.98648  ,  138.65187  , ...,  135.27888  ,
          135.96191  ,  136.64113  ],
        [ 124.77108  ,  124.96145  ,  125.15015  , ...,  124.190674 ,
          124.38561  ,  124.57911  ]]], dtype=float32)

Lets convert both datasets to a netcdf format for easier calculations.¶

Converting pyshbundle processed data into netcdf format using xarray, to ds_pysh

In [4]:

gs=1;
lon = np.arange(-180,180,gs)
lat = np.arange(89,-91,-gs)
ds_pysh = xr.Dataset(
    data_vars=dict(
        tws=(["time","lat", "lon"], var2)
    ),
    coords = {
        "time":(('time'),temp.time.data),
        "lat":lat,
        "lon":lon },
        attrs=dict(units="mm",),
)

Same with shubundle processed data, to ds_msh

In [5]:

ds_msh = xr.Dataset(
    data_vars=dict(
        tws=(["time","lat", "lon"], var1)
    ),
    coords = {
        "time":(('time'),temp.time.data),
        "lat":lat,
        "lon":lon },
    attrs=dict(units="mm",),
)

1. Gridwise RMSE calculation¶

Before finding the grid-wise RMSE values we need to ignore the data for the missing months.

In [8]:

# Calculate the difference between the two datasets
diff = ds_msh['tws'].dropna(dim='time').values - ds_pysh['tws'].dropna(dim='time').values    # dropna is used to remove nan values, the dates where the GRACE data is missing

# Calculate the squared difference
squared_diff = diff**2

# Calculate the mean squared difference along the time axis
mean_squared_diff = np.mean(squared_diff, axis=0)

# Calculate the root mean squared error (RMSE)
gridwise_rmse = np.sqrt(mean_squared_diff);

In [10]:

fig, (ax1) = plt.subplots(figsize=(6, 3), dpi=150, layout="constrained")    #gridspec_kw={'wspace': 0.01, 'hspace':0.2} 

# Set ticks and labels for the first subplot (ax1)
ax1.set_xticks(np.linspace(0, 359, 7)) 
ax1.set_xticklabels(['180W', '120W', '60W', '0', '60E', '120E', '180E'])
ax1.set_yticks(np.linspace(0, 179, 5))
ax1.set_yticklabels(['90N', '45N', '0', '45S', '90S'])

im1 = ax1.imshow(gridwise_rmse, cmap='jet')
cbar1 = fig.colorbar(im1, ax=ax1, pad=0.0, shrink=0.9)
ax1.set_title('RMSE between shbundle & pyshbundle generated TWS', fontsize=10, y=1.0)


plt.savefig('./figs/02_error_rmse.png', dpi=300, bbox_inches='tight')
plt.show()

In [11]:

if np.all(gridwise_rmse < 1e-3):
    print('The two datasets are almost identical')

The two datasets are almost identical

2. Gridwise NRMSE¶

In [12]:

# Calculate the normalized root mean squared error (NRMSE)
gridwise_nrmse = gridwise_rmse/np.std(ds_msh['tws'].dropna(dim='time').values, axis=0)

Out[12]:

array([[2.60043166e-08, 2.83013224e-08, 2.73803129e-08, ...,
        2.38515174e-08, 2.78223459e-08, 2.74585735e-08],
       [2.45173028e-08, 2.85903603e-08, 2.61927912e-08, ...,
        2.38695757e-08, 2.64437658e-08, 2.92895277e-08],
       [2.77640184e-08, 2.53198445e-08, 2.56853287e-08, ...,
        2.84399610e-08, 2.72859504e-08, 2.70947129e-08],
       ...,
       [2.97354718e-08, 3.33065197e-08, 3.47337847e-08, ...,
        3.32066023e-08, 3.68647979e-08, 3.25539337e-08],
       [3.43433379e-08, 3.48153238e-08, 3.33528178e-08, ...,
        3.62858782e-08, 3.30633574e-08, 3.46478197e-08],
       [3.42314674e-08, 2.61329397e-08, 3.27963439e-08, ...,
        3.12500385e-08, 3.01157924e-08, 3.60173407e-08]])

In [14]:

fig, (ax2) = plt.subplots(figsize=(6, 3), dpi=150, layout="constrained")    #gridspec_kw={'wspace': 0.01, 'hspace':0.2} 

# Set ticks and labels for the first subplot (ax2)
ax2.set_xticks(np.linspace(0, 359, 7))  # Adjust as needed
ax2.set_xticklabels(['180W', '120W', '60W', '0', '60E', '120E', '180E'])  # Adjust as needed
ax2.set_yticks(np.linspace(0, 179, 5))  # Adjust as needed
ax2.set_yticklabels(['90N', '45N', '0', '45S', '90S'])  # Adjust as needed

im1 = ax2.imshow(gridwise_nrmse, cmap='jet')
cbar1 = fig.colorbar(im1, ax=ax2, pad=0.0, shrink=0.9)
ax2.set_title('NRMSE between shbundle & pyshbundle generated TWS', fontsize=10, y=1)


plt.savefig('./figs/02_error_nrmse.png', dpi=300, bbox_inches='tight')
plt.show()

3. Global area weighted water budget closure¶

In [15]:

# Area of grids
from pyshbundle.hydro import area_weighting
global_grid_area=area_weighting(1)
global_grid_area_sum = np.sum(global_grid_area)
print('global surface area in m\u00b2:', global_grid_area_sum)

Initializing PySHbundle v0.3.0
global surface area in m²: 511207893395810.94

Calculate the global area weighted water budget closure error

In [16]:

# Create a copy of the datasets
ds_msh_area_weighted, ds_pysh_area_weighted = ds_msh.copy(), ds_pysh.copy()

# Area weight with the global grid area and calculate the sum over lat and lon
ds_msh_area_weighted = ds_msh['tws']*global_grid_area / global_grid_area_sum
ds_msh_area_weighted = ds_msh_area_weighted.sum(dim=['lat', 'lon'])

# Same for the pyshbundle dataset
ds_pysh_area_weighted = ds_pysh['tws']*global_grid_area / global_grid_area_sum
ds_pysh_area_weighted = ds_pysh_area_weighted.sum(dim=['lat', 'lon'])

diff_global = ds_msh_area_weighted - ds_pysh_area_weighted

# Reinsert the NaN values where the GRACE data is missing and the time coordinate
ds_msh_area_weighted=ds_msh_area_weighted.where(~np.isnan(ds_pysh['tws'][:,0,0]), np.nan)
ds_pysh_area_weighted=ds_pysh_area_weighted.where(~np.isnan(ds_pysh['tws'][:,0,0]), np.nan)
diff_global=diff_global.where(~np.isnan(ds_pysh['tws'][:,0,0]), np.nan)

In [17]:

# Create a figure with three subplots
fig, axes = plt.subplots(3, 1, figsize=(12, 15),)
fig.suptitle('Global water budget closure', fontsize=16)

# Plot 1
ds_msh_area_weighted.plot(ax=axes[0], color='#807dba', label='SHbundle')
axes[0].set_title('a) Water Budget Closure with SHbundle', fontsize=14)
axes[0].set_ylabel('Area Weighted TWSA (in cm)', fontsize=9)
axes[0].legend()
axes[0].grid(True, alpha=0.7)

# Plot 2
ds_pysh_area_weighted.plot(ax=axes[1], color='#4292c6', label='PySHbundle')
axes[1].set_title('b) Water Budget Closure with PySHbundle', fontsize=14)
axes[1].set_ylabel('Area Weighted TWSA (in cm)', fontsize=9)
axes[1].legend()
axes[1].grid(True, alpha=0.7)

# Plot 3
diff_global.plot(ax=axes[2], color='#525252', label='Difference')
axes[2].set_title('Difference between SHbundle and PySHbundle (a-b)', fontsize=14)
axes[2].set_ylabel('Area Weighted TWSA (in cm)', fontsize=9)
axes[2].legend()
axes[2].grid(True, alpha=0.7)

# Adjust layout for better spacing
plt.tight_layout()
plt.show()

4. Difference in basin-average Time Series¶

In [18]:

import geopandas as gpd

In [19]:

path_shapefile = '../pyshbundle/data/mrb_shapefiles/mrb_basins.shp'
shp = gpd.read_file(path_shapefile)
shp.plot(figsize=(8, 4))  

Out[19]:

<Axes: >

In [20]:

basin_name='KRISHNA'
shp_basin=shp[shp['RIVER_BASI']==basin_name]
print(shp_basin.head(), '\n')
shp_basin.plot()
basin_area=np.float64(shp_basin['SUM_SUB_AR'].values[0])*1e6          # basin area already in m^2
print('Basin area is :', basin_area, 'm\u00b2');

     MRBID RIVER_BASI CONTINENT            SEA         OCEAN  SUM_SUB_AR  \
115   2312    KRISHNA      Asia  Bay of Bengal  Indian Ocean    258785.1   

     Shape_Leng  Shape_Area                                           geometry  
115   36.078865   21.832582  POLYGON ((73.99167 19.37917, 73.99202 19.37674...   

Basin area is : 258785100000.0 m²

In [21]:

from pyshbundle.hydro import Basinaverage
_, basin_avg_tws_msh = Basinaverage(ds_msh, gs, shp_basin, basin_area)
_, basin_avg_tws_pysh = Basinaverage(ds_pysh, gs, shp_basin, basin_area)

In [22]:

new_dates=pd.date_range(start=basin_avg_tws_msh.time[0].values, 
                        end=basin_avg_tws_msh.time[-1].values, freq='ME',)
new_dates

Out[22]:

DatetimeIndex(['2002-04-30', '2002-05-31', '2002-06-30', '2002-07-31',
               '2002-08-31', '2002-09-30', '2002-10-31', '2002-11-30',
               '2002-12-31', '2003-01-31',
               ...
               '2021-09-30', '2021-10-31', '2021-11-30', '2021-12-31',
               '2022-01-31', '2022-02-28', '2022-03-31', '2022-04-30',
               '2022-05-31', '2022-06-30'],
              dtype='datetime64[ns]', length=243, freq='ME')

In [23]:

# Empty dataset for the gapped data, msh
basin_avg_tws_gapped_msh = xr.Dataset(
        data_vars = dict(   tws=(["time"], np.nan*np.arange(len(new_dates)))),
        coords=dict(time=new_dates),);
# Empty dataset for the gapped data, pysh
basin_avg_tws_gapped_pysh = xr.Dataset(
        data_vars = dict(   tws=(["time"], np.nan*np.arange(len(new_dates)))),
        coords=dict(time=new_dates),);

In [24]:

#
basin_avg_tws_gapped_msh['tws'] = basin_avg_tws_msh['tws'].where(
    basin_avg_tws_msh['time'].isin(basin_avg_tws_gapped_msh['time']),)
#
basin_avg_tws_gapped_pysh['tws'] = basin_avg_tws_pysh['tws'].where(
    basin_avg_tws_pysh['time'].isin(basin_avg_tws_gapped_pysh['time']),)

diff_global_basin = basin_avg_tws_gapped_msh['tws'] - basin_avg_tws_gapped_pysh['tws']

Out[24]:

<xarray.DataArray 'tws' (time: 243)> Size: 2kB
array([-9.69094529e-07,  2.36812042e-07,             nan,             nan,
        8.72978152e-08, -1.22263451e-07,  4.77928275e-08,  6.80909089e-08,
       -1.37740276e-07,  1.01072999e-08, -7.46400161e-07,  9.89647361e-07,
        5.14377902e-07, -2.70129902e-07,             nan, -6.01012488e-07,
       -3.09933981e-08,  6.28074375e-07, -5.55916685e-07,  3.76842447e-07,
       -5.70394079e-08,  2.17199208e-07,  3.39967215e-07,  8.08937116e-07,
        1.29339418e-06,  5.10734992e-07, -9.80370345e-08,             nan,
                   nan,             nan,             nan,  1.06307496e-07,
        8.68036309e-09, -7.27863338e-07,  2.53009333e-07, -1.59850345e-07,
       -2.06998760e-06,  8.39926741e-07, -2.76245459e-07,  9.98629304e-08,
        7.68413599e-07,  5.36044709e-07, -6.23643089e-07, -3.18224664e-07,
       -5.31668150e-08, -1.18423031e-08, -7.10703212e-08, -5.07394873e-08,
        3.89335582e-07,  1.04053072e-07,  4.57618640e-07, -4.59484362e-09,
        2.88781294e-07, -1.25456124e-06,  7.81874121e-09,  1.07865930e-07,
       -3.85457177e-07,  2.34572131e-08, -2.31180174e-07, -3.60701080e-07,
        2.07377070e-07, -1.33825765e-07, -5.05677363e-07,  4.24036536e-07,
        1.33390401e-07,  1.30523659e-07,  9.86716245e-07, -2.10345803e-08,
        7.35452730e-08, -8.95017749e-09, -9.91731675e-08,  1.29662244e-07,
       -1.78431570e-07, -5.90120649e-07,  1.17968489e-06, -3.21184501e-08,
       -1.75152451e-07,  6.74823468e-07,  6.38817767e-07, -3.48368417e-07,
...
       -1.13959693e-07,  7.35506944e-08,             nan, -5.70835113e-08,
                   nan, -1.35550025e-06, -2.58927741e-08,  1.28360909e-07,
       -1.60221802e-07,             nan,             nan, -7.22614857e-10,
        1.00813597e-07,  8.52889684e-07,             nan,  2.22721461e-07,
        9.06335572e-07, -2.43086902e-07,             nan,             nan,
                   nan,             nan,             nan,             nan,
                   nan,             nan,             nan,             nan,
                   nan,             nan,  8.60203485e-07,  4.81101566e-07,
                   nan,             nan,  4.96292998e-08, -3.21344373e-08,
       -1.32220993e-07, -1.40416368e-07,             nan,  3.48717322e-08,
       -5.65086509e-07,  1.99339013e-07,  1.42135383e-06, -1.75112262e-07,
        1.64046334e-07,  4.76149893e-07, -6.60354260e-08, -1.97107710e-07,
        3.18719103e-07, -2.36057595e-07, -2.61367182e-08,  1.02312743e-07,
        1.44168624e-08, -2.41245061e-07,  1.09676563e-07,  1.71292669e-08,
       -4.57230357e-07,  3.95559994e-07, -2.18310092e-06, -4.80689124e-07,
       -8.21612119e-07, -6.15576880e-07, -1.20758642e-07, -3.08953387e-08,
        1.51928063e-07,  4.51007708e-07,  4.07144896e-09,  3.63877753e-07,
       -2.96939973e-07,  6.93324552e-07,  8.75828334e-08, -6.22073145e-07,
        6.41594710e-07, -5.43634883e-07,  6.60334599e-07, -3.24580469e-08,
       -2.22809993e-09,  1.79623873e-07, -1.27318494e-08])
Coordinates:
  * time     (time) datetime64[ns] 2kB 2002-04-30 2002-05-31 ... 2022-06-30

xarray.DataArray

'tws'

time: 243

-9.691e-07 2.368e-07 nan nan ... -2.228e-09 1.796e-07 -1.273e-08

array([-9.69094529e-07,  2.36812042e-07,             nan,             nan,
        8.72978152e-08, -1.22263451e-07,  4.77928275e-08,  6.80909089e-08,
       -1.37740276e-07,  1.01072999e-08, -7.46400161e-07,  9.89647361e-07,
        5.14377902e-07, -2.70129902e-07,             nan, -6.01012488e-07,
       -3.09933981e-08,  6.28074375e-07, -5.55916685e-07,  3.76842447e-07,
       -5.70394079e-08,  2.17199208e-07,  3.39967215e-07,  8.08937116e-07,
        1.29339418e-06,  5.10734992e-07, -9.80370345e-08,             nan,
                   nan,             nan,             nan,  1.06307496e-07,
        8.68036309e-09, -7.27863338e-07,  2.53009333e-07, -1.59850345e-07,
       -2.06998760e-06,  8.39926741e-07, -2.76245459e-07,  9.98629304e-08,
        7.68413599e-07,  5.36044709e-07, -6.23643089e-07, -3.18224664e-07,
       -5.31668150e-08, -1.18423031e-08, -7.10703212e-08, -5.07394873e-08,
        3.89335582e-07,  1.04053072e-07,  4.57618640e-07, -4.59484362e-09,
        2.88781294e-07, -1.25456124e-06,  7.81874121e-09,  1.07865930e-07,
       -3.85457177e-07,  2.34572131e-08, -2.31180174e-07, -3.60701080e-07,
        2.07377070e-07, -1.33825765e-07, -5.05677363e-07,  4.24036536e-07,
        1.33390401e-07,  1.30523659e-07,  9.86716245e-07, -2.10345803e-08,
        7.35452730e-08, -8.95017749e-09, -9.91731675e-08,  1.29662244e-07,
       -1.78431570e-07, -5.90120649e-07,  1.17968489e-06, -3.21184501e-08,
       -1.75152451e-07,  6.74823468e-07,  6.38817767e-07, -3.48368417e-07,
...
       -1.13959693e-07,  7.35506944e-08,             nan, -5.70835113e-08,
                   nan, -1.35550025e-06, -2.58927741e-08,  1.28360909e-07,
       -1.60221802e-07,             nan,             nan, -7.22614857e-10,
        1.00813597e-07,  8.52889684e-07,             nan,  2.22721461e-07,
        9.06335572e-07, -2.43086902e-07,             nan,             nan,
                   nan,             nan,             nan,             nan,
                   nan,             nan,             nan,             nan,
                   nan,             nan,  8.60203485e-07,  4.81101566e-07,
                   nan,             nan,  4.96292998e-08, -3.21344373e-08,
       -1.32220993e-07, -1.40416368e-07,             nan,  3.48717322e-08,
       -5.65086509e-07,  1.99339013e-07,  1.42135383e-06, -1.75112262e-07,
        1.64046334e-07,  4.76149893e-07, -6.60354260e-08, -1.97107710e-07,
        3.18719103e-07, -2.36057595e-07, -2.61367182e-08,  1.02312743e-07,
        1.44168624e-08, -2.41245061e-07,  1.09676563e-07,  1.71292669e-08,
       -4.57230357e-07,  3.95559994e-07, -2.18310092e-06, -4.80689124e-07,
       -8.21612119e-07, -6.15576880e-07, -1.20758642e-07, -3.08953387e-08,
        1.51928063e-07,  4.51007708e-07,  4.07144896e-09,  3.63877753e-07,
       -2.96939973e-07,  6.93324552e-07,  8.75828334e-08, -6.22073145e-07,
        6.41594710e-07, -5.43634883e-07,  6.60334599e-07, -3.24580469e-08,
       -2.22809993e-09,  1.79623873e-07, -1.27318494e-08])

Coordinates: (1)

time

(time)

datetime64[ns]

2002-04-30 ... 2022-06-30

array(['2002-04-30T00:00:00.000000000', '2002-05-31T00:00:00.000000000',
       '2002-06-30T00:00:00.000000000', ..., '2022-04-30T00:00:00.000000000',
       '2022-05-31T00:00:00.000000000', '2022-06-30T00:00:00.000000000'],
      dtype='datetime64[ns]')

Indexes: (1)

time

PandasIndex

PandasIndex(DatetimeIndex(['2002-04-30', '2002-05-31', '2002-06-30', '2002-07-31',
               '2002-08-31', '2002-09-30', '2002-10-31', '2002-11-30',
               '2002-12-31', '2003-01-31',
               ...
               '2021-09-30', '2021-10-31', '2021-11-30', '2021-12-31',
               '2022-01-31', '2022-02-28', '2022-03-31', '2022-04-30',
               '2022-05-31', '2022-06-30'],
              dtype='datetime64[ns]', name='time', length=243, freq='ME'))

Attributes: (0)

In [25]:

# Create a figure with three subplots
fig, axes = plt.subplots(3, 1, figsize=(12, 15),)
fig.suptitle('Basin average comparision', fontsize=16)

# Plot 1
basin_avg_tws_gapped_msh['tws'].plot(ax=axes[0], color='#807dba', label='SHbundle')
axes[0].set_title('a) Basin average with SHbundle', fontsize=14)
axes[0].set_ylabel('TWSA (in cm)', fontsize=9)
axes[0].legend()
axes[0].grid(True, alpha=0.7)

# Plot 2
basin_avg_tws_gapped_pysh['tws'].plot(ax=axes[1], color='#4292c6', label='PySHbundle')
axes[1].set_title('b) Basin average with with PySHbundle', fontsize=14)
axes[1].set_ylabel('TWSA (in cm)', fontsize=9)
axes[1].legend()
axes[1].grid(True, alpha=0.7)

# Plot 3
diff_global_basin.plot(ax=axes[2], color='#525252', label='Difference')
axes[2].set_title('Difference between SHbundle and PySHbundle (a-b)', fontsize=14)
axes[2].set_ylabel('Area Weighted TWSA (in cm)', fontsize=9)
axes[2].legend()
axes[2].grid(True, alpha=0.7)

# Adjust layout for better spacing
plt.tight_layout()
plt.show()

Thank you for reading!