On Friday, August 26, 1983, Bilbao was celebrating its Aste Nagusia or Great Week, the main annual festivity in the city, when it and other municipalities of the Basque Country, Burgos, and Cantabria suffered devastating flooding due to heavy rains. In 24 hours, the volume of water registered 600 liters per square meter. Across all the affected areas, the weather service recorded 1.5 billion tons of water. In areas of Bilbao, the water reached a height of 5 meters (15 feet). Transportation, electricity and gas services, drinking water, food, telephone, and many other basic services were severely affected. 32 people died in Biscay, 4 people died in Cantabria, 2 people died in Alava, and 2 people died Burgos. 5 more people went missing.
How often does it rain 50 mm in 1 day? What about 100 mm in 1 day? How big is a “once-in-a-century event”?
Let’s examine Bilbao’s daily rainfall (mm), between 1947 to 2021
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns
sns.set_theme(style="ticks", font_scale=1.5)
import urllib.request
from matplotlib.dates import DateFormatter
import matplotlib.dates as mdates
import altair as alt
alt.data_transformers.disable_max_rows()
from scipy.stats import genextremedef download_data(station_name, station_code):
url_daily = 'https://www.ncei.noaa.gov/data/global-historical-climatology-network-daily/access/'
url_monthly = 'https://www.ncei.noaa.gov/data/gsom/access/'
# download daily data - uncomment to make this work
urllib.request.urlretrieve(url_daily + station_code + '.csv',
station_name + '_daily.csv')
# download monthly data
urllib.request.urlretrieve(url_monthly + station_code + '.csv',
station_name + '_monthly.csv')
download_data('BILBAO', 'SPE00120611')df = pd.read_csv('BILBAO_daily.csv',
sep=",",
parse_dates=['DATE'],
index_col='DATE')
# IMPORTANT!! daily precipitation data is in tenths of mm, divide by 10 to get it in mm.
df['PRCP'] = df['PRCP'] / 10
fig, ax = plt.subplots(figsize=(10,7))
ax.plot(df['PRCP'])
ax.set(xlabel="date",
ylabel="daily rainfall (mm)",
title="Bilbao, Spain, 1947--2021"
)
ax.annotate("26 August 1983",
xy=(pd.to_datetime('1983-08-26'), 2500), xycoords='data',
xytext=(0.7, 0.95), textcoords='axes fraction',
fontsize=16, va="center",
arrowprops=dict(facecolor='black', shrink=0.05));
On the week of 22-28 August 1983, Bilbao’s weather station measured 45 cm of rainfall!
fig, ax = plt.subplots(figsize=(10,7))
one_week = df.loc['1983-08-22':'1983-08-28', 'PRCP']
bars = ax.bar(one_week.index, one_week)
ax.set_xlabel("date")
ax.set_ylabel("daily rainfall (mm)")
ax.set_title("Bilbao, Spain, August 1983")
# write daily rainfall
for i in range(len(one_week)):
ax.text(one_week.index[i], one_week.iloc[i], f"{one_week.iloc[i]:.0f}", ha="center", fontsize=16)
# ax.text(one_week.index[i], one_week[i], f"{one_week[i]:.0f}", ha="center", fontsize=16);
ax.text(0.1, 0.8, f"Total rainfall during this week:\n{one_week.sum():.0f} mm",
transform=ax.transAxes, fontsize=16)
# Define the date format
# https://docs.python.org/2/library/datetime.html#strftime-and-strptime-behavior
date_form = DateFormatter("%b-%d")
ax.xaxis.set_major_formatter(date_form)
# Ensure a major tick for each day using (interval=1)
# https://matplotlib.org/stable/api/dates_api.html#date-tickers
ax.xaxis.set_major_locator(mdates.DayLocator(interval=1))
Let’s analyze this data and find out how rare such events are. First we need to find the annual maximum for each hydrological year. Do we see any seasonal patterns with our eyes? Play with the widget below.
# Altair only recognizes column data; it ignores index values.
# You can plot the index data by first resetting the index
# I know that I've just made 'DATE' the index, but I want to have this here nonetheless so I can refer to this in the future
df_new = df.reset_index()
source = df_new[['DATE', 'PRCP']]
brush = alt.selection_interval(encodings=['x'])
base = alt.Chart(source).mark_line().encode(
x = 'DATE:T',
y = 'PRCP:Q'
).properties(
width=600,
height=200
)
upper = base.encode(
alt.X('DATE:T', scale=alt.Scale(domain=brush)),
alt.Y('PRCP:Q', scale=alt.Scale(domain=(0,100)))
)
lower = base.properties(
height=60
).add_params(brush)
alt.vconcat(upper, lower)