{ "cells": [ { "cell_type": "markdown", "id": "0", "metadata": {}, "source": [ "# Lesson: timeseries analysis and extreme values" ] }, { "cell_type": "markdown", "id": "1", "metadata": {}, "source": [ "Over the past two weeks, you’ve learned how to work with *NetCDF* files, exploring their complex multidimensional structures, plotting global maps, and summarizing climate projection datasets. Today, we’re shifting gears slightly. Instead of working with gridded data, we’ll (temporarily) return to something more familiar perhaps: time series in tabular form. While this may seem simpler at first, handling time series data comes with its own set of challenges!\n", "\n", "To get started, let’s import the libraries we’ll be using. This week, we’re introducing two key newcomers:\n", "- [pandas](https://pandas.pydata.org/), replacing xarray for now: pandas is the backbone of data analysis in Python. It actually predates xarray and serves as the foundation for many of its core functions. If you’ve grown comfortable with xarray, you’ll likely find pandas’ syntax quite familiar.\n", "- [SciPy](https://scipy.org/): the Swiss Army knife of statistical analysis in Python. Whether you’re performing interpolation, fitting models, or conducting statistical tests, SciPy is an essential tool for scientific computing.\n", "\n", "*Copyright notice: parts of this class is inpired by the excellent [climatematch tutorials](https://comptools.climatematch.io/tutorials/W2D3_ExtremesandVariability/student/W2D3_Tutorial2.html). I really recommend you to check them out!*\n", "\n", "
\n", " This a dense section. I know! This is why you are receiving this before the consolidation week - I encourage you to take your time, and spend next week to revisit the fundamentals of both pandas and xarray. Using these resources, but also my recommended resources. I recommend Climatematch and Project Pythia.\n", "
" ] }, { "cell_type": "code", "execution_count": null, "id": "2", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import pandas as pd # This is new\n", "from scipy import stats # This also" ] }, { "cell_type": "markdown", "id": "3", "metadata": {}, "source": [ "## Daily meteorological observations at Heathrow" ] }, { "cell_type": "markdown", "id": "4", "metadata": {}, "source": [ "Go to the [download page](https://fabienmaussion.info/climate_risks/ready/03-download.html). Download the Heathrow meteorological data and put it in your `data` folder. To keep things clean, I decided to put mine in a `csv` folder within `data`. \"CSV\" stands for comma-separated values, and I'm sure you already heard of it? If not, a quick web search is in order.\n", "\n", "The `gsod-heathrow.csv` data file is a text file I downloaded from the \"[Global Surface Summary of the Day](https://www.ncei.noaa.gov/access/metadata/landing-page/bin/iso?id=gov.noaa.ncdc:C00516)\" database, managed by NOAA. Since we can't trust the current US adminstration to keep this infrastructure running for the duration of the class, I mirrored the GSOD readme on our server - you'll find it [here](https://cluster.klima.uni-bremen.de/~fmaussion/teaching/qcr/csv/readme.txt).\n", "\n", "Let's open the file and display its content:" ] }, { "cell_type": "code", "execution_count": null, "id": "5", "metadata": {}, "outputs": [], "source": [ "dfh = pd.read_csv('../data/csv/gsod-heathrow.csv', # Read the CSV file using Pandas\n", " index_col=1, # what's this doing?\n", " parse_dates=True # and this?\n", " )" ] }, { "cell_type": "code", "execution_count": null, "id": "6", "metadata": {}, "outputs": [], "source": [ "dfh" ] }, { "cell_type": "markdown", "id": "7", "metadata": {}, "source": [ "Like xarray objects, pandas objects are tailored for data exploration. In xarray, we manipulated objects called \"Datasets\" and \"DataArrays\". In pandas, the equivalents are \"Dataframes\" and \"Series\":" ] }, { "cell_type": "code", "execution_count": null, "id": "8", "metadata": {}, "outputs": [], "source": [ "type(dfh)" ] }, { "cell_type": "code", "execution_count": null, "id": "9", "metadata": {}, "outputs": [], "source": [ "type(dfh['TEMP'])" ] }, { "cell_type": "code", "execution_count": null, "id": "10", "metadata": {}, "outputs": [], "source": [ "dfh['TEMP']" ] }, { "cell_type": "markdown", "id": "11", "metadata": {}, "source": [ "All \"Series\" in a \"Dataframe\" share the same coordinate, the \"index\":" ] }, { "cell_type": "code", "execution_count": null, "id": "12", "metadata": {}, "outputs": [], "source": [ "dfh.index" ] }, { "cell_type": "markdown", "id": "13", "metadata": {}, "source": [ "Let's plot the data:" ] }, { "cell_type": "code", "execution_count": null, "id": "14", "metadata": {}, "outputs": [], "source": [ "dfh['TEMP'].plot();" ] }, { "cell_type": "markdown", "id": "15", "metadata": {}, "source": [ "**Q: spend some time exploring the dataset. Plot one or two variables, including `MAX` and `PRCP`. What can you notice?**" ] }, { "cell_type": "code", "execution_count": null, "id": "16", "metadata": {}, "outputs": [], "source": [ "# Your answer here" ] }, { "cell_type": "markdown", "id": "17", "metadata": {}, "source": [ "Welcome to the messy world of real-world observations! Historical data is always incomplete or imperfect, and we have to deal with it. This file, by the way, is considerably easier to deal with and has much less missing data than the Bristol or Cardiff stations data.\n", "\n", "This is an unfortunate bit of working with in-situ observations and non-model data: they always need a bit of \"cleaning\", and I won't lie: it's the boring bit and we have to go through it.\n", "\n", "## Data preparation\n", "\n", "
\n", " Data preparation is an annoying but indispensable process of any data analysis workflow. In your real-world job as an environmental consultant or analyst, you'll find that a lot of the work has to do with data preparation. In the few cells below I guide you through the process I did to get your data ready for analysis. If you really can't be bothered, fast forward to the next section. If you want to learn more, stay with me!\n", "
\n", "\n", "### Missing values\n", "\n", "Let's start by noting that the [README](https://cluster.klima.uni-bremen.de/~fmaussion/teaching/qcr/csv/readme.txt) mentions missing values. GSOD chose to replace them with specific numbers, which are messing with our plots and statistics. Let's replace these with `np.nan`, which is short for \"Not a Number\":" ] }, { "cell_type": "code", "execution_count": null, "id": "18", "metadata": {}, "outputs": [], "source": [ "dfh = dfh.replace([99.99, 999.9, 9999.9], np.nan)" ] }, { "cell_type": "markdown", "id": "19", "metadata": {}, "source": [ "**Q: compare your new plot of `MAX` and `PRCP` with the ones you've plotted above. Notice the difference?** \n", "\n", "Note that while the missing values are not messing with the plot anymore, they are still missing. Somewhere in the plot, there are \"holes\". We can actually count missing values each year with a simple trick:" ] }, { "cell_type": "code", "execution_count": null, "id": "20", "metadata": {}, "outputs": [], "source": [ "(dfh['PRCP'].isnull().resample('YS').mean() * 100).plot(); # We'll explain this below\n", "\n", "plt.title('Missing values per calendar year for the variable PRCP');\n", "plt.ylabel('Missing values (%)');" ] }, { "cell_type": "markdown", "id": "21", "metadata": {}, "source": [ "OK - let's take a step back and see what we've done. \n", "\n", "**Q: step by step, decompose the one-liner above. First, check what `dfh['PRCP'].isnull()` does. Then remember when you've see `.resample()` last time, and note that it works exactly the same!**" ] }, { "cell_type": "code", "execution_count": null, "id": "22", "metadata": {}, "outputs": [], "source": [ "# Your answer here" ] }, { "cell_type": "markdown", "id": "23", "metadata": {}, "source": [ "The final trick is to note that averaging booleans (\"`False`\" and \"`True`\") is just averaging 0s and 1s:" ] }, { "cell_type": "code", "execution_count": null, "id": "24", "metadata": {}, "outputs": [], "source": [ "np.mean([False, True, False])" ] }, { "cell_type": "markdown", "id": "25", "metadata": {}, "source": [ "If the above is not clear yet, please ask!\n", "\n", "Now, let's analyse the plot. We have about 30% of missing PRCP data in the first year of the series, followed by a blank period in the 1960s. To complicate the matter a bit, it turns out that in that period we have no data at all: no timestamps, and no measurements. This is actually hiding the fact that we have 100% data missing here!\n", "\n", "One more step of data preparation is to notice that we expect a certain format for these files. Since we know the frequency of the measurements (daily) as well as the start and end dates, we can use pandas to create a timestamp index of the length (in days) of the full time period:" ] }, { "cell_type": "code", "execution_count": null, "id": "26", "metadata": {}, "outputs": [], "source": [ "full_dates = pd.date_range(dfh.index[0], # Initial date\n", " dfh.index[-1], # Last date\n", " freq='D' # Frequency (daily)\n", " )\n", "\n", "len(dfh), len(full_dates) # Compare the lengths" ] }, { "cell_type": "markdown", "id": "27", "metadata": {}, "source": [ "We see that we are missing approx. 5000 records in the Dataframe. Fortunately, it's fairly easy to format the data appropriately:" ] }, { "cell_type": "code", "execution_count": null, "id": "28", "metadata": {}, "outputs": [], "source": [ "dfh = dfh.reindex(full_dates)" ] }, { "cell_type": "markdown", "id": "29", "metadata": {}, "source": [ "Let's reevaluate our missing values plot from before:" ] }, { "cell_type": "code", "execution_count": null, "id": "30", "metadata": {}, "outputs": [], "source": [ "(dfh['PRCP'].isnull().resample('YS').mean() * 100).plot(); # We'll explain this below\n", "\n", "plt.title('Missing values per calendar year for the variable PRCP');\n", "plt.ylabel('Missing values (%)');" ] }, { "cell_type": "markdown", "id": "31", "metadata": {}, "source": [ "That's more like it. For over a decade, 100% of the data is missing. This makes the record a bit useless before 1973. Also, since 2025 is not complete yet, we might as well get rid of it! \n", "\n", "We now select the time period we'll focus on for further analysis:" ] }, { "cell_type": "code", "execution_count": null, "id": "32", "metadata": {}, "outputs": [], "source": [ "dfh = dfh.loc['1973':'2024']" ] }, { "cell_type": "markdown", "id": "33", "metadata": {}, "source": [ "### Outliers / artifacts" ] }, { "cell_type": "markdown", "id": "34", "metadata": {}, "source": [ "Now, we're done! Or, are we really? Let's have a final look at precipitation:" ] }, { "cell_type": "code", "execution_count": null, "id": "35", "metadata": {}, "outputs": [], "source": [ "dfh['PRCP'].plot();" ] }, { "cell_type": "code", "execution_count": null, "id": "36", "metadata": {}, "outputs": [], "source": [ "(dfh['PRCP'] > 10).sum()" ] }, { "cell_type": "markdown", "id": "37", "metadata": {}, "source": [ "Mh, there are 4 days with over 10 inches (254 millimeters) of recorded rain in 24 hours. While this is not strictly impossible, these values would crush many confirmed UK record precipitation at weather stations ([which are around 50-60 mm daily rain](https://rmets.onlinelibrary.wiley.com/doi/10.1002/wea.3910)) and would represent more than a quarter of the annual rain in a single day. While I can' be 100% sure, I'm quite confifent that this represent true outliers and should be removed.\n", "\n", "Let's do that:" ] }, { "cell_type": "code", "execution_count": null, "id": "38", "metadata": {}, "outputs": [], "source": [ "dfh['PRCP'] = dfh['PRCP'].where(dfh['PRCP'] < 10) # Replace values above 10 inches with NaN\n", "\n", "dfh['PRCP'].plot();" ] }, { "cell_type": "markdown", "id": "39", "metadata": {}, "source": [ "This just looks more plausible now. The remaining most visible extreme value (4 inches, i.e. ~ 100 mm in a day) is still very high but more plausible. It happened on 1994-02-24, and might or might not be correct:" ] }, { "cell_type": "code", "execution_count": null, "id": "40", "metadata": {}, "outputs": [], "source": [ "dfh['PRCP'].max(), dfh['PRCP'].idxmax()" ] }, { "cell_type": "markdown", "id": "41", "metadata": {}, "source": [ "### Units" ] }, { "cell_type": "markdown", "id": "42", "metadata": {}, "source": [ "The [README](https://cluster.klima.uni-bremen.de/~fmaussion/teaching/qcr/csv/readme.txt) tells us the units of our variables. They are (horror!) in US units. Fortunately, it is easy to convert)" ] }, { "cell_type": "code", "execution_count": null, "id": "43", "metadata": {}, "outputs": [], "source": [ "dfh['TEMP_c'] = (dfh['TEMP'] - 32) * 5/9\n", "dfh['TEMP_c'].plot();" ] }, { "cell_type": "markdown", "id": "44", "metadata": {}, "source": [ "Good job! Now we’re ready to move forward with the analysis. 🚀\n", "\n", "
\n", " Key Steps in Data Preparation:\n", " \n", "
\n" ] }, { "cell_type": "markdown", "id": "45", "metadata": {}, "source": [ "## Data exploration and descriptive statistics " ] }, { "cell_type": "markdown", "id": "46", "metadata": {}, "source": [ "Let's repeat the data preparation steps in a single cell to make sure we are on the same page:" ] }, { "cell_type": "code", "execution_count": null, "id": "47", "metadata": {}, "outputs": [], "source": [ "# Read the data\n", "dfh = pd.read_csv('../data/csv/gsod-heathrow.csv', index_col=1, parse_dates=True)\n", "\n", "# Missing values\n", "dfh = dfh.replace([99.99, 999.9, 9999.9], np.nan)\n", "dfh = dfh.reindex(pd.date_range(dfh.index[0], dfh.index[-1], freq='D'))\n", "\n", "# Period selection\n", "dfh = dfh.loc['1973':'2024']\n", "\n", "# Outlier filtering\n", "dfh['PRCP'] = dfh['PRCP'].where(dfh['PRCP'] < 10) \n", "\n", "# Keep the variables we want\n", "dfh = dfh[['TEMP', 'MAX', 'MIN', 'PRCP']].copy()\n", "\n", "# Convert units\n", "dfh[['TEMP', 'MAX', 'MIN']] = (dfh[['TEMP', 'MAX', 'MIN']] - 32) * 5/9\n", "dfh['PRCP'] = dfh['PRCP'] * 25.4" ] }, { "cell_type": "code", "execution_count": null, "id": "48", "metadata": {}, "outputs": [], "source": [ "dfh" ] }, { "cell_type": "markdown", "id": "49", "metadata": {}, "source": [ "### Descriptive statistics" ] }, { "cell_type": "markdown", "id": "50", "metadata": {}, "source": [ "Like xarray, pandas allows to compute simple statics, such as the average temperature and standard deviation of daily temperature:" ] }, { "cell_type": "code", "execution_count": null, "id": "51", "metadata": {}, "outputs": [], "source": [ "mean = dfh['TEMP'].mean()\n", "avg = dfh['TEMP'].std()\n", "\n", "f\"The temp average is {mean:.2f}°C and the daily standard deviation is {avg:.2f}°C\"" ] }, { "cell_type": "markdown", "id": "52", "metadata": {}, "source": [ "It also allows to compute the quantiles of the data:" ] }, { "cell_type": "code", "execution_count": null, "id": "53", "metadata": {}, "outputs": [], "source": [ "median = dfh['TEMP'].median()\n", "q75 = dfh['TEMP'].quantile(0.75)\n", "q25 = dfh['TEMP'].quantile(0.25)\n", "\n", "f\"The temp median is {median:.2f}°C and the interquartile range is [{q25:.2f}-{q75:.2f}]°C\"" ] }, { "cell_type": "markdown", "id": "54", "metadata": {}, "source": [ "**E: compute the mean, median, and 10 and 90% quantiles of daily precipitation.**" ] }, { "cell_type": "code", "execution_count": null, "id": "55", "metadata": {}, "outputs": [], "source": [ "# Your answer here" ] }, { "cell_type": "markdown", "id": "56", "metadata": {}, "source": [ "### Period selection" ] }, { "cell_type": "markdown", "id": "57", "metadata": {}, "source": [ "It's equally easy to select only parts of the index. For example, to plot the daily temperatures over the year 1987:" ] }, { "cell_type": "code", "execution_count": null, "id": "58", "metadata": {}, "outputs": [], "source": [ "dfh_87 = dfh.loc['1987'] # Similar to \".sel()\" in xarray" ] }, { "cell_type": "code", "execution_count": null, "id": "59", "metadata": {}, "outputs": [], "source": [ "dfh_87['TEMP'].plot(color='k', linewidth=2, label='Daily avg')\n", "dfh_87['MAX'].plot(color='C3', label='Daily max')\n", "dfh_87['MIN'].plot(color='C0', label='Daily min')\n", "plt.legend(); plt.title('Daily temperature statistics at LHR in 1987'); plt.ylabel('°C');" ] }, { "cell_type": "markdown", "id": "60", "metadata": {}, "source": [ "**E: Now repeat the plot with the year 2022 (the hottest day on record), and zoom in the months of June to September 2022.** *Hint: to select a time period, you will need to use `.loc[date1:date2]`*" ] }, { "cell_type": "code", "execution_count": null, "id": "61", "metadata": {}, "outputs": [], "source": [ "# Your answer here" ] }, { "cell_type": "markdown", "id": "62", "metadata": {}, "source": [ "### Temporal aggregation: resampling" ] }, { "cell_type": "markdown", "id": "63", "metadata": {}, "source": [ "Like xarray, pandas allows to compute statistics for sepcific time average. For example monthly:" ] }, { "cell_type": "code", "execution_count": null, "id": "64", "metadata": {}, "outputs": [], "source": [ "dfh.resample('MS').mean()['TEMP'].plot();" ] }, { "cell_type": "markdown", "id": "65", "metadata": {}, "source": [ "Or annual:" ] }, { "cell_type": "code", "execution_count": null, "id": "66", "metadata": {}, "outputs": [], "source": [ "dfh.resample('YS').mean()['TEMP'].plot();" ] }, { "cell_type": "markdown", "id": "67", "metadata": {}, "source": [ "**E: compute the difference between the average temperature over the last decade (2015-2024) and the first decade in the record (1973-1982). Repeat with maximum, minimum temperature and precipitation.** (*Hint: you can avoid repetition by noticing that one can do statistics on entire dataframes.*\n", "\n", "**Q: Did the maximum or mean temperatures change more during the two period? Was it on average wetter or drier over the past decade?**" ] }, { "cell_type": "code", "execution_count": null, "id": "68", "metadata": {}, "outputs": [], "source": [ "# Your answer here" ] }, { "cell_type": "markdown", "id": "69", "metadata": {}, "source": [ "### Temporal aggregation: rolling window " ] }, { "cell_type": "markdown", "id": "70", "metadata": {}, "source": [ "\"Rolling window\" statistics is another way to compute temporal averages. As the name suggests, for each timestamp, statistics are computed over a window before or centered on that timestamp.\n", "\n", "Let's compute, for example, the decadal average of temperature change over the annual timeseries:" ] }, { "cell_type": "code", "execution_count": null, "id": "71", "metadata": {}, "outputs": [], "source": [ "at.rolling(11, center=True).mean()" ] }, { "cell_type": "code", "execution_count": null, "id": "72", "metadata": {}, "outputs": [], "source": [ "at = dfh['TEMP'].resample('YS').mean()\n", "roll = at.rolling(11, center=True).mean()\n", "\n", "at.plot(c='k', label='Annual temperature');\n", "roll.plot(c='C0', linewidth=3, label='11-year rolling mean');\n", "plt.legend();" ] }, { "cell_type": "markdown", "id": "73", "metadata": {}, "source": [ "**Q: interpret the plot. Why do you think the rolling timeseries are not covering the full period? Now repet the plot without the `center=True` argument.**" ] }, { "cell_type": "code", "execution_count": null, "id": "74", "metadata": {}, "outputs": [], "source": [ "# Your answer here" ] }, { "cell_type": "markdown", "id": "75", "metadata": {}, "source": [ "### Selection by condition" ] }, { "cell_type": "markdown", "id": "76", "metadata": {}, "source": [ "You may find it useful to select parts of the data only, for example all August days, or all days above 0 mm precipitation. This is done with comparison operators (e.g. `==`, `>`, `<=`, etc.).\n", "\n", "For example, let's select all days where precip is non zero:" ] }, { "cell_type": "code", "execution_count": null, "id": "77", "metadata": {}, "outputs": [], "source": [ "dfhp = dfh.loc[dfh['PRCP'] > 0]" ] }, { "cell_type": "markdown", "id": "78", "metadata": {}, "source": [ "**Q: do you understand what `dfh['PRCP'] > 0` is doing? Perhaps it's worth exploring it step by step, for example by printing the output of `dfh['PRCP']`, then remember what `>` does.**" ] }, { "cell_type": "code", "execution_count": null, "id": "79", "metadata": {}, "outputs": [], "source": [ "# Your answer here" ] }, { "cell_type": "markdown", "id": "80", "metadata": {}, "source": [ "**E: count the number of rainy days at Heathrow. Compute the ratio of rainy days over all days. Now repeat with heavy rain days, with more than > 10 mm precipitation per day. What's the ratio of heavy rain days over all rainy days?**" ] }, { "cell_type": "code", "execution_count": null, "id": "81", "metadata": {}, "outputs": [], "source": [ "# Your answer here" ] }, { "cell_type": "markdown", "id": "82", "metadata": {}, "source": [ "We can use the same strategy to pick days withing a given month:" ] }, { "cell_type": "code", "execution_count": null, "id": "83", "metadata": {}, "outputs": [], "source": [ "dfh_aug = dfh.loc[dfh.index.month == 8]" ] }, { "cell_type": "markdown", "id": "84", "metadata": {}, "source": [ "**Q: do you understand what `dfh.index.month == 8` is doing? Perhaps it's worth exploring it step by step, for example by printing the output of `dfh.index.month`, then remember what `==` does.**" ] }, { "cell_type": "code", "execution_count": null, "id": "85", "metadata": {}, "outputs": [], "source": [ "# Your answer here" ] }, { "cell_type": "markdown", "id": "86", "metadata": {}, "source": [ "**E: OK, now compute the average daily temperature in August. Repeat for the month of January.**" ] }, { "cell_type": "code", "execution_count": null, "id": "87", "metadata": {}, "outputs": [], "source": [ "# Your answer here" ] }, { "cell_type": "markdown", "id": "88", "metadata": {}, "source": [ "Bonus: selection for multiple months:" ] }, { "cell_type": "code", "execution_count": null, "id": "89", "metadata": {}, "outputs": [], "source": [ "dfh_jja = dfh.loc[(dfh.index.month == 6) | (dfh.index.month == 7) | (dfh.index.month == 8)] # OK\n", "dfh_jja = dfh.loc[dfh.index.month.isin([6, 7, 8])] # better" ] }, { "cell_type": "markdown", "id": "90", "metadata": {}, "source": [ "### Selection by category: Groupby" ] }, { "cell_type": "markdown", "id": "91", "metadata": {}, "source": [ "What we've done with two months above can be generalized for all various categories. For example, it is possible to compute summary statistics organised by monthss:" ] }, { "cell_type": "code", "execution_count": null, "id": "92", "metadata": {}, "outputs": [], "source": [ "dfhm = dfh.groupby(dfh.index.month).mean()\n", "dfhm['TEMP'].plot();" ] }, { "cell_type": "markdown", "id": "93", "metadata": {}, "source": [ "**Q: do you understand what `dfh.groupby(dfh.index.month)` is doing? Perhaps it's worth exploring it step by step, for example by printing the output of `dfh.index.month`, then remember what `groupby` does, and what `.mean()` does.**" ] }, { "cell_type": "code", "execution_count": null, "id": "94", "metadata": {}, "outputs": [], "source": [ "# Your answer here" ] }, { "cell_type": "markdown", "id": "95", "metadata": {}, "source": [ "**Q: check that the monthly averages for January and August you computed above match the plot.**" ] }, { "cell_type": "markdown", "id": "96", "metadata": {}, "source": [ "**E: now plot (and interpret) the monthly average precipitation chart (instead of temperature):**" ] }, { "cell_type": "code", "execution_count": null, "id": "97", "metadata": {}, "outputs": [], "source": [ "# Your answer here" ] }, { "cell_type": "markdown", "id": "98", "metadata": {}, "source": [ "Groupby is very useful. We could also group by rain categories, but this is perhaps for another time. We can however easily adapt \"groupby\" to use days instead of months:" ] }, { "cell_type": "code", "execution_count": null, "id": "99", "metadata": {}, "outputs": [], "source": [ "dfhd = dfh.groupby(dfh.index.dayofyear) # We have only grouped! No statistics yet\n", "\n", "dfhd['TEMP'].mean().plot(color='k');\n", "dfhd['MAX'].mean().plot(color='C3');\n", "dfhd['MAX'].max().plot(style='+', color='C3', alpha=0.5);" ] }, { "cell_type": "markdown", "id": "100", "metadata": {}, "source": [ "**E: Add proper labels to the curves above. Do you understand what's plotted there? Now also add the equivalent for the MIN temperatures.** *Hint: the goal is to reproduce a plot similar to the [Meteoblue climate diagrams](https://www.meteoblue.com/en/weather/historyclimate/climateobserved/bristol_united-kingdom_2654675).*" ] }, { "cell_type": "code", "execution_count": null, "id": "101", "metadata": {}, "outputs": [], "source": [ "# Your answer here" ] }, { "cell_type": "markdown", "id": "102", "metadata": {}, "source": [ "## Histograms " ] }, { "cell_type": "markdown", "id": "103", "metadata": {}, "source": [ "Let's start by selecting all daily values in the month of July:" ] }, { "cell_type": "code", "execution_count": null, "id": "104", "metadata": {}, "outputs": [], "source": [ "dfh_jul = dfh.loc[dfh.index.month == 7]" ] }, { "cell_type": "markdown", "id": "105", "metadata": {}, "source": [ "These represent 1612 values altogether. We are talking about a statistical sample of a population. The histogram of this sample can be plotted with pandas easily, here with the example TEMP:" ] }, { "cell_type": "code", "execution_count": null, "id": "106", "metadata": {}, "outputs": [], "source": [ "sample = dfh_jul['TEMP']\n", "sample.plot.hist(edgecolor='k'); \n", "plt.title('Histogram of daily temperature in July');" ] }, { "cell_type": "markdown", "id": "107", "metadata": {}, "source": [ "The Y axis represents the number of events by bins in the sample (\"frequency\" is a bit misleading). Usually we like to work in density (percentages). \n", "\n", "Pandas decides on the bins automatically. It's possible to impose a number of bins, or bin size, or bin ranges:" ] }, { "cell_type": "code", "execution_count": null, "id": "108", "metadata": {}, "outputs": [], "source": [ "sample.plot.hist(edgecolor='k', density=True, bins=51); \n", "plt.title('Histogram of daily temperature in July');" ] }, { "cell_type": "code", "execution_count": null, "id": "109", "metadata": {}, "outputs": [], "source": [ "sample.plot.hist(edgecolor='k', density=True, \n", " bins=np.arange(10, 31, 1)); # Note: this implies you know the range of values you want to cover\n", "plt.title('Histogram of daily temperature in July');" ] }, { "cell_type": "markdown", "id": "110", "metadata": {}, "source": [ "The histogram values can also be computed for further analysis instead of plotting:" ] }, { "cell_type": "code", "execution_count": null, "id": "111", "metadata": {}, "outputs": [], "source": [ "# Compute histogram values (raw counts)\n", "bins = np.arange(10, 31, 2)\n", "hist_values, bin_edges = np.histogram(sample, bins=bins)" ] }, { "cell_type": "markdown", "id": "112", "metadata": {}, "source": [ "**E: Explore the `hist_values` and `bin_edges` variables. Can you understand what is done?**" ] }, { "cell_type": "code", "execution_count": null, "id": "113", "metadata": {}, "outputs": [], "source": [ "# Your answer here" ] }, { "cell_type": "markdown", "id": "114", "metadata": {}, "source": [ "**E: repeat the analysis for precipitation, first by letting pandas decide on the bins, then by defining bins yourself starting at 0 mm up to 45 mm (plot A). Now do another plot (plot B) where bins start at 1 mm instead of 0 mm. Discuss the differences in the plot.**" ] }, { "cell_type": "code", "execution_count": null, "id": "115", "metadata": {}, "outputs": [], "source": [ "# Your answer here" ] }, { "cell_type": "markdown", "id": "116", "metadata": {}, "source": [ "## Probability density functions" ] }, { "cell_type": "markdown", "id": "117", "metadata": {}, "source": [ "A population can be approximated mathematically with probabilty density functions (PDFs). The most common one is the gaussian PDF, which is defined by only two parameters, the mean and standard deviation." ] }, { "cell_type": "code", "execution_count": null, "id": "118", "metadata": {}, "outputs": [], "source": [ "mean_temp = sample.mean()\n", "std_temp = sample.std()\n", "\n", "ax = sample.plot.hist(edgecolor='k', density=True, bins=np.arange(10, 31, 1));\n", "plt.title('Histogram of daily temperature in July');\n", "\n", "# add in vertical line at mean\n", "ylim = ax.get_ylim()\n", "ax.vlines(mean_temp, ymin=ylim[0], ymax=ylim[1], color=\"C3\", lw=3, label=\"mean\")\n", "\n", "# add PDF\n", "x_pdf = np.arange(7, 31, 0.1)\n", "y_pdf = stats.norm.pdf(x_pdf, mean_temp, std_temp)\n", "ax.plot(x_pdf, y_pdf, c=\"k\", lw=3, label='Gaussian PDF');\n", "plt.legend();" ] }, { "cell_type": "markdown", "id": "119", "metadata": {}, "source": [ "The gaussian PDF approximate the data quite well, but there are also quite some deviations from the Gaussian. Can you spot them?\n", "\n", "**Now is the time to refresh your knowledge about normality tests and Q-Q plots.** I wont cover the fundamentals here, (we're not that much interested in this topic here), but the Climatematch tutorials offer excellent code examples and explanations: [climatematch Extremes and Variability](https://comptools.climatematch.io/tutorials/W2D3_ExtremesandVariability/chapter_title.html).\n", "\n", "For now, let's assume that a Gaussian is a good approximation of our population. From the histogram and the PDF, we can compute the probability of occurrence of a certain even. For example, the probability of a day to be hotter than 25°C on average.\n", "\n", "The *empirical* probability this event is simply given by the data. We count the number of events and evaluate their occurence in the past:" ] }, { "cell_type": "code", "execution_count": null, "id": "120", "metadata": {}, "outputs": [], "source": [ "n_events = np.sum(sample > 25) # Equivalent: len(sample.loc[sample > 25])\n", "n_sample = len(sample)\n", "f'Empirical probability of daily avg temperature > 25°C: {n_events / (n_sample + 1) * 100:.1f}%'" ] }, { "cell_type": "markdown", "id": "121", "metadata": {}, "source": [ "Note that we used (n + 1) to compute the probability. The reason for this will be explained below. For the theoretical probability based on the approximated (fitted) Gaussian distribution, we have to refer to the cumulative PDF. Remember?" ] }, { "cell_type": "code", "execution_count": null, "id": "122", "metadata": {}, "outputs": [], "source": [ "cum_y_pdf = np.cumsum(y_pdf) / np.sum(y_pdf)\n", "plt.plot(x_pdf, cum_y_pdf, c=\"k\", lw=3, label='Gaussian cumulative PDF');\n", "plt.legend(); plt.grid();" ] }, { "cell_type": "markdown", "id": "123", "metadata": {}, "source": [ "The probability of exceedance is 1 - the value of the cumulative PDF at 25°C. Can you try to see it from the plot? It's quite hard for small probabilities isn't it? Therefore, scipy gives us a way to compute it with cenrtainty:" ] }, { "cell_type": "code", "execution_count": null, "id": "124", "metadata": {}, "outputs": [], "source": [ "prob_exceedance = 1 - stats.norm.cdf(25, loc=mean_temp, scale=std_temp)\n", "f'Theoretical probability of daily avg temperature > 25°C: {prob_exceedance * 100:.1f}%'" ] }, { "cell_type": "markdown", "id": "125", "metadata": {}, "source": [ "**Q: discuss the differences between the empirical and theoretical probability of exceedance by looking at the histogram plot overlaid with the PDF. Does this make sense to you?**" ] }, { "cell_type": "code", "execution_count": null, "id": "126", "metadata": {}, "outputs": [], "source": [ "# Your answer here" ] }, { "cell_type": "markdown", "id": "127", "metadata": {}, "source": [ "## Empirical return levels" ] }, { "cell_type": "markdown", "id": "128", "metadata": {}, "source": [ "Let's continue our exploration of extreme events. We could for example have a look at the maximum recorded temperature each year.\n", "\n", "We can use resample for this:" ] }, { "cell_type": "code", "execution_count": null, "id": "129", "metadata": {}, "outputs": [], "source": [ "maxt = dfh['MAX'].resample('YS').max() # .max() is the important bit here!\n", "maxt.plot(style='-o');\n", "plt.grid(); plt.title('Annual max temperature at Heathrow'); plt.ylabel('°C');" ] }, { "cell_type": "markdown", "id": "130", "metadata": {}, "source": [ "The sample size is now much smaller! Just the number of years in the time series. It becomes trickier to do stats on this sample:" ] }, { "cell_type": "code", "execution_count": null, "id": "131", "metadata": {}, "outputs": [], "source": [ "maxt.plot.hist(edgecolor='k', bins=np.arange(26, 42, 1));" ] }, { "cell_type": "markdown", "id": "132", "metadata": {}, "source": [ "The shape of the ditribution is far from Gaussian. This is a problem we will tackle at another time.\n", "\n", "What we are interested here is to compute the probability of a given temperature to occur each year, based on the data only. Intuitively, the probability of annual maximum temperature exceeding 40°C can be computed easily:" ] }, { "cell_type": "code", "execution_count": null, "id": "133", "metadata": {}, "outputs": [], "source": [ "prob_over_40 = 1 / len(maxt)\n", "f'The probability of this event is {prob_over_40 * 100:.1f}%'" ] }, { "cell_type": "markdown", "id": "134", "metadata": {}, "source": [ "This is because this event occurred only once in the 52 years of the observational record. In other terms, the empirical return period of this event if 52 years. This is easier to understand than 1.9%, yet it's exactly the same concept. Mathematically:\n", "\n", "$R = 1 / p$\n", "\n", "or\n", "\n", "\n", "$p = 1 / R$\n", "\n", "With R the return period and p the probability of the event.\n", "\n", "Now, let's try to work more systematically on the classification of extreme events. Let's start by noticing that the *year* at which an event occurred is irrelevant. What matters is how \"extreme\" an event it. Therefore, we will work with another organisation of the data. We create a new table, where we will organise the events by their \"rank\", from most extreme to least extreme:" ] }, { "cell_type": "code", "execution_count": null, "id": "135", "metadata": {}, "outputs": [], "source": [ "# Create an empty dataframe indexed by the number of years\n", "edf = pd.DataFrame(index=np.arange(1, len(maxt)+1))\n", "\n", "# Add the sorted temperature data\n", "edf['sorted'] = maxt.sort_values(ascending=False).values # Careful! Note the use of .values here" ] }, { "cell_type": "markdown", "id": "136", "metadata": {}, "source": [ "**E: Explore the edf dataframe. What is the maximum temperature? What is the smallest maximum temperature?**\n", "\n", "**E: note that we used `.values` here to assign the values from one dataframe to the other. That's because pandas wants to align indices where possible. Check for example what happens if you remove `.values`. The role of `.values` is to use the numpy arrays (the numbers) to fill the `edf` dataframe, not the combination of index + values.**" ] }, { "cell_type": "code", "execution_count": null, "id": "137", "metadata": {}, "outputs": [], "source": [ "# Your answer here" ] }, { "cell_type": "markdown", "id": "138", "metadata": {}, "source": [ "Good. Now we can notice that the index of the Dataframe represents the \"rank\" of the event, from most extreme to least extreme. Let's verify that events of equal intensity are not given similar ranks. This occurs if we rely on categorical data:" ] }, { "cell_type": "code", "execution_count": null, "id": "139", "metadata": {}, "outputs": [], "source": [ "len(edf), len(np.unique(edf['sorted']))" ] }, { "cell_type": "markdown", "id": "140", "metadata": {}, "source": [ "**E: explore `edf` again. Can you spot the identical values? Can you figure out how this would happen in the data records?**\n", "\n", "In nature, extreme events are truly unique, but our observational systems and the storage of numbers with computers leads to this effect (discretisation).\n", "\n", "OK well, then we have to rank the data differently, since we dont want to assign the same return level with a different rank. Let's use scipy instead (rankdata):" ] }, { "cell_type": "code", "execution_count": null, "id": "141", "metadata": {}, "outputs": [], "source": [ "# rank via scipy instead to deal with duplicate values\n", "edf['ranks'] = np.sort(stats.rankdata(-edf['sorted'])) # Note the \"-\" to rank from max to min" ] }, { "cell_type": "markdown", "id": "142", "metadata": {}, "source": [ "**E: explore edf again. Can you detect what the ranking function did to the data?**\n", "\n", "Let's make sure we have the same number of unique ranks as unique values:" ] }, { "cell_type": "code", "execution_count": null, "id": "143", "metadata": {}, "outputs": [], "source": [ "len(np.unique(edf['sorted'])), len(np.unique(edf['ranks']))" ] }, { "cell_type": "markdown", "id": "144", "metadata": {}, "source": [ "OK all good now. " ] }, { "cell_type": "markdown", "id": "145", "metadata": {}, "source": [ "Next, we compute the empirical probability of exceedance by dividing the rank (r) by the total number of values (n) plus 1:\n", "\n", "**Important! This formula is slightly different from the one used earlier — we divide by n+1 instead of n.**\n", "\n", "Why? Without this adjustment, the probability of exceeding the lowest maximum value (19.05°C) would be 1, meaning the probability of observing a maximum temperature below 19°C in any given year would be zero. This is not realistic.\n", "\n", "Adding n+1 accounts for the fact that our dataset is just a sample, and we haven’t observed all possible values in the true population. See also [Weibull plotting position](https://glossary.ametsoc.org/wiki/Weibull_plotting_position)." ] }, { "cell_type": "code", "execution_count": null, "id": "146", "metadata": {}, "outputs": [], "source": [ "# find exceedance probability\n", "n = len(edf)\n", "edf[\"exceedance\"] = edf[\"ranks\"] / (n + 1)" ] }, { "cell_type": "markdown", "id": "147", "metadata": {}, "source": [ "**E: as an excercise, compute the probability of exceedance by dividing by n instead of n+1, and verify the claim aboce. Does it make sense to use n+1 now?**" ] }, { "cell_type": "code", "execution_count": null, "id": "148", "metadata": {}, "outputs": [], "source": [ "# Your answer here" ] }, { "cell_type": "markdown", "id": "149", "metadata": {}, "source": [ "Despite this small difference, we can verify that the probability of exceedance of the >40°C event matches our simple calculations from above. Now however, we can compute the probabilities of exceedance (and therefore the return values) of all extreme events in the timeseries. Let's go:" ] }, { "cell_type": "code", "execution_count": null, "id": "150", "metadata": {}, "outputs": [], "source": [ "# compute the return period\n", "edf[\"period\"] = 1 / edf[\"exceedance\"]" ] }, { "cell_type": "markdown", "id": "151", "metadata": {}, "source": [ "**E: explore edf again and get aquainted with the new data.**" ] }, { "cell_type": "code", "execution_count": null, "id": "152", "metadata": {}, "outputs": [], "source": [ "# Your answer here" ] }, { "cell_type": "markdown", "id": "153", "metadata": {}, "source": [ "Let's plot (finally!) our long awaited return level plot:" ] }, { "cell_type": "code", "execution_count": null, "id": "154", "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots()\n", "ax.plot(edf['period'], edf['sorted'], \"o\")\n", "ax.grid(True)\n", "ax.set_title('Return periods/levels of annual maximum temperature at Heathrow')\n", "ax.set_xlabel(\"Return Period (years)\")\n", "ax.set_ylabel(\"Return Level (°C)\")\n", "ax.set_xscale(\"linear\") # notice the xscale" ] }, { "cell_type": "markdown", "id": "155", "metadata": {}, "source": [ "For readability and ease of use, we often discuss return periods using a logarithmic scale, referring to events with return periods of 1 year, 10 years, and 100 years. Let’s modify the plot above by changing the x-axis to a logarithmic scale." ] }, { "cell_type": "code", "execution_count": null, "id": "156", "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots()\n", "ax.plot(edf['period'], edf['sorted'], \"o\")\n", "ax.grid(True)\n", "ax.set_title('Return periods/levels of annual maximum temperature at Heathrow (log scale)')\n", "ax.set_xlabel(\"Return Period (years)\")\n", "ax.set_ylabel(\"Return Level (°C)\")\n", "ax.set_xscale(\"log\") # change the xscale" ] }, { "cell_type": "markdown", "id": "157", "metadata": {}, "source": [ "**What is the approximative empirical return period of a 34°C max temperature? And of a 38°C max temperature?**" ] }, { "cell_type": "code", "execution_count": null, "id": "158", "metadata": {}, "outputs": [], "source": [ "# Your answer here" ] }, { "cell_type": "markdown", "id": "159", "metadata": {}, "source": [ "**E: now repeat the steps above but with precipitation instead. Compute the return period of a 40 mm/day event.**" ] }, { "cell_type": "code", "execution_count": null, "id": "160", "metadata": {}, "outputs": [], "source": [ "# Your answer here" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.12.8" } }, "nbformat": 4, "nbformat_minor": 5 }