Calculating quantiles of timeseries data and then create a "fanplot"

Question

My code is working properly, but I am looking for better approach in the calculation of quantiles and the finding of the data in the dataframe.

import datetime
import pandas
import pandas as pd
import numpy as np 
from datetime import datetime, timedelta
import seaborn as sns
import matplotlib.pyplot as plt
colnames = ['Date', 'Energy']
df_1 = pd.read_csv('WTPV_Spring_2016.csv', names = colnames , encoding="utf8", delimiter=";")
start_date = datetime(2022, 9, 27, 1, 0)
end_date = datetime(2022, 9, 28, 1, 0)
def daterange(start_date, end_date):
 delta = timedelta(hours = 1)
 while start_date < end_date:
 yield start_date
 start_date += delta
 
df = pd.DataFrame(columns = 
['Time','pct0.1','pct0.2','pct0.3','pct0.4','pct0.5','pct0.6','pct0.7','pct0.8','pct0.9'])
for single_date in daterange(start_date, end_date):
 df.loc[single_date, ['Time']] = single_date.strftime("%H:%M:%S")
x = []
for index in df.index:
 y = df_1.loc[df_1['Date'].str.contains(df['Time'][index])]
 for i in np.arange(1, 10, 1)/10:
 x.append(y.quantile(i))
pct = pd.DataFrame(x, columns = ['Energy'])
df['pct0.1'] = pct.loc[0.1].values
df['pct0.2'] = pct.loc[0.2].values
df['pct0.3'] = pct.loc[0.3].values
df['pct0.4'] = pct.loc[0.4].values
df['pct0.5'] = pct.loc[0.5].values
df['pct0.6'] = pct.loc[0.6].values
df['pct0.7'] = pct.loc[0.7].values
df['pct0.8'] = pct.loc[0.8].values
df['pct0.9'] = pct.loc[0.9].values
sns.set(font_scale = 1.5, style = "white")
fig, ax = plt.subplots(figsize = (8, 6))
xs = np.arange(len(df))
colors = plt.cm.Greens(np.linspace(0.1, 0.6, 5))
for lower, upper, color in zip([f'pct0.{i}' for i in range(1, 5)], [f'pct0.{i}' for i in range(9, 5, -1)], colors):
 ax.fill_between(xs, df[lower], df[upper], color = color, label = lower + '-' + upper)
ax.plot(xs, df['pct0.5'], color = 'black' , lw = 1.5, label ='Median')
ax.set_xticks(xs)
ax.set_xticklabels(df['Time'], fontsize = 11)
ax.margins(x = 0)
ax.set_ylim(ymin = 0)
for sp in ['top', 'right']:
 ax.spines[sp].set_visible(True)
plt.xticks(rotation = 90)
plt.yticks(fontsize = 11)
plt.title('WT & PV Spring 2016', fontsize = 15)
plt.ylabel('Energy MWh', fontsize = 12)
plt.xlabel('Time', fontsize = 12)
plt.ylim([0, 3200])
plt.show()

My data have a format like this:

2016年09月01日 00:00:00;0
2016年09月01日 01:00:00;0
2016年09月01日 02:00:00;0
2016年09月01日 03:00:00;0
2016年09月01日 04:00:00;0
2016年09月01日 05:00:00;0
2016年09月01日 06:00:00;0
2016年09月01日 07:00:00;0
2016年09月01日 08:00:00;103
2016年09月01日 09:00:00;420
2016年09月01日 10:00:00;837
2016年09月01日 11:00:00;1213
2016年09月01日 12:00:00;1470
2016年09月01日 13:00:00;1594
2016年09月01日 14:00:00;1609
2016年09月01日 15:00:00;1506
2016年09月01日 16:00:00;1288
2016年09月01日 17:00:00;994
2016年09月01日 18:00:00;645
2016年09月01日 19:00:00;309
2016年09月01日 20:00:00;72
2016年09月01日 21:00:00;0
2016年09月01日 22:00:00;0
2016年09月01日 23:00:00;0
2016年09月02日 00:00:00;0
2016年09月02日 01:00:00;0
2016年09月02日 02:00:00;0
2016年09月02日 03:00:00;0
2016年09月02日 04:00:00;0
2016年09月02日 05:00:00;0
2016年09月02日 06:00:00;0
2016年09月02日 07:00:00;0
2016年09月02日 08:00:00;94
2016年09月02日 09:00:00;376
2016年09月02日 10:00:00;745
2016年09月02日 11:00:00;1085
2016年09月02日 12:00:00;1341
2016年09月02日 13:00:00;1496
2016年09月02日 14:00:00;1527
2016年09月02日 15:00:00;1457
2016年09月02日 16:00:00;1286
2016年09月02日 17:00:00;1005
2016年09月02日 18:00:00;659
2016年09月02日 19:00:00;326
2016年09月02日 20:00:00;76
2016年09月02日 21:00:00;0
2016年09月02日 22:00:00;0
2016年09月02日 23:00:00;0
2016年09月03日 00:00:00;0
2016年09月03日 01:00:00;0
2016年09月03日 02:00:00;0
2016年09月03日 03:00:00;0
2016年09月03日 04:00:00;0
2016年09月03日 05:00:00;0
2016年09月03日 06:00:00;0
2016年09月03日 07:00:00;0
2016年09月03日 08:00:00;100
2016年09月03日 09:00:00;418
2016年09月03日 10:00:00;840
2016年09月03日 11:00:00;1220
2016年09月03日 12:00:00;1483
2016年09月03日 13:00:00;1622
2016年09月03日 14:00:00;1648
2016年09月03日 15:00:00;1568
2016年09月03日 16:00:00;1374
2016年09月03日 17:00:00;1086
2016年09月03日 18:00:00;726
2016年09月03日 19:00:00;355
2016年09月03日 20:00:00;78
2016年09月03日 21:00:00;0
2016年09月03日 22:00:00;0
2016年09月03日 23:00:00;0
2016年09月04日 00:00:00;0
2016年09月04日 01:00:00;0
2016年09月04日 02:00:00;0
2016年09月04日 03:00:00;0
2016年09月04日 04:00:00;0
2016年09月04日 05:00:00;0
2016年09月04日 06:00:00;0
2016年09月04日 07:00:00;0
2016年09月04日 08:00:00;100
2016年09月04日 09:00:00;430
2016年09月04日 10:00:00;861
2016年09月04日 11:00:00;1253
2016年09月04日 12:00:00;1526
2016年09月04日 13:00:00;1664
2016年09月04日 14:00:00;1688
2016年09月04日 15:00:00;1610
2016年09月04日 16:00:00;1425
2016年09月04日 17:00:00;1138
2016年09月04日 18:00:00;762
2016年09月04日 19:00:00;371
2016年09月04日 20:00:00;79
2016年09月04日 21:00:00;0
2016年09月04日 22:00:00;0
2016年09月04日 23:00:00;0
2016年09月05日 00:00:00;0
2016年09月05日 01:00:00;0
2016年09月05日 02:00:00;0
2016年09月05日 03:00:00;0
2016年09月05日 04:00:00;0
2016年09月05日 05:00:00;0
2016年09月05日 06:00:00;0
2016年09月05日 07:00:00;0
2016年09月05日 08:00:00;98
2016年09月05日 09:00:00;430
2016年09月05日 10:00:00;865
2016年09月05日 11:00:00;1263
2016年09月05日 12:00:00;1539
2016年09月05日 13:00:00;1673
2016年09月05日 14:00:00;1697
2016年09月05日 15:00:00;1603
2016年09月05日 16:00:00;1384
2016年09月05日 17:00:00;1078
2016年09月05日 18:00:00;696
2016年09月05日 19:00:00;320
2016年09月05日 20:00:00;63
2016年09月05日 21:00:00;0
2016年09月05日 22:00:00;0
2016年09月05日 23:00:00;0
2016年09月06日 00:00:00;0
2016年09月06日 01:00:00;0
2016年09月06日 02:00:00;0
2016年09月06日 03:00:00;0
2016年09月06日 04:00:00;0
2016年09月06日 05:00:00;0
2016年09月06日 06:00:00;0
2016年09月06日 07:00:00;0
2016年09月06日 08:00:00;64
2016年09月06日 09:00:00;252
2016年09月06日 10:00:00;491
2016年09月06日 11:00:00;718
2016年09月06日 12:00:00;855
2016年09月06日 13:00:00;920
2016年09月06日 14:00:00;937
2016年09月06日 15:00:00;831
2016年09月06日 16:00:00;650
2016年09月06日 17:00:00;478
2016年09月06日 18:00:00;288
2016年09月06日 19:00:00;128
2016年09月06日 20:00:00;26
2016年09月06日 21:00:00;0
2016年09月06日 22:00:00;0
2016年09月06日 23:00:00;0
2016年09月07日 00:00:00;0
2016年09月07日 01:00:00;0
2016年09月07日 02:00:00;0
2016年09月07日 03:00:00;0
2016年09月07日 04:00:00;0
2016年09月07日 05:00:00;0
2016年09月07日 06:00:00;0
2016年09月07日 07:00:00;0
2016年09月07日 08:00:00;43
2016年09月07日 09:00:00;172
2016年09月07日 10:00:00;353
2016年09月07日 11:00:00;537
2016年09月07日 12:00:00;708
2016年09月07日 13:00:00;913
2016年09月07日 14:00:00;1011
2016年09月07日 15:00:00;1002
2016年09月07日 16:00:00;969
2016年09月07日 17:00:00;797
2016年09月07日 18:00:00;508
2016年09月07日 19:00:00;238
2016年09月07日 20:00:00;47
2016年09月07日 21:00:00;0
2016年09月07日 22:00:00;0
2016年09月07日 23:00:00;0
2016年09月08日 00:00:00;0
2016年09月08日 01:00:00;0
2016年09月08日 02:00:00;0
2016年09月08日 03:00:00;0
2016年09月08日 04:00:00;0
2016年09月08日 05:00:00;0
2016年09月08日 06:00:00;0
2016年09月08日 07:00:00;0
2016年09月08日 08:00:00;60
2016年09月08日 09:00:00;247
2016年09月08日 10:00:00;477
2016年09月08日 11:00:00;692
2016年09月08日 12:00:00;870
2016年09月08日 13:00:00;1040
2016年09月08日 14:00:00;1105
2016年09月08日 15:00:00;1033
2016年09月08日 16:00:00;905
2016年09月08日 17:00:00;701
2016年09月08日 18:00:00;444
2016年09月08日 19:00:00;215
2016年09月08日 20:00:00;42
2016年09月08日 21:00:00;0
2016年09月08日 22:00:00;0
2016年09月08日 23:00:00;0
2016年09月09日 00:00:00;0
2016年09月09日 01:00:00;0
2016年09月09日 02:00:00;0
2016年09月09日 03:00:00;0
2016年09月09日 04:00:00;0
2016年09月09日 05:00:00;0
2016年09月09日 06:00:00;0
2016年09月09日 07:00:00;0
2016年09月09日 08:00:00;60
2016年09月09日 09:00:00;276
2016年09月09日 10:00:00;555
2016年09月09日 11:00:00;809
2016年09月09日 12:00:00;995
2016年09月09日 13:00:00;1124
2016年09月09日 14:00:00;1160
2016年09月09日 15:00:00;1103
2016年09月09日 16:00:00;979
2016年09月09日 17:00:00;769
2016年09月09日 18:00:00;487
2016年09月09日 19:00:00;220
2016年09月09日 20:00:00;39
2016年09月09日 21:00:00;0
2016年09月09日 22:00:00;0
2016年09月09日 23:00:00;0
2016年09月10日 00:00:00;0
2016年09月10日 01:00:00;0
2016年09月10日 02:00:00;0
2016年09月10日 03:00:00;0
2016年09月10日 04:00:00;0
2016年09月10日 05:00:00;0
2016年09月10日 06:00:00;0
2016年09月10日 07:00:00;0
2016年09月10日 08:00:00;65
2016年09月10日 09:00:00;310
2016年09月10日 10:00:00;629
2016年09月10日 11:00:00;925
2016年09月10日 12:00:00;1152
2016年09月10日 13:00:00;1310
2016年09月10日 14:00:00;1349
2016年09月10日 15:00:00;1265
2016年09月10日 16:00:00;1096
2016年09月10日 17:00:00;845
2016年09月10日 18:00:00;531
2016年09月10日 19:00:00;240
2016年09月10日 20:00:00;41
2016年09月10日 21:00:00;0
2016年09月10日 22:00:00;0
2016年09月10日 23:00:00;0
2016年09月11日 00:00:00;0
2016年09月11日 01:00:00;0
2016年09月11日 02:00:00;0
2016年09月11日 03:00:00;0
2016年09月11日 04:00:00;0
2016年09月11日 05:00:00;0
2016年09月11日 06:00:00;0
2016年09月11日 07:00:00;0
2016年09月11日 08:00:00;73
2016年09月11日 09:00:00;359
2016年09月11日 10:00:00;726
2016年09月11日 11:00:00;1066
2016年09月11日 12:00:00;1305
2016年09月11日 13:00:00;1443
2016年09月11日 14:00:00;1472
2016年09月11日 15:00:00;1386
2016年09月11日 16:00:00;1197
2016年09月11日 17:00:00;914
2016年09月11日 18:00:00;569
2016年09月11日 19:00:00;254
2016年09月11日 20:00:00;41
2016年09月11日 21:00:00;0
2016年09月11日 22:00:00;0
2016年09月11日 23:00:00;0
2016年09月12日 00:00:00;0
2016年09月12日 01:00:00;0
2016年09月12日 02:00:00;0
2016年09月12日 03:00:00;0
2016年09月12日 04:00:00;0
2016年09月12日 05:00:00;0
2016年09月12日 06:00:00;0
2016年09月12日 07:00:00;0
2016年09月12日 08:00:00;66
2016年09月12日 09:00:00;337
2016年09月12日 10:00:00;695
2016年09月12日 11:00:00;1041
2016年09月12日 12:00:00;1310
2016年09月12日 13:00:00;1478
2016年09月12日 14:00:00;1510
2016年09月12日 15:00:00;1426
2016年09月12日 16:00:00;1241
2016年09月12日 17:00:00;952
2016年09月12日 18:00:00;595
2016年09月12日 19:00:00;265
2016年09月12日 20:00:00;42
2016年09月12日 21:00:00;0
2016年09月12日 22:00:00;0
2016年09月12日 23:00:00;0
2016年09月13日 00:00:00;0
2016年09月13日 01:00:00;0

• Reinderien: Do your actual data have more than one day's worth of readings? Why do you replace dates in 2021 with dates in 2022?
• cebep27: Yes they have more than one day of readings, I was trying to read every hour of the day, and calculate the quantiles for each hour seperately.
• Reinderien: Do you have a link to the full data? I ask because currently you have so little sample data that your quantile spines are invisible.
• cebep27: I have edited the question in order to include more data, I hope this makes things clearer

Answer 1

You don't strictly need Seaborn for what you're doing.

Prefer to express your column names as an immutable tuple.

You should tell read_csv to do actual date parsing right off the hop. Your daterange method and your start_date and end_date should not be necessary. As much as possible, use "real" Numpy/Pandas/Matplotlib datetime support, rather than stringy manipulation (contains, etc.). There are some annoying gaps that have to be hacked around, but overall this can be dragged into a workable state.

Add PEP484 type hints.

You should not need to hard-code an entire list of percentiles (i.e. pct0.1, etc.) and this should be done in a loop. Once properly abstracted, you'll find that you aren't limited to ten quantiles - you could just as easily use 20 or more. This should be explained in your graph with a side colour bar.

Don't arange/10; use linspace instead.

Don't call quantile() once for each of your quantiles; instead call it once passing the entire array for this operation to vectorise.

I don't think it's appropriate to set xs, the horizontal axis of your plot, to be an arange over the entire dataframe: instead, you should apply grouping over the time.

Avoid calling plt.* functions when more direct and explicit ax.* and fig.* methods are available.

I find the automatic y-limits to fit somewhat better than your hard-coded 3200.

This may sound radical, but you should replace your Greens with a perceptually uniform colour map. I've shown a common one - viridis - but there are others. If you want to learn more, this is a fairly deep rabbit hole but fascinating reading material. Long story short, such a colour map is a higher-fidelity human-computer interface.

Suggested

import numpy as np
from datetime import datetime, date, time
import matplotlib.pyplot as plt
import pandas as pd
from matplotlib.cm import ScalarMappable, get_cmap
from matplotlib.dates import DateFormatter
from pandas.core.groupby import SeriesGroupBy
def load_data(
 filename: str = 'WTPV_Spring_2016.csv',
) -> pd.DataFrame:
 """
 Load the indicated CSV, parsing dates.
 """
 return pd.read_csv(
 filepath_or_buffer=filename,
 names=('Date', 'Energy'), parse_dates=[0],
 encoding='utf8', delimiter=';',
 )
def make_quantiles(
 df: pd.DataFrame,
 quantiles: np.ndarray,
) -> pd.Series:
 """
 Apply the given quantiles over time groups. The input quantiles are assumed
 to include 0 and 1 but these are not included in the calculation.
 The output is a series with a multi-level index: time, quantile. For plot
 compatibility the time is represented as a time of day on Jan 1 1970.
 """
 fake_date = date(1970, 1, 1)
 def dt_to_time(t: time) -> datetime:
 return datetime.combine(fake_date, t)
 times = df.Date.dt.time.apply(dt_to_time)
 df.set_index(times, inplace=True)
 by_time: SeriesGroupBy = df.Energy.groupby(level=0)
 bands = by_time.quantile(quantiles[1:-1])
 bands.index.names = ('Time', 'Quantile')
 return bands
def plot(quantiles: np.ndarray, bands: pd.Series) -> plt.Figure:
 """
 Plot the given quantile bands as filled regions with an associated colour
 bar. The colour bar shows only the first half of the quantile range; the
 second half is symmetric and implied.
 """
 fig, ax = plt.subplots()
 ax.set_title('WT & PV Spring 2016')
 ax.set_xlabel('Time')
 ax.set_ylabel('Energy (MWh)')
 fig.autofmt_xdate()
 ax.xaxis.set_major_formatter(DateFormatter('%H:%M'))
 map = ScalarMappable(
 cmap=get_cmap('viridis'),
 norm=plt.Normalize(vmin=0, vmax=0.5),
 )
 ax.set_facecolor(map.to_rgba(0))
 fig.colorbar(map, label='quantile')
 counterposed_quantiles = np.vstack((
 quantiles[1: len(quantiles)//2],
 quantiles[-2: len(quantiles)//2: -1],
 )).T
 for q0, q1 in counterposed_quantiles:
 y0 = bands.loc[:, q0]
 x = y0.index
 y1 = bands.loc[x, q1]
 ax.fill_between(x, y0, y1, color=map.to_rgba(q0))
 q = 0.5
 y = bands.loc[:, q]
 x = y.index
 ax.plot(x, y, color=map.to_rgba(q))
 return fig
def main() -> None:
 quantiles = np.linspace(0, 1, 41)
 data = load_data()
 bands = make_quantiles(data, quantiles)
 plot(quantiles, bands)
 plt.show()
if __name__ == '__main__':
 main()

graph