TopoPyScale Benchmarking Against Local Observation¶

S. Filhol

In [32]:

import requests, os, glob
import pandas as pd
import matplotlib.pyplot as plt
%matplotlib
import numpy as np
import xarray as xr
import seaborn as sns
plt.style.use('style_1')
import rasterio
from scipy import stats

Using matplotlib backend: Qt5Agg

In [20]:

df_stn = pd.read_csv('inputs/dem/station_list.csv')
ds_stn = xr.open_dataset('outputs/stn_downscaled.nc')

In [21]:

df_obs = pd.concat(map(pd.read_pickle, glob.glob(os.path.join('', "inputs/obs/metno*.pckl"))))

In [22]:

df_obs.elementId.unique()

Out[22]:

array(['air_temperature', 'sum(precipitation_amount PT12H)', 'wind_speed',
       'wind_from_direction', 'relative_humidity',
       'air_pressure_at_sea_level'], dtype=object)

In [38]:

dem = rasterio.open('inputs/dem/ASTER_Finse_Kris.tif')
plt.figure()
show(dem)
plt.scatter(df_stn.x, df_stn.y, s=100)
for i,row in df_stn.iterrows():
    plt.text(row.x,row.y,row.Name)

In [24]:

# Convert pandas dataframe to xarray dataset
df = pd.pivot_table(df_obs, columns=['elementId'],values=['value'], index=['sourceId', 'referenceTime'])
ds = df.xs('value', axis=1, level=0).to_xarray()      
ds['referenceTime']=pd.to_datetime(ds.referenceTime)

In [61]:

stn_name='SN25830:0'

stn = pd.DataFrame()
stn['time'] = ds.referenceTime
stn['Downscaled'] = (ds_stn.t.sel(point_id=df_stn.loc[df_stn.stn_number==stn_name[:-2]].index)-273.15).copy().drop('point_id').sel(time=slice('2018-01-01', '2021-08-30'))[0].values
stn['Observed'] = ds.air_temperature.sel(sourceId=stn_name).dropna(dim='referenceTime').sel(referenceTime=slice('2018-01-01', '2021-08-30')).copy().values

Compare Air Temperature¶

In [49]:

# plot temperature for each station
for stn_name in df_stn.stn_number:
    stn = ds.air_temperature.sel(sourceId=stn_name+':0').dropna(dim='referenceTime').copy().to_dataframe().resample('1H').mean()
    stn['Downscaled'] = (ds_stn.t.sel(point_id=df_stn.loc[df_stn.stn_number==stn_name].index)-273.15).copy().drop('point_id').sel(time=stn.index)[0].values
    stn.rename(columns={'air_temperature':'Observed', 'sourceId':'stn_id'}, inplace=True)
    #stn = stn.dropna(axis=0)
    stn['var_diff'] = stn.Downscaled - stn.Observed

    plt.figure()
    ax = sns.histplot(x=stn.Observed, y=stn.Downscaled, cmap='magma')
    ax.figure.colorbar(ax.collections[0])
    plt.plot([-20,20], [-20,20], c='r', label='1:1 line')
    reg = stats.linregress(stn.Observed, stn.Downscaled)
    plt.plot([-20,20], np.array([-20,20])*reg.slope + reg.intercept, c='b', label='regression')
    plt.ylabel('Downscaled')
    plt.xlabel('Observed')
    plt.title('Air Temperature [$^{o}C$] at ' + df_stn.Name.loc[df_stn.stn_number==stn_name].values[0])
    plt.legend()
    
    # Plot time series 
    fig, ax = plt.subplots(2,1, sharex=True, gridspec_kw={'height_ratios': [3, 1]})
    # absolute value
    stn.Observed.plot(ax=ax[0], label='Obs.')
    stn.Downscaled.plot(ax=ax[0], label='TopoPyScale')
    ax[0].legend()
    ax[0].set_ylabel('Tair [$^{o}C$]')
    ax[0].set_title(df_stn.Name.loc[df_stn.stn_number==stn_name].values[0])
    # timeseries difference
    stn.var_diff.plot(ax=ax[1])
    (stn.Observed*0).plot(color='k', linestyle='--')
    ax[1].fill_between(stn.index, (stn.Observed*0)+stn.var_diff.quantile(0.05),(stn.Observed*0)+stn.var_diff.quantile(0.95), alpha=.3, color='r', zorder=5, label='90% values')
    ((stn.Observed*0)+stn.var_diff.median()).plot(ax=ax[1], label='median', c='r', linestyle='-', zorder=10)
    ax[1].legend()
    ax[1].set_ylabel('Tair diff [$^{o}C$]')

/tmp/ipykernel_13568/1664442083.py:8: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).
  plt.figure()

Compare Wind Speed¶

In [96]:

# plot temperature for each station
for stn_name in df_stn.stn_number:
    try:
        stn = ds.wind_speed.sel(sourceId=stn_name+':0').dropna(dim='referenceTime').copy().to_dataframe().resample('1H').mean()
        stn['Downscaled'] = (ds_stn.ws.sel(point_id=df_stn.loc[df_stn.stn_number==stn_name].index)).copy().drop('point_id').sel(time=stn.index)[0].values
        stn.rename(columns={'wind_speed':'Observed', 'sourceId':'stn_id'}, inplace=True)
        #stn = stn.dropna(axis=0)
        stn['var_diff'] = stn.Downscaled - stn.Observed

        plt.figure()
        ax = sns.histplot(x=stn.Observed, y=stn.Downscaled, cmap='magma')
        ax.figure.colorbar(ax.collections[0])
        plt.plot([0,20], [0,20], c='r', label='1:1 line')
        reg = stats.linregress(stn.Observed, stn.Downscaled)
        plt.plot([0,20], np.array([0,20])*reg.slope + reg.intercept, c='b', label='regression')
        plt.ylabel('Downscaled')
        plt.xlabel('Observed')
        plt.title('Wind Speed [m.s$^{-1}$] at ' + df_stn.Name.loc[df_stn.stn_number==stn_name].values[0])
        plt.legend()

        # Plot time series 
        fig, ax = plt.subplots(2,1, sharex=True, gridspec_kw={'height_ratios': [3, 1]})
        # absolute value
        stn.Observed.plot(ax=ax[0], label='Obs.')
        stn.Downscaled.plot(ax=ax[0], label='TopoPyScale')
        ax[0].legend()
        ax[0].set_ylabel('Wind Speed [m.s$^{-1}$]')
        ax[0].set_title(df_stn.Name.loc[df_stn.stn_number==stn_name].values[0])
        # timeseries difference
        stn.var_diff.plot(ax=ax[1])
        (stn.Observed*0).plot(color='k', linestyle='--')
        ax[1].fill_between(stn.index, (stn.Observed*0)+stn.var_diff.quantile(0.05),(stn.Observed*0)+stn.var_diff.quantile(0.95), alpha=.3, color='r', zorder=5, label='90% values')
        ((stn.Observed*0)+stn.var_diff.median()).plot(ax=ax[1], label='median', c='r', linestyle='-', zorder=10)
        ax[1].legend()
        ax[1].set_ylabel('Ws diff [m.s$^{-1}$]')
    except Exception as e: 
        print(e)
        

/tmp/ipykernel_13568/2219099889.py:10: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).
  plt.figure()

list index out of range

In [97]:

list(ds)

Out[97]:

['air_pressure_at_sea_level',
 'air_temperature',
 'relative_humidity',
 'sum(precipitation_amount PT12H)',
 'wind_from_direction',
 'wind_speed']

In [100]:

ds['sum(precipitation_amount PT12H)'].sel(sourceId=stn_name+':0').dropna(dim='referenceTime')

Out[100]:

<xarray.DataArray 'sum(precipitation_amount PT12H)' (referenceTime: 2666)>
array([7.3, 4.8, 2.8, ..., 0. , 0. , 0. ])
Coordinates:
    sourceId       <U9 'SN53480:0'
  * referenceTime  (referenceTime) datetime64[ns] 2018-01-01T06:00:00 ... 202...

In [144]:

# plot weekly precip for each station
from matplotlib.colors import LogNorm
for stn_name in df_stn.stn_number:
    try:
        stn = ds['sum(precipitation_amount PT12H)'].sel(sourceId=stn_name+':0').dropna(dim='referenceTime').copy().to_dataframe().resample('7D').sum()
        stn['Downscaled'] = (ds_stn.tp.sel(point_id=df_stn.loc[df_stn.stn_number==stn_name].index)).copy().drop('point_id').to_dataframe().droplevel(0).resample('7D').sum()
        stn.rename(columns={'sum(precipitation_amount PT12H)':'Observed', 'sourceId':'stn_id'}, inplace=True)
        #stn = stn.dropna(axis=0)
        stn['var_diff'] = stn.Downscaled - stn.Observed

        plt.figure()
        
        ax = sns.lmplot(x='Observed', y='Downscaled', data = stn)
        plt.plot([0,140], [0,140], c='r', label='1:1 line')
        reg = stats.linregress(stn.Observed, stn.Downscaled)
        plt.plot([0,140], np.array([0,140])*reg.slope + reg.intercept, c='b', label='regression')
        plt.ylabel('Downscaled')
        plt.xlabel('Observed')
        plt.title('7 Days Precip at ' + df_stn.Name.loc[df_stn.stn_number==stn_name].values[0])
        plt.legend()

        # Plot time series 
        fig, ax = plt.subplots(2,1, sharex=True, gridspec_kw={'height_ratios': [3, 1]})
        # absolute value
        stn.Observed.plot(ax=ax[0], label='Obs.')
        stn.Downscaled.plot(ax=ax[0], label='TopoPyScale')
        ax[0].legend()
        ax[0].set_ylabel('precip [mm]')
        ax[0].set_title(df_stn.Name.loc[df_stn.stn_number==stn_name].values[0])
        # timeseries difference
        stn.var_diff.plot(ax=ax[1])
        (stn.Observed*0).plot(color='k', linestyle='--')
        ax[1].fill_between(stn.index, (stn.Observed*0)+stn.var_diff.quantile(0.05),(stn.Observed*0)+stn.var_diff.quantile(0.95), alpha=.3, color='r', zorder=5, label='90% values')
        ((stn.Observed*0)+stn.var_diff.median()).plot(ax=ax[1], label='median', c='r', linestyle='-', zorder=10)
        ax[1].legend()
        ax[1].set_ylabel('Precip diff [mm]')
    except Exception as e: 
        print(e)
        

Wind Direction¶

In [147]:

from windrose import WindroseAxes

In [201]:

for stn_name in df_stn.stn_number:
    try:
        stn = ds.wind_speed.sel(sourceId=stn_name+':0').dropna(dim='referenceTime').copy().to_dataframe().resample('1H').mean()
        stn['Downscaled_ws'] = (ds_stn.ws.sel(point_id=df_stn.loc[df_stn.stn_number==stn_name].index)).copy().drop('point_id').sel(time=stn.index)[0].values
        stn.rename(columns={'wind_speed':'Observed_ws', 'sourceId':'stn_id'}, inplace=True)
        #stn = stn.dropna(axis=0)
        ds['cos_dir'] = np.cos(np.deg2rad(ds.wind_from_direction))
        ds['sin_dir'] = np.sin(np.deg2rad(ds.wind_from_direction))

        stn['cos_dir'] = ds.cos_dir.sel(sourceId=stn_name+':0').dropna(dim='referenceTime').copy().to_dataframe().resample('1H').mean()
        stn['sin_dir'] = ds.sin_dir.sel(sourceId=stn_name+':0').dropna(dim='referenceTime').copy().to_dataframe().resample('1H').mean()
        stn['Observed_wd'] = (np.rad2deg(np.arctan2(stn.sin_dir, stn.cos_dir ))+360)%360
        stn['Downscaled_wd'] = np.rad2deg(ds_stn.wd.sel(point_id=df_stn.loc[df_stn.stn_number==stn_name].index)).copy().drop('point_id').sel(time=stn.index)[0].values
        
        fig, ax = plt.subplots(2,2, gridspec_kw={'height_ratios': [3, 1]})
        ax[0,0].hist(stn.Downscaled_wd, bins=100, density=True, label='Downscaled')
        ax[0,0].hist(stn.Observed_wd, bins=100, alpha=.5, density=True, label='Observed')
        ax[0,0].legend()
        ax[0,0].set_ylabel('Density')
        #ax[0].set_xlabel('Wind direction [$^{o}$]')
        ax[0,0].set_title(df_stn.Name.loc[df_stn.stn_number==stn_name].values[0])
        
        ax[1,0].hist(stn.Downscaled_wd.loc[stn.Downscaled_ws>=stn.Downscaled_ws.quantile(0.5)], bins=100, density=True, label='Downscaled')
        ax[1,0].hist(stn.Observed_wd.loc[stn.Observed_ws>=stn.Observed_ws.quantile(0.5)], bins=100, alpha=.5, density=True, label='Observed')
        ax[1,0].legend()
        ax[1,0].set_ylabel('Density')
        ax[1,0].set_xlabel('Wind direction [$^{o}$]')
        #ax[1,0].set_title(df_stn.Name.loc[df_stn.stn_number==stn_name].values[0])
        
        ax[0,1].hist(stn.Downscaled_ws, bins=100, density=True, label='Downscaled')
        ax[0,1].hist(stn.Observed_ws, bins=100, alpha=.5, density=True, label='Observed')
        ax[0,1].legend()
        #ax[0,1].set_ylabel('Density')
        #ax[0,1].set_xlabel('Wind speed [$m.s^{-1}$]')
        ax[0,1].set_xlim((0,20))
        
        ax[1,1].hist(stn.Downscaled_ws, bins=100, density=True, label='Downscaled', cumulative=True)
        ax[1,1].hist(stn.Observed_ws, bins=100, alpha=.5, density=True, label='Observed', cumulative=True)
        ax[1,1].legend()
        #ax[1,1].set_ylabel('Density')
        ax[1,1].set_xlabel('Wind speed [$m.s^{-1}$]')
        ax[1,1].plot([0,20], [0.5,0.5], c='r', linestyle='--')
        ax[1,1].set_xlim((0,20))
        
    except Exception as e: 
        print(e)      

Wrong number of items passed 2, placement implies 1

In [180]:

ax = WindroseAxes.from_ax()
ax.contourf(np.rad2deg(stn.Downscaled_wd), stn.Downscaled_ws, bins=np.arange(0, 10, 1), normed=True, cmap=plt.cm.hot)
ax.set_legend()

Out[180]:

<matplotlib.legend.Legend at 0x7f367beb4340>

In [194]:

plt.figure()
plt.hist(stn.Downscaled_ws, bins=100, density=True, label='Downscaled', cumulative=True)
plt.hist(stn.Observed_ws, bins=100, alpha=.5, density=True, label='Observed', cumulative=True)
plt.legend()
plt.ylabel('Density')
plt.xlabel('Wind speed [$m.s^{-1}$]')

Out[194]:

Text(0.5, 0, 'Wind speed [$m.s^{-1}$]')

In [193]:

stn.Observed_ws.quantile(0.7)

Out[193]:

5.283333333333333

In [39]:

# plot wind direction for each station
# TODO
for stn_name in df_stn.stn_number:
    try:
        stn = ds.wind_speed.sel(sourceId=stn_name+':0').dropna(dim='referenceTime').copy().to_dataframe().resample('1H').mean()
        stn['Downscaled'] = (ds_stn.ws.sel(point_id=df_stn.loc[df_stn.stn_number==stn_name].index)).copy().drop('point_id').sel(time=stn.index)[0].values
        stn.rename(columns={'wind_speed':'Observed', 'sourceId':'stn_id'}, inplace=True)
        #stn = stn.dropna(axis=0)
        stn['var_diff'] = stn.Downscaled - stn.Observed

        plt.figure()
        ax = sns.histplot(x=stn.Observed, y=stn.Downscaled, cmap='magma')
        ax.figure.colorbar(ax.collections[0])
        plt.plot([0,20], [0,20], c='r', label='1:1 line')
        reg = stats.linregress(stn.Observed, stn.Downscaled)
        plt.plot([0,20], np.array([0,20])*reg.slope + reg.intercept, c='b', label='regression')
        plt.ylabel('Downscaled')
        plt.xlabel('Observed')
        plt.title('Wind Speed [m.s$^{-1}$] at ' + df_stn.Name.loc[df_stn.stn_number==stn_name].values[0])
        plt.legend()

        # Plot time series 
        fig, ax = plt.subplots(2,1, sharex=True, gridspec_kw={'height_ratios': [3, 1]})
        # absolute value
        stn.Observed.plot(ax=ax[0], label='Obs.')
        stn.Downscaled.plot(ax=ax[0], label='TopoPyScale')
        ax[0].legend()
        ax[0].set_ylabel('Wind Speed [m.s$^{-1}$]')
        ax[0].set_title(df_stn.Name.loc[df_stn.stn_number==stn_name].values[0])
        # timeseries difference
        stn.var_diff.plot(ax=ax[1])
        (stn.Observed*0).plot(color='k', linestyle='--')
        ax[1].fill_between(stn.index, (stn.Observed*0)+stn.var_diff.quantile(0.05),(stn.Observed*0)+stn.var_diff.quantile(0.95), alpha=.3, color='r', zorder=5, label='90% values')
        ((stn.Observed*0)+stn.var_diff.median()).plot(ax=ax[1], label='median', c='r', linestyle='-', zorder=10)
        ax[1].legend()
        ax[1].set_ylabel('Ws diff [m.s$^{-1}$]')
    except Exception as e: 
        print(e)
        

In [ ]: