{ "cells": [ { "cell_type": "markdown", "id": "8da35269", "metadata": {}, "source": [ "# 08: Working with tabular data in Pandas\n", "\n", "" ] }, { "cell_type": "markdown", "id": "adea9fbe", "metadata": {}, "source": [ "## Not that type of Panda -- Python's Pandas package\n", "\n", "Pandas is a powerful, flexible and easy to use open source data analysis and manipulation tool. Pandas is commonly used for operations that would normally be done in a spreadsheet environment and includes powerful data analysis and manipulation tools.\n", "\n", "Let's begin by importing the libraries and setting our data path " ] }, { "cell_type": "code", "execution_count": null, "id": "3b885ca2", "metadata": {}, "outputs": [], "source": [ "from pathlib import Path\n", "import numpy as np\n", "import matplotlib as mpl\n", "import matplotlib.pyplot as plt\n", "import pandas as pd\n", "import statsmodels.api as sm\n", "from dataretrieval import nwis\n", "from matplotlib.backends.backend_pdf import PdfPages\n", "from scipy.signal import detrend\n", "\n", "data_path = Path(\"..\", \"data\", \"pandas\")" ] }, { "cell_type": "markdown", "id": "f32f58cb", "metadata": {}, "source": [ "## Get site information data for part of the Russian River near Guerneville, CA from NWIS " ] }, { "cell_type": "code", "execution_count": null, "id": "6d6d2451-5c5d-43ea-81d4-7daf3e757ad5", "metadata": {}, "outputs": [], "source": [ "bbox = (-123.10, 38.45, -122.90, 38.55)\n", "info, metadata = nwis.get_info(bBox=[str(i) for i in bbox])\n", "if 'geometry' in list(info):\n", " # dataretreival now returns as geopandas dataframe, drop \"geometry\" column for this example\n", " info.drop(columns=[\"geometry\"], inplace=True)" ] }, { "cell_type": "code", "execution_count": null, "id": "97e48826-2eff-4cbb-ac69-f93f2011c0ec", "metadata": {}, "outputs": [], "source": [ "info.to_csv(data_path / \"site_info.csv\", index=False)\n", "info.head()" ] }, { "cell_type": "markdown", "id": "0bc3b4d3", "metadata": {}, "source": [ "Wow that's a lot of gages in this location, let's info on all of them and then use this data to learn about pandas" ] }, { "cell_type": "markdown", "id": "3740b558", "metadata": {}, "source": [ "The data returned to us from our NWIS query `info` was returned to us as pandas `DataFrame`. We'll be working with this to start learning about the basics of pandas." ] }, { "cell_type": "markdown", "id": "ab43af55", "metadata": {}, "source": [ "## Viewing data in pandas\n", "\n", "Pandas has built in methods to inspect `DataFrame` objects. We'll look at a few handy methods:\n", "\n", " - `.head()`: inspect the first few rows of data\n", " - `.tail()`: inspect the last few rows of data\n", " - `.index()`: show the row indexes\n", " - `.columns()`: show the column names\n", " - `.describe()`: statistically describe the data" ] }, { "cell_type": "code", "execution_count": null, "id": "12384282", "metadata": {}, "outputs": [], "source": [ "info.head()" ] }, { "cell_type": "code", "execution_count": null, "id": "fd42e690", "metadata": {}, "outputs": [], "source": [ "info.tail()" ] }, { "cell_type": "code", "execution_count": null, "id": "40ec8e5e", "metadata": {}, "outputs": [], "source": [ "# print the index names\n", "print(info.index)" ] }, { "cell_type": "code", "execution_count": null, "id": "7a1c0d8a", "metadata": {}, "outputs": [], "source": [ "# print the column names\n", "print(info.columns)\n", "\n", "# print the column names as a list\n", "print(list(info))" ] }, { "cell_type": "markdown", "id": "2bff1bcf", "metadata": {}, "source": [ "The `describe()` method is only useful for numerical data. This example is does not have a lot of useful data for `describe()`, however we'll use it again later" ] }, { "cell_type": "code", "execution_count": null, "id": "f86033a1", "metadata": {}, "outputs": [], "source": [ "info.describe()" ] }, { "cell_type": "markdown", "id": "09d73a52", "metadata": {}, "source": [ "## Getting data from a pandas dataframe\n", "\n", "There are multiple methods to get data out of a pandas dataframe as either a \"series\", numpy array, or a list\n", "\n", "Let's start by getting data as a series using a few methods" ] }, { "cell_type": "code", "execution_count": null, "id": "1e4d13e4", "metadata": {}, "outputs": [], "source": [ "# get a series of site numbers by key\n", "info[\"site_no\"]" ] }, { "cell_type": "code", "execution_count": null, "id": "daff4fea", "metadata": {}, "outputs": [], "source": [ "# get station names by attribute\n", "info.station_nm" ] }, { "cell_type": "markdown", "id": "5bbba79a", "metadata": {}, "source": [ "getting data from a dataframe as a numpy array can be accomplished by using `.values`" ] }, { "cell_type": "code", "execution_count": null, "id": "d98d409c", "metadata": {}, "outputs": [], "source": [ "info.station_nm.values" ] }, { "cell_type": "markdown", "id": "31def645", "metadata": {}, "source": [ "### Selection by position\n", "\n", "We can get data by position in the dataframe using the `.iloc` attribute " ] }, { "cell_type": "code", "execution_count": null, "id": "730501c2", "metadata": {}, "outputs": [], "source": [ "info.iloc[0:2, 1:4]" ] }, { "cell_type": "markdown", "id": "b4d7c291", "metadata": {}, "source": [ "### Selection by label\n", "\n", "Pandas allows the user to get data from the dataframe by index and column labels" ] }, { "cell_type": "code", "execution_count": null, "id": "b430c660", "metadata": {}, "outputs": [], "source": [ "info.loc[0:3, [\"site_no\", \"station_nm\", \"site_tp_cd\"]]" ] }, { "cell_type": "markdown", "id": "43700bbe", "metadata": {}, "source": [ "### Boolean indexing\n", "\n", "pandas dataframes supports boolean indexing that allows a user to create a new dataframe with only the data that meets a boolean condition defined by the user.\n", "\n", "Let's get a dataframe of only groundwater sites from the `info` dataframe" ] }, { "cell_type": "code", "execution_count": null, "id": "24ea2d22", "metadata": {}, "outputs": [], "source": [ "dfgw = info[info[\"site_tp_cd\"] == \"GW\"]\n", "dfgw" ] }, { "cell_type": "markdown", "id": "b55d5649", "metadata": {}, "source": [ "### Reading and writing data to `.csv` files\n", "\n", "Pandas has support to both read and write many types of files. For this example we are focusing on `.csv` files. For information on other file types that are supported see the [ten minutes to pandas](https://pandas.pydata.org/docs/user_guide/10min.html#importing-and-exporting-data) tutorial documentation" ] }, { "cell_type": "markdown", "id": "8d69feb7", "metadata": {}, "source": [ "For this part we'll write a new `.csv` file of the groundwater sites that we found in NWIS using `to_csv()`.\n", "\n", "`to_csv()` has a bunch of handy options for writing to file. For this example, I'm going to drop the index column while writing by passing `index=False`" ] }, { "cell_type": "code", "execution_count": null, "id": "ecb2281a", "metadata": {}, "outputs": [], "source": [ "csv_file = data_path / \"RussianRiverGWsites.csv\"\n", "dfgw.to_csv(csv_file, index=False)" ] }, { "cell_type": "markdown", "id": "485540c8", "metadata": {}, "source": [ "Now we can load the csv file back into a pandas dataframe with the `read_csv()` method." ] }, { "cell_type": "code", "execution_count": null, "id": "da14892b", "metadata": {}, "outputs": [], "source": [ "dfgw2 = pd.read_csv(csv_file)\n", "dfgw2" ] }, { "cell_type": "markdown", "id": "5e5db256", "metadata": {}, "source": [ "## Class exercise 1 \n", "\n", "Using the methods presented in this notebook create a `DataFrame` of surface water sites from the `info` dataframe, write it to a csv file named `\"RussianRiverSWsites.csv\"`, and read it back in as a new `DataFrame`" ] }, { "cell_type": "code", "execution_count": null, "id": "126bba4e", "metadata": {}, "outputs": [], "source": [ "csv_name = data_path / \"RussianRiverSWsites.csv\"\n", "dfsw = info[info.site_tp_cd == \"ST\"]\n", "dfsw.to_csv(csv_name, index=False)\n", "\n", "dfsw = pd.read_csv(csv_name)\n", "dfsw.head()" ] }, { "cell_type": "markdown", "id": "2182d335", "metadata": {}, "source": [ "## Get timeseries data from NWIS for surface water sites\n", "\n", "Now to start working with stream gage data from NWIS. We're going to send all surface water sites to NWIS and get only sites with daily discharge measurements.\n", "\n", "([11467002](https://waterdata.usgs.gov/nwis/inventory?agency_code=USGS&site_no=11467002))" ] }, { "cell_type": "code", "execution_count": null, "id": "edd1afbc", "metadata": {}, "outputs": [], "source": [ "# filter for only SW sites\n", "dfsw = info[info.site_tp_cd == \"ST\"]\n", "\n", "# create query\n", "sites = dfsw.site_no.to_list()\n", "pcode = \"00060\"\n", "start_date = \"1939-09-22\"\n", "end_date = \"2023-12-31\"\n", "df, metadata = nwis.get_dv(\n", " sites=sites, \n", " start=start_date, \n", " end=end_date, \n", " parameterCd=pcode, \n", " multi_index=False\n", ")\n", "df.head()" ] }, { "cell_type": "markdown", "id": "9f00ba06", "metadata": {}, "source": [ "Awesome! There are two gages in this study area that have daily discharge data!\n", "\n", "### Inspect the data using `.plot()`" ] }, { "cell_type": "code", "execution_count": null, "id": "e2970f2e", "metadata": {}, "outputs": [], "source": [ "ax = df.plot()\n", "ax.set_ylabel(r\"cubic meter per second\")\n", "ax.set_yscale(\"log\")" ] }, { "cell_type": "markdown", "id": "fa2c35dd", "metadata": {}, "source": [ "### Updating column names\n", "\n", "Let's update column names for these gages based on their locations instead of their USGS gage id. First we'll get name information from the `info` class and then we'll use the station name information to remap the column names.\n", "\n", "the `rename()` method accepts a dictionary that is formatted `{current_col_name: new_col_name}`" ] }, { "cell_type": "code", "execution_count": null, "id": "b800fbbc", "metadata": {}, "outputs": [], "source": [ "df" ] }, { "cell_type": "code", "execution_count": null, "id": "88b7dd45", "metadata": {}, "outputs": [], "source": [ "df2 = df[['site_no', '00060_Mean']].pivot(columns=['site_no'])\n", "df2.columns = df2.columns.get_level_values(1)\n", "\n", "site_names = {site_no: '_'.join(name.split()[:2]).lower() for site_no, name \n", " in info.groupby('site_no')['station_nm'].first().loc[df2.columns].items()}\n", "df2.rename(columns=site_names, inplace=True)\n", "df2.plot()" ] }, { "cell_type": "markdown", "id": "f2c20720", "metadata": {}, "source": [ "### Adding a new column of data to the dataframe\n", "\n", "Adding new columns of data is relatively simple in pandas. The call signature is similar to a dictionary where df[\"key\"] = \"some_array_of_values\"\n", "\n", "Let's add a year column (water year) to the dataframe by getting the year from the index and adjusting it using `np.where()`" ] }, { "cell_type": "code", "execution_count": null, "id": "39ebd1e7", "metadata": {}, "outputs": [], "source": [ "df2['year'] = np.where(df2.index.month >=10, df2.index.year, df2.index.year - 1)" ] }, { "cell_type": "markdown", "id": "5ccdd9e2", "metadata": {}, "source": [ "### Manipulating data\n", "\n", "Columns in a pandas dataframe can be manipulated similarly to working with a dictionary of numpy arrays. The discharge data units are $\\frac{ft^3}{s}$, let's accumulate this to daily discharge in $\\frac {ft^3}{d}$.\n", "\n", "The conversion for this is:\n", "\n", "## $ \\frac{1 ft^3}{s} \\times \\frac{60s}{min} \\times \\frac{60min}{hour} \\times \\frac{24hours}{day} \\rightarrow \\frac{ft^3}{day} $\n", "\n", "and then convert that to acre-feet (1 acre-ft = 43559.9)" ] }, { "cell_type": "code", "execution_count": null, "id": "b294a0ac", "metadata": {}, "outputs": [], "source": [ "df2" ] }, { "cell_type": "code", "execution_count": null, "id": "81980161", "metadata": {}, "outputs": [], "source": [ "conv = 60 * 60 * 24\n", "df2[\"russian_r_cfd\"] = df2.russian_r * conv\n", "df2[\"russian_r_af\"] = df2.russian_r_cfd / 43559.9\n", "df2\n", "\n", "# now let's plot up discharge from 2018 - 2020\n", "dfplt = df2[(df2[\"year\"] >= 2015) & (df2[\"year\"] <= 2020)]\n", "ax = dfplt[\"russian_r_af\"].plot()\n", "ax.set_yscale(\"log\")\n" ] }, { "cell_type": "markdown", "id": "112292ee", "metadata": {}, "source": [ "## Class exercise 2:\n", "\n", "Using the methods described in the notebook. Convert the discharge for `\"austin_c\"` from $\\frac{m^3}{s}$ acre-ft by adding two additional fields to `df2` named `\"austin_c_cfs\"` and `\"austin_c_af\"`. \n", "\n", "After these two fields have been added to df2 use boolean indexing to create a new dataframe from 2015 through 2020. Finally plot the discharge data (in acre-ft) for \"austin_cr\".\n", "\n", "**bonus exercise** try to plot both `\"russian_r_af\"` and `\"austin_c_af\"` on the same plot" ] }, { "cell_type": "code", "execution_count": null, "id": "e92c3bc7", "metadata": {}, "outputs": [], "source": [ "conv = 60 * 60 * 24\n", "df2[\"austin_c_cfd\"] = df2.austin_c * conv\n", "df2[\"austin_c_af\"] = df2.austin_c_cfd / 43559.9\n", "df2\n", "\n", "# now let's plot up discharge from 2018 - 2020\n", "dfplt = df2[(df2[\"year\"] >= 2015) & (df2[\"year\"] <= 2020)]\n", "ax = dfplt[\"austin_c_af\"].plot()\n", "ax.set_yscale(\"log\")\n" ] }, { "cell_type": "code", "execution_count": null, "id": "6f349f19", "metadata": {}, "outputs": [], "source": [ "# now let's plot up discharge from 2018 - 2020\n", "dfplt = df2[(df2[\"year\"] >= 2015) & (df2[\"year\"] <= 2020)]\n", "ax = dfplt[[\"russian_r_af\",\"austin_c_af\"]].plot()\n", "ax.set_yscale(\"log\")" ] }, { "cell_type": "markdown", "id": "28e18569", "metadata": {}, "source": [ "### groupby: grouping data and performing mathematical operations on it\n", "\n", "`groupby()` is a powerful method that allows for performing statistical operations over a groups of \"common\" data within a dataframe. \n", "\n", "For this example we'll use it to get mean daily flows for the watershed. Pandas will group all common days of the year with each other and then calculate the mean value of these via the function `.mean()`. `groupby()` also supports other operations such as `.median()`, `.sum()`, `max()`, `min()`, `std()`, and other functions. " ] }, { "cell_type": "code", "execution_count": null, "id": "61938c67", "metadata": {}, "outputs": [], "source": [ "df2[\"day_of_year\"] = df2.index.day_of_year\n", "\n", "df_mean_day = df2.groupby(by=[\"day_of_year\"], as_index=False)[[\"russian_r_af\", \"austin_c_af\"]].mean()\n", "ax = df_mean_day[[\"russian_r_af\", \"austin_c_af\"]].plot()\n", "ax.set_yscale(\"log\")" ] }, { "cell_type": "markdown", "id": "c0b36148", "metadata": {}, "source": [ "We can see that around March flow starts decreasing heavily and doesn't recover until sometime in october. What's going on here? Let's load some climate data and see what's happening." ] }, { "cell_type": "code", "execution_count": null, "id": "bcc86071", "metadata": {}, "outputs": [], "source": [ "cimis_file = data_path / \"santa_rosa_CIMIS_83.csv\"\n", "df_cimis = pd.read_csv(cimis_file)\n", "drop_list = [\"qc\"] + [f\"qc.{i}\" for i in range(1, 22)]\n", "df_cimis.drop(columns=drop_list, inplace=True)\n", "df_cimis" ] }, { "cell_type": "code", "execution_count": null, "id": "072f50f4", "metadata": {}, "outputs": [], "source": [ "df_mean_clim = df_cimis.groupby(by=[\"Jul\"], as_index=False)[[\"ETo (in)\", \"Precip (in)\"]].mean()\n", "df_mean_clim = df_mean_clim.rename(columns={\"Jul\": \"day_of_year\"})\n", "df_mean_clim.head()" ] }, { "cell_type": "markdown", "id": "bb5edf3c", "metadata": {}, "source": [ "### Joining two `DataFrames` with an inner join\n", "\n", "Inner joins can be made in pandas with the `merge()` method. As inner join, joins two dataframes on a common key or values, when a key or value exists in one dataframe, but not the other that row of data will be excluded from the joined dataframe. \n", "\n", "There are a number of other ways to join dataframes in pandas. Detailed examples and discussion of the different merging methods can be found [here](https://realpython.com/pandas-merge-join-and-concat/)" ] }, { "cell_type": "code", "execution_count": null, "id": "d2b055b6", "metadata": {}, "outputs": [], "source": [ "df_merged = pd.merge(df_mean_day, df_mean_clim, on=[\"day_of_year\",])\n", "df_merged.head()" ] }, { "cell_type": "code", "execution_count": null, "id": "1915d1b9", "metadata": {}, "outputs": [], "source": [ "cols = [i for i in list(df_merged) if i != \"day_of_year\"]\n", "ax = df_merged[cols].plot()\n", "ax.set_yscale(\"log\")" ] }, { "cell_type": "markdown", "id": "3ed8ec21", "metadata": {}, "source": [ "This starts to make a little more sense! We can see that this area gets most of it's precipitation in the winter and early spring. The evapotranspiration curve shows what we'd expect: that there is more evapotranspiration in the summer than the winter. " ] }, { "cell_type": "markdown", "id": "47d9ca7e-fafe-4280-a2a1-7a87fb1ef3b5", "metadata": {}, "source": [ "## Stepping back for a minute: exploring Merge and Concat\n", "\n", "Merge and concat are very powerful and important methods for combining pandas DataFrames (and later geopandas GeoDataFrames). This section takes a moment, as an aside to explore some of the common funtions and options for performing merges and concatenation. " ] }, { "cell_type": "markdown", "id": "893e3dab-72c4-4702-a50b-d9ee832af85f", "metadata": {}, "source": [ "### The `merge` method\n", "\n", "Pandas `merge` method allows the user to join two dataframes based on column values. There are a few methods that can be used to control the type of merge and the result.\n", "\n", " - `inner`: operates as a intersection of the data and keeps rows from both dataframes if the row has a matching key in both dataframes.\n", " - `outer`: operates as a union and joins all rows from both dataframes, fills columns with nan where data is not available for a row.\n", " - `left`: keeps all rows from the 'left' dataframe and joins columns from the right\n", " - `right`: keeps all rows from the 'right' dataframe and joins columns from the left" ] }, { "cell_type": "code", "execution_count": null, "id": "1472f91c-6cd7-429a-8203-9ba2648c5aae", "metadata": {}, "outputs": [], "source": [ "ex0 = {\n", " \"name\": [], \n", " \"childhood_dream_job\": []\n", "}\n", "ex1 = {\n", " \"name\": [],\n", " \"age\": []\n", "}" ] }, { "cell_type": "code", "execution_count": null, "id": "e873668d-b867-498e-bf4a-66c465bc90b3", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "id": "4ac4b420-0be1-400c-b725-4d914a761889", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "id": "97d66c72-fe0a-42d0-a6cd-a8613db2da75", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "4c998ffd-f62b-4b0a-8f16-55f814ea15b7", "metadata": {}, "source": [ "### The `concat` method\n", "\n", "Pandas `concat` method allows the user to concatenate dataframes together" ] }, { "cell_type": "code", "execution_count": null, "id": "aa7f3542-38eb-4082-a6f3-e8baa5bdac29", "metadata": {}, "outputs": [], "source": [ "gw = info[info.site_tp_cd == \"GW\"]\n", "sw = info[info.site_tp_cd == \"ST\"]\n", "sw = sw[list(sw)[0:8]]\n", "sw.head()" ] }, { "cell_type": "code", "execution_count": null, "id": "fc81aead-4d2c-48b0-a868-46c546c0901c", "metadata": {}, "outputs": [], "source": [ "newdf = pd.concat((gw, sw), ignore_index=True)\n", "newdf" ] }, { "cell_type": "markdown", "id": "835dd7c5-bce8-434b-889e-5ef837fac152", "metadata": {}, "source": [ "## Back to working with our climate and stream data" ] }, { "cell_type": "markdown", "id": "0b486148", "metadata": {}, "source": [ "Is there anything else we can look at that might give us more insights into this basin? \n", "\n", "Let's try looking at long term, yearly trends to see what's there." ] }, { "cell_type": "code", "execution_count": null, "id": "e62109ef", "metadata": {}, "outputs": [], "source": [ "year = []\n", "for dstr in df_cimis.Date.values:\n", " mo, d, yr = [int(i) for i in dstr.split(\"/\")]\n", " if mo < 10:\n", " year.append(yr - 1)\n", " else:\n", " year.append(yr)\n", "df_cimis[\"year\"] = year\n", "df_cimis.head()" ] }, { "cell_type": "code", "execution_count": null, "id": "826df4a9", "metadata": {}, "outputs": [], "source": [ "df2.head()" ] }, { "cell_type": "markdown", "id": "49a97ec7", "metadata": {}, "source": [ "### `.aggregate()`\n", "\n", "The pandas `.aggregate()` can be used with `.groupy()` to perform multiple statistical operations.\n", "\n", "Example usage could be:\n", "\n", "```python\n", " \n", " agdf = df2.groupby(by=[\"day_of_year\"], as_index=False)[\"austin_c_af\"].aggregate([np.min, np.max, np.mean, np.std])\n", "```\n", "\n", "And would return a dataframe grouped by the day of year and take the min, max, mean, and standard deviation of these data." ] }, { "cell_type": "markdown", "id": "ffd30f76", "metadata": {}, "source": [ "## Class exercise 3: groupy, aggregate, and merge\n", "\n", "For this exercise we need to produce yearly mean and standard deviation flows for the `russian_r_af` variable and merge the data with yearly mean precipitation and evapotranspiration \n", "\n", "**Hints:** \n", " - `df2` and `df_cimis` are the dataframes these operations should be performed on \n", " - use `.groupby()` to group by the year and do `precip` and `ETo` in separate groupby routines\n", " - rename the columns in the ETo aggregated dataframe and the Precip dataframe to `mean_et`, `std_et`, `mean_precip`, `std_precip`\n", " - `np.mean` and `np.std` can be used to calculate mean and std flows \n", " \n", "Name your final joined dataframe `df_yearly`" ] }, { "cell_type": "code", "execution_count": null, "id": "3e0ac248", "metadata": {}, "outputs": [], "source": [ "dfg0 = df2.groupby(by=[\"year\"], as_index=False)[\"russian_r_af\"].aggregate(['mean', 'std'])\n", "dfg1 = df_cimis.groupby(by=[\"year\"], as_index=False)[\"ETo (in)\"].aggregate(['mean', 'std'])\n", "dfg2 = df_cimis.groupby(by=[\"year\"], as_index=False)[\"Precip (in)\"].aggregate(['mean', 'std'])\n", "\n", "dfm1 = pd.merge(dfg0, dfg1, on=[\"year\"], suffixes=(None, \"_et\"))\n", "df_yearly = pd.merge(dfm1, dfg2, on=[\"year\"], suffixes=(None, \"_precip\"))\n", "\n", "df_yearly" ] }, { "cell_type": "markdown", "id": "c1f8ccef", "metadata": {}, "source": [ "Lets examine the long term flow record for the Russian River" ] }, { "cell_type": "code", "execution_count": null, "id": "31bcbf39", "metadata": {}, "outputs": [], "source": [ "df_y_flow = df2.groupby(by=[\"year\"], as_index=False)[\"russian_r_af\"].aggregate(['mean', 'std'])\n", "\n", "fig = plt.figure(figsize=(12,4))\n", "\n", "lower_ci = df_y_flow[\"mean\"] - 2 * df_y_flow['std']\n", "lower_ci = np.where(lower_ci < 0, 0, lower_ci)\n", "upper_ci = df_y_flow[\"mean\"] + 2 * df_y_flow['std']\n", "ax = df_y_flow[\"mean\"].plot(style=\"b.-\")\n", "ax.fill_between(df_y_flow.index, lower_ci, upper_ci, color=\"b\", alpha=0.5);" ] }, { "cell_type": "markdown", "id": "0f91c518", "metadata": {}, "source": [ "From this plot it looks like the long term trend in yearly discharge in the Russian River has been decreasing. We will revisit and test this later. \n", "\n", "First let's see if there are any relationships between yearly discharge and climate" ] }, { "cell_type": "code", "execution_count": null, "id": "06406208", "metadata": {}, "outputs": [], "source": [ "dfg0 = df2.groupby(by=[\"year\"], as_index=False)[\"russian_r_af\"].aggregate(['mean', 'std', 'sum'])\n", "dfg1 = df_cimis.groupby(by=[\"year\"], as_index=False)[\"ETo (in)\"].aggregate(['mean', 'std', 'sum'])\n", "dfg2 = df_cimis.groupby(by=[\"year\"], as_index=False)[\"Precip (in)\"].aggregate(['mean', 'std', 'sum'])\n", "dfm1 = pd.merge(dfg0, dfg1, on=[\"year\"], suffixes=(None, \"_et\"))\n", "df_yearly = pd.merge(dfm1, dfg2, on=[\"year\"], suffixes=(None, \"_precip\"))\n", "\n", "fig = plt.figure(figsize=(12,4))\n", "\n", "lower_ci = df_yearly[\"mean\"] - 2 * df_yearly['std']\n", "lower_ci = np.where(lower_ci < 0, 0, lower_ci)\n", "upper_ci = df_yearly[\"mean\"] + 2 * df_yearly['std']\n", "ax = df_yearly[\"mean\"].plot(style=\"b.-\", label=\"flow_af\")\n", "ax.fill_between(df_yearly.index, lower_ci, upper_ci, color=\"b\", alpha=0.5)\n", "ax2 = ax.twinx()\n", "ax2.plot(df_yearly.index, df_yearly.sum_precip, \"k--\", lw=1.5, label=\"Precip\")\n", "ax2.plot(df_yearly.index, df_yearly.sum_et, \"r.--\", lw=1.5, label=\"ET\")\n", "plt.legend()\n", "plt.show();" ] }, { "cell_type": "markdown", "id": "7a0f73bb", "metadata": {}, "source": [ "As expected, precipitation amount is the main driver of the yearly discharge regime here" ] }, { "cell_type": "markdown", "id": "5834f46d", "metadata": {}, "source": [ "### Baseflow separation\n", "\n", "Baseflow separation is a method to separate the quick response hydrograph (storm runoff) from the long term flow. We're going to go back to the complete daily dataset to perform this operation. The following cell contains a low pass filtration method that is used for digital baseflow separation. We'll use these in our analysis" ] }, { "cell_type": "code", "execution_count": null, "id": "21356a99", "metadata": {}, "outputs": [], "source": [ "def _baseflow_low_pass_filter(arr, beta, enforce):\n", " \"\"\"\n", " Private method to apply digital baseflow separation filter\n", " (Lyne & Hollick, 1979; Nathan & Mcmahon, 1990;\n", " Boughton, 1993; Chapman & Maxwell, 1996).\n", "\n", " This method should not be called by the user!\n", "\n", " Method removes \"spikes\" which would be consistent with storm flow\n", " events from data.\n", "\n", " Parameters\n", " ----------\n", " arr : np.array\n", " streamflow or municipal pumping time series\n", "\n", " beta : float\n", " baseflow filtering parameter that ranges from 0 - 1\n", " values in 0.8 - 0.95 range used in literature for\n", " streamflow baseflow separation\n", "\n", " enforce : bool\n", " enforce physical constraint of baseflow less than measured flow\n", "\n", " Returns\n", " -------\n", " np.array of filtered data\n", " \"\"\"\n", " # prepend 10 records to data for initial spin up\n", " # these records will be dropped before returning data to user\n", " qt = np.zeros((arr.size + 10,), dtype=float)\n", " qt[0:10] = arr[0:10]\n", " qt[10:] = arr[:]\n", "\n", " qdt = np.zeros((arr.size + 10,), dtype=float)\n", " qbf = np.zeros((arr.size + 10,), dtype=float)\n", "\n", " y = (1. + beta) / 2.\n", "\n", " for ix in range(qdt.size):\n", " if ix == 0:\n", " qbf[ix] = qt[ix]\n", " continue\n", "\n", " x = beta * qdt[ix - 1]\n", " z = qt[ix] - qt[ix - 1]\n", " qdt[ix] = x + (y * z)\n", "\n", " qb = qt[ix] - qdt[ix]\n", " if enforce:\n", " if qb > qt[ix]:\n", " qbf[ix] = qt[ix]\n", " else:\n", " qbf[ix] = qb\n", "\n", " else:\n", " qbf[ix] = qb\n", "\n", " return qbf[10:]\n", "\n", "\n", "def baseflow_low_pass_filter(arr, beta=0.9, T=1, enforce=True):\n", " \"\"\"\n", " User method to apply digital baseflow separation filter\n", " (Lyne & Hollick, 1979; Nathan & Mcmahon, 1990;\n", " Boughton, 1993; Chapman & Maxwell, 1996).\n", "\n", " Method removes \"spikes\" which would be consistent with storm flow\n", " events from data.\n", "\n", " Parameters\n", " ----------\n", " arr : np.array\n", " streamflow or municipal pumping time series\n", "\n", " beta : float\n", " baseflow filtering parameter that ranges from 0 - 1\n", " values in 0.8 - 0.95 range used in literature for\n", " streamflow baseflow separation\n", "\n", " T : int\n", " number of recursive filtering passes to apply to the data\n", "\n", " enforce : bool\n", " enforce physical constraint of baseflow less than measured flow\n", "\n", " Returns\n", " -------\n", " np.array of baseflow\n", " \"\"\"\n", " for _ in range(T):\n", " arr = _baseflow_low_pass_filter(arr, beta, enforce)\n", "\n", " return arr\n" ] }, { "cell_type": "markdown", "id": "17636f07", "metadata": {}, "source": [ "## Class exercise 4: baseflow separation\n", "\n", "Using the full dataframe discharge dataframe, `df2`, get the baseflows for the Russian River (in acre-ft) and add these to the dataframe as a new column named `russian_bf_af`.\n", "\n", "The function `baseflow_low_pass_filter()` will be used to perform baseflow. This function takes a numpy array of stream discharge and runs a digital low pass filtration method to calculate baseflow. \n", "\n", "After performing baseflow separation plot both the baseflow and the total discharge for the years 2015 - 2017.\n", "\n", "**Hints:**\n", " - make sure to use the `.values` property when you get discharge data from the pandas dataframe\n", " - feel free to play with the input parameters `beta` and `T`. \n", " - `beta` is commonly between 0.8 - 0.95 for hydrologic problems.\n", " - `T` is the number of filter passes over the data, therefore increasing `T` will create a smoother profile. I recommend starting with a `T` value of around 5 and adjusting it to see how it changes the baseflow profile." ] }, { "cell_type": "code", "execution_count": null, "id": "99dee893", "metadata": {}, "outputs": [], "source": [ "baseflow = baseflow_low_pass_filter(df2[\"russian_r_af\"].values, beta=0.90, T=5)\n", "df2[\"russian_bf_af\"] = baseflow\n", "\n", "fig = plt.figure(figsize=(15, 10))\n", "tdf = df2[(df2.year >= 2015) & (df2.year <= 2017)]\n", "ax = tdf[[\"russian_r_af\", \"russian_bf_af\"]].plot()\n", "ax.set_yscale(\"log\");" ] }, { "cell_type": "markdown", "id": "fc78258f", "metadata": {}, "source": [ "## Looking at trends and seasonal cycles in the data with `statsmodels`\n", "\n", "`statsmodels` is a Python module that provides classes and functions for the estimation of many different statistical models, as well as for conducting statistical tests, and statistical data exploration. For more information on `statsmodels` and it's many features visit [the getting started guide](https://www.statsmodels.org/dev/gettingstarted.html)\n", "\n", "For this example we'll be using time series analysis to perform a seasonal decomposition using moving averages. The `statsmodels` function `seasonal_decompose` will perform this analysis and return the data trend, seasonal fluctuations, and residuals." ] }, { "cell_type": "code", "execution_count": null, "id": "3be15037", "metadata": {}, "outputs": [], "source": [ "df2[\"russian_bf_af\"] = baseflow_low_pass_filter(df2[\"russian_r_af\"].values, beta=0.9, T=5)\n", "# drop nan values prior to running decomposition\n", "dfsm = df2[df2['russian_bf_af'].notna()]\n", "\n", "dfsm.head()\n", "# decompose\n", "decomposition = sm.tsa.seasonal_decompose(dfsm[\"russian_bf_af\"], period=365, model=\"additive\")\n", " \n", "decomposition.plot()\n", "mpl.rcParams['figure.figsize'] = [14, 14];" ] }, { "cell_type": "markdown", "id": "6f80c43e", "metadata": {}, "source": [ "Now let's plot our baseflow trend against the baseflow we calculated " ] }, { "cell_type": "code", "execution_count": null, "id": "5aea97b2", "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(figsize=(8, 8))\n", "\n", "ax = df2[\"russian_r_af\"].plot(ax=ax)\n", "ax.plot(decomposition.trend.index, decomposition.trend, \"r-\", lw=1.5)\n", "ax.set_yscale('log')\n", "plt.savefig(\"russian_river.png\")\n" ] }, { "cell_type": "markdown", "id": "3b253563", "metadata": {}, "source": [ "What's going on at the end of the 1970's? Why is the baseflow trend so low and peak flows are also low? \n", "\n", "From [\"California's Most Significant Droughts\"](https://cawaterlibrary.net/wp-content/uploads/2017/05/CalSignficantDroughts_v10_int.pdf), Chapter 3 \n", "\n", "The setting for the 1976-77 drought differed significantly from the dry times of the 1920s-30s. Although\n", "only a two-year event, its hydrology was severe.\n", "Based on 114 years of computed statewide runoff,\n", "1977 occupies rank 114 (driest year) and 1976 is in\n", "rank 104. The drought was notable for the impacts\n", "experienced by water agencies that were unprepared\n", "for such conditions. One reason for the lack of\n", "preparedness was the perception of relatively ample\n", "water supplies in most areas of the state. " ] }, { "cell_type": "markdown", "id": "ec0d85d6", "metadata": {}, "source": [ "We can also see the analysis that there is an overall reduction in baseflow from 1940's to the 2020's" ] }, { "cell_type": "markdown", "id": "019341c3", "metadata": {}, "source": [ "### Writing yearly plots of data to a PDF document\n", "\n", "We can write the plots we make with pandas to a multipage pdf using matplotlib's `PdfPages`. For this example we'll make yearly plots of all of the discharge, baseflow, and baseflow trend and bundle them into a single pdf document." ] }, { "cell_type": "code", "execution_count": null, "id": "a430d0ac", "metadata": {}, "outputs": [], "source": [ "# add the decomposition trend to the dataframe to make plotting easy\n", "df2[\"baseflow_trend\"] = decomposition.trend\n", "\n", "pdf_file = data_path / \"all_years_russian.pdf\"\n", "with PdfPages(pdf_file) as outpdf:\n", " for year in sorted(df2.year.unique()):\n", " plt.figure()\n", " tdf = df2[df2.year == year]\n", " tdf[[\"russian_r_af\", \"russian_bf_af\", \"baseflow_trend\"]].plot()\n", " plt.xlabel(\"date\")\n", " plt.ylabel(\"acre-feet\")\n", " plt.ylim([1, 300000])\n", " plt.yscale(\"log\")\n", " plt.tight_layout()\n", " outpdf.savefig()\n", " plt.close('all')" ] }, { "cell_type": "markdown", "id": "f2761901", "metadata": {}, "source": [ "Go into the \"data/pandas\" and open up the file \"all_years_russian.pdf\" to check out the plots." ] }, { "cell_type": "markdown", "id": "94ffa931", "metadata": {}, "source": [ "## Back to the basics:\n", "\n", "Now that we've done some analysis on real data, let's go back to the basic and briefly discuss how to adjust data values in pandas\n", "\n", "\n", "### Adjusting a whole column of data values\n", "\n", "Because pandas series objects (columns) are built off of numpy, we can perform mathematical operations in place on a whole column similarly to numpy arrays\n" ] }, { "cell_type": "code", "execution_count": null, "id": "45e51e85", "metadata": {}, "outputs": [], "source": [ "df3 = df2.copy()\n", "# Adding values\n", "df3[\"russian_r\"] += 10000\n", "\n", "# Subtracting values\n", "df3[\"russian_r\"] -= 10000\n", "\n", "# multiplication\n", "df3[\"russian_r\"] *= 5\n", "\n", "# division\n", "df3[\"russian_r\"] /= 5\n", "\n", "# more complex operations\n", "df3[\"russian_r\"] = np.log10(df3[\"russian_r\"].values)\n", "df3.head()" ] }, { "cell_type": "markdown", "id": "ee819e35", "metadata": {}, "source": [ "### updating values by position\n", "\n", "We can update single or multiple values by their position in the dataframe using `.iloc`" ] }, { "cell_type": "code", "execution_count": null, "id": "ba429e93", "metadata": {}, "outputs": [], "source": [ "df3.iloc[0, 1] = 999\n", "df3.head()" ] }, { "cell_type": "markdown", "id": "d37f204f", "metadata": {}, "source": [ "### updating values based on location\n", "\n", "We can update values in the dataframe based on their index and column headers too using `.loc`" ] }, { "cell_type": "code", "execution_count": null, "id": "d2df90e4", "metadata": {}, "outputs": [], "source": [ "df3.loc[df3.index[0], 'russian_r'] *= 100\n", "df3.head()" ] }, { "cell_type": "markdown", "id": "b6ee8610", "metadata": {}, "source": [ "## Class exercise 5: setting data\n", "\n", "The russian river discharge in cubic feet per day highly variable and represented by large numbers. To make this data easier to read and plot, update it to a log10 representation.\n", "\n", "Then use either a location based method to change the value in austin_c on 10-4-1939 to 0.\n", "\n", "Finally display the first 5 records in the notebook." ] }, { "cell_type": "code", "execution_count": null, "id": "ff4b6eda", "metadata": {}, "outputs": [], "source": [ "df3[\"russian_r_cfd\"] = np.log10(df3[\"russian_r_cfd\"])\n", "df3.loc[df3.index[3], \"austin_c\"] = 0\n", "df3.head()" ] }, { "cell_type": "markdown", "id": "f1407f4b", "metadata": {}, "source": [ "## Now for some powerful analysis using fourier analysis to extract the data's signal.\n", "\n", "Some background for your evening viewing: https://youtu.be/spUNpyF58BY\n", "\n", "We have to detrend the data for fast fourier transforms to work properly. Here's a discussion on why: \n", "\n", "https://groups.google.com/forum/#!topic/comp.dsp/kYDZqetr_TU\n", "\n", "Fortunately we can easily do this in python using scipy!\n", "\n", "https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.detrend.html" ] }, { "cell_type": "markdown", "id": "442b8004", "metadata": {}, "source": [ "Let's create a new dataframe with only the data we're interested in analysing." ] }, { "cell_type": "code", "execution_count": null, "id": "8cfa9ee7", "metadata": {}, "outputs": [], "source": [ "df_data = df2[[\"russian_r_af\", \"year\"]]\n", "df_data = df_data.rename(columns={\"russian_r_af\": \"Q\"})\n", "df_data.describe()" ] }, { "cell_type": "markdown", "id": "7a5a2210", "metadata": {}, "source": [ "#### Now we're going to detrend the data\n", "\n", "But first let's also set the water year by date on this dataframe and drop the single nan value from the dataframe" ] }, { "cell_type": "code", "execution_count": null, "id": "1520c21c", "metadata": {}, "outputs": [], "source": [ "df_data['water_year'] = df_data.index.shift(30+31+31,freq='d').year\n", "df_data.dropna(inplace=True)\n", "df_data[\"detrended\"] = detrend(df_data.Q.values)" ] }, { "cell_type": "code", "execution_count": null, "id": "e8104b3e", "metadata": {}, "outputs": [], "source": [ "fig = plt.figure(figsize=(12,6))\n", "ax = plt.subplot(3,1,1)\n", "plt.plot(df_data.detrended)\n", "plt.title('detrended signal')\n", "\n", "plt.subplot(3,1,2)\n", "plt.plot(df_data.index, df_data.Q)\n", "plt.title('Raw Signal')\n", "\n", "plt.subplot(3,1,3)\n", "plt.plot(df_data.index, df_data.Q - df_data.detrended)\n", "plt.title('Difference')\n", "\n", "plt.tight_layout()" ] }, { "cell_type": "markdown", "id": "af16ca60", "metadata": {}, "source": [ "### Evaluate and plot the Period Spectrum to see timing of recurrence\n", "\n", "Here we've created a function that performs fast fourier transforms and then plot's the spectrum for signals of various lengths." ] }, { "cell_type": "code", "execution_count": null, "id": "d2b002d4", "metadata": {}, "outputs": [], "source": [ "def fft_and_plot(df, plot_dominant_periods=4):\n", " N = len(df)\n", " yf = np.fft.fft(df.detrended)\n", " yf = np.abs(yf[:int(N / 2)])\n", " \n", " # get the right frequency \n", " # https://docs.scipy.org/doc/numpy/reference/generated/numpy.fft.fftfreq.html#numpy.fft.fftfreq\n", " d = 1. # day\n", " f = np.fft.fftfreq(N, d)[:int(N/2)]\n", " f[0] = .00001\n", " per = 1./f # days\n", " \n", " fig = plt.figure(figsize=(12,6))\n", " ax = plt.subplot(2,1,1)\n", " plt.plot(per, yf)\n", " plt.xscale('log')\n", "\n", " top=np.argsort(yf)[-plot_dominant_periods:]\n", " j=(10-plot_dominant_periods)/10\n", " for i in top:\n", " plt.plot([per[i], per[i]], [0,np.max(yf)], 'r:')\n", " plt.text(per[i], j*np.max(yf), f\"{per[i] :.2f}\")\n", " j+=0.1\n", "\n", " plt.title('Period Spectrum')\n", " plt.grid()\n", " ax.set_xlabel('Period (days)')\n", " plt.xlim([1, 1e4])\n", "\n", " plt.subplot(2,1,2)\n", " plt.plot(df.index,df.Q)\n", " plt.title('Raw Signal')\n", " plt.tight_layout()" ] }, { "cell_type": "markdown", "id": "89c462a3", "metadata": {}, "source": [ "### Now let's look at the whole signal" ] }, { "cell_type": "code", "execution_count": null, "id": "ae061a46", "metadata": {}, "outputs": [], "source": [ "fft_and_plot(df_data)" ] }, { "cell_type": "markdown", "id": "35e769e4", "metadata": {}, "source": [ "We see a strong annual frequency that corresponds to spring snowmelt from the surrounding mountain ranges. The second strongest frequency is biannually, and is likely due to the onset of fall in winter rains" ] }, { "cell_type": "markdown", "id": "92936f20", "metadata": {}, "source": [ "### Okay, what about years before dams were installed on the Russian River.\n", "\n", "Let's get in the wayback machine and look at only the years prior to 1954!" ] }, { "cell_type": "code", "execution_count": null, "id": "ed01ace9", "metadata": {}, "outputs": [], "source": [ "fft_and_plot(df_data[df_data.index.year < 1954])" ] }, { "cell_type": "markdown", "id": "5ae5807f", "metadata": {}, "source": [ "The binanual and annual signals have a much wider spread, which suggests that peak runoff was a little more variable in it's timing compared to the entire data record." ] }, { "cell_type": "markdown", "id": "238d70b2", "metadata": {}, "source": [ "### What if we focus in on after 1954?" ] }, { "cell_type": "code", "execution_count": null, "id": "b71e84e3", "metadata": {}, "outputs": [], "source": [ "fft_and_plot(df_data[(df_data.index.year > 1954)])" ] }, { "cell_type": "markdown", "id": "e2e13b0a", "metadata": {}, "source": [ "It looks like flows are more controlled with a tighter periods.\n", "\n", "So post dam construction, flows are much more controlled with regular release schedules. " ] } ], "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.11" } }, "nbformat": 4, "nbformat_minor": 5 }