08: Working with tabular data in Pandas
Not that type of Panda – Python’s Pandas package
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.
Let’s begin by importing the libraries and setting our data path
[1]:
from pathlib import Path
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import pandas as pd
import statsmodels.api as sm
from dataretrieval import nwis
from matplotlib.backends.backend_pdf import PdfPages
from scipy.signal import detrend
data_path = Path("..", "data", "pandas")
Get site information data for part of the Russian River near Guerneville, CA from NWIS
[2]:
bbox = (-123.10, 38.45, -122.90, 38.55)
info, metadata = nwis.get_info(bBox=[str(i) for i in bbox])
if 'geometry' in list(info):
# dataretreival now returns as geopandas dataframe, drop "geometry" column for this example
info.drop(columns=["geometry"], inplace=True)
---------------------------------------------------------------------------
SSLEOFError Traceback (most recent call last)
File ~/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/urllib3/connectionpool.py:464, in HTTPConnectionPool._make_request(self, conn, method, url, body, headers, retries, timeout, chunked, response_conn, preload_content, decode_content, enforce_content_length)
463 try:
--> 464 self._validate_conn(conn)
465 except (SocketTimeout, BaseSSLError) as e:
File ~/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/urllib3/connectionpool.py:1093, in HTTPSConnectionPool._validate_conn(self, conn)
1092 if conn.is_closed:
-> 1093 conn.connect()
1095 # TODO revise this, see https://github.com/urllib3/urllib3/issues/2791
File ~/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/urllib3/connection.py:790, in HTTPSConnection.connect(self)
788 server_hostname_rm_dot = server_hostname.rstrip(".")
--> 790 sock_and_verified = _ssl_wrap_socket_and_match_hostname(
791 sock=sock,
792 cert_reqs=self.cert_reqs,
793 ssl_version=self.ssl_version,
794 ssl_minimum_version=self.ssl_minimum_version,
795 ssl_maximum_version=self.ssl_maximum_version,
796 ca_certs=self.ca_certs,
797 ca_cert_dir=self.ca_cert_dir,
798 ca_cert_data=self.ca_cert_data,
799 cert_file=self.cert_file,
800 key_file=self.key_file,
801 key_password=self.key_password,
802 server_hostname=server_hostname_rm_dot,
803 ssl_context=self.ssl_context,
804 tls_in_tls=tls_in_tls,
805 assert_hostname=self.assert_hostname,
806 assert_fingerprint=self.assert_fingerprint,
807 )
808 self.sock = sock_and_verified.socket
File ~/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/urllib3/connection.py:969, in _ssl_wrap_socket_and_match_hostname(sock, cert_reqs, ssl_version, ssl_minimum_version, ssl_maximum_version, cert_file, key_file, key_password, ca_certs, ca_cert_dir, ca_cert_data, assert_hostname, assert_fingerprint, server_hostname, ssl_context, tls_in_tls)
967 server_hostname = normalized
--> 969 ssl_sock = ssl_wrap_socket(
970 sock=sock,
971 keyfile=key_file,
972 certfile=cert_file,
973 key_password=key_password,
974 ca_certs=ca_certs,
975 ca_cert_dir=ca_cert_dir,
976 ca_cert_data=ca_cert_data,
977 server_hostname=server_hostname,
978 ssl_context=context,
979 tls_in_tls=tls_in_tls,
980 )
982 try:
File ~/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/urllib3/util/ssl_.py:480, in ssl_wrap_socket(sock, keyfile, certfile, cert_reqs, ca_certs, server_hostname, ssl_version, ciphers, ssl_context, ca_cert_dir, key_password, ca_cert_data, tls_in_tls)
478 context.set_alpn_protocols(ALPN_PROTOCOLS)
--> 480 ssl_sock = _ssl_wrap_socket_impl(sock, context, tls_in_tls, server_hostname)
481 return ssl_sock
File ~/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/urllib3/util/ssl_.py:524, in _ssl_wrap_socket_impl(sock, ssl_context, tls_in_tls, server_hostname)
522 return SSLTransport(sock, ssl_context, server_hostname)
--> 524 return ssl_context.wrap_socket(sock, server_hostname=server_hostname)
File ~/micromamba/envs/pyclass-docs/lib/python3.11/ssl.py:517, in SSLContext.wrap_socket(self, sock, server_side, do_handshake_on_connect, suppress_ragged_eofs, server_hostname, session)
511 def wrap_socket(self, sock, server_side=False,
512 do_handshake_on_connect=True,
513 suppress_ragged_eofs=True,
514 server_hostname=None, session=None):
515 # SSLSocket class handles server_hostname encoding before it calls
516 # ctx._wrap_socket()
--> 517 return self.sslsocket_class._create(
518 sock=sock,
519 server_side=server_side,
520 do_handshake_on_connect=do_handshake_on_connect,
521 suppress_ragged_eofs=suppress_ragged_eofs,
522 server_hostname=server_hostname,
523 context=self,
524 session=session
525 )
File ~/micromamba/envs/pyclass-docs/lib/python3.11/ssl.py:1104, in SSLSocket._create(cls, sock, server_side, do_handshake_on_connect, suppress_ragged_eofs, server_hostname, context, session)
1103 raise ValueError("do_handshake_on_connect should not be specified for non-blocking sockets")
-> 1104 self.do_handshake()
1105 except:
File ~/micromamba/envs/pyclass-docs/lib/python3.11/ssl.py:1382, in SSLSocket.do_handshake(self, block)
1381 self.settimeout(None)
-> 1382 self._sslobj.do_handshake()
1383 finally:
SSLEOFError: [SSL: UNEXPECTED_EOF_WHILE_READING] EOF occurred in violation of protocol (_ssl.c:1016)
During handling of the above exception, another exception occurred:
SSLError Traceback (most recent call last)
File ~/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/urllib3/connectionpool.py:787, in HTTPConnectionPool.urlopen(self, method, url, body, headers, retries, redirect, assert_same_host, timeout, pool_timeout, release_conn, chunked, body_pos, preload_content, decode_content, **response_kw)
786 # Make the request on the HTTPConnection object
--> 787 response = self._make_request(
788 conn,
789 method,
790 url,
791 timeout=timeout_obj,
792 body=body,
793 headers=headers,
794 chunked=chunked,
795 retries=retries,
796 response_conn=response_conn,
797 preload_content=preload_content,
798 decode_content=decode_content,
799 **response_kw,
800 )
802 # Everything went great!
File ~/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/urllib3/connectionpool.py:488, in HTTPConnectionPool._make_request(self, conn, method, url, body, headers, retries, timeout, chunked, response_conn, preload_content, decode_content, enforce_content_length)
487 new_e = _wrap_proxy_error(new_e, conn.proxy.scheme)
--> 488 raise new_e
490 # conn.request() calls http.client.*.request, not the method in
491 # urllib3.request. It also calls makefile (recv) on the socket.
SSLError: [SSL: UNEXPECTED_EOF_WHILE_READING] EOF occurred in violation of protocol (_ssl.c:1016)
The above exception was the direct cause of the following exception:
MaxRetryError Traceback (most recent call last)
File ~/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/requests/adapters.py:644, in HTTPAdapter.send(self, request, stream, timeout, verify, cert, proxies)
643 try:
--> 644 resp = conn.urlopen(
645 method=request.method,
646 url=url,
647 body=request.body,
648 headers=request.headers,
649 redirect=False,
650 assert_same_host=False,
651 preload_content=False,
652 decode_content=False,
653 retries=self.max_retries,
654 timeout=timeout,
655 chunked=chunked,
656 )
658 except (ProtocolError, OSError) as err:
File ~/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/urllib3/connectionpool.py:841, in HTTPConnectionPool.urlopen(self, method, url, body, headers, retries, redirect, assert_same_host, timeout, pool_timeout, release_conn, chunked, body_pos, preload_content, decode_content, **response_kw)
839 new_e = ProtocolError("Connection aborted.", new_e)
--> 841 retries = retries.increment(
842 method, url, error=new_e, _pool=self, _stacktrace=sys.exc_info()[2]
843 )
844 retries.sleep()
File ~/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/urllib3/util/retry.py:519, in Retry.increment(self, method, url, response, error, _pool, _stacktrace)
518 reason = error or ResponseError(cause)
--> 519 raise MaxRetryError(_pool, url, reason) from reason # type: ignore[arg-type]
521 log.debug("Incremented Retry for (url='%s'): %r", url, new_retry)
MaxRetryError: HTTPSConnectionPool(host='waterservices.usgs.gov', port=443): Max retries exceeded with url: /nwis/site?bBox=-123.1%2C38.45%2C-122.9%2C38.55&siteOutput=Expanded&format=rdb (Caused by SSLError(SSLEOFError(8, '[SSL: UNEXPECTED_EOF_WHILE_READING] EOF occurred in violation of protocol (_ssl.c:1016)')))
During handling of the above exception, another exception occurred:
SSLError Traceback (most recent call last)
Cell In[2], line 2
1 bbox = (-123.10, 38.45, -122.90, 38.55)
----> 2 info, metadata = nwis.get_info(bBox=[str(i) for i in bbox])
3 if 'geometry' in list(info):
4 # dataretreival now returns as geopandas dataframe, drop "geometry" column for this example
5 info.drop(columns=["geometry"], inplace=True)
File ~/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/dataretrieval/nwis.py:816, in get_info(ssl_check, **kwargs)
812 else:
813 # cannot have both seriesCatalogOutput and the expanded format
814 kwargs["siteOutput"] = "Expanded"
--> 816 response = query_waterservices("site", ssl_check=ssl_check, **kwargs)
818 return _read_rdb(response.text), NWIS_Metadata(response, **kwargs)
File ~/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/dataretrieval/nwis.py:633, in query_waterservices(service, ssl_check, **kwargs)
629 kwargs["format"] = "rdb"
631 url = WATERSERVICE_URL + service
--> 633 return query(url, payload=kwargs, ssl_check=ssl_check)
File ~/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/dataretrieval/utils.py:190, in query(url, payload, delimiter, ssl_check)
183 # for index in range(len(payload)):
184 # key, value = payload[index]
185 # payload[index] = (key, to_str(value))
186
187 # define the user agent for the query
188 user_agent = {"user-agent": f"python-dataretrieval/{dataretrieval.__version__}"}
--> 190 response = requests.get(url, params=payload, headers=user_agent, verify=ssl_check)
192 if response.status_code == 400:
193 raise ValueError(
194 f"Bad Request, check that your parameters are correct. URL: {response.url}"
195 )
File ~/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/requests/api.py:73, in get(url, params, **kwargs)
62 def get(url, params=None, **kwargs):
63 r"""Sends a GET request.
64
65 :param url: URL for the new :class:`Request` object.
(...) 70 :rtype: requests.Response
71 """
---> 73 return request("get", url, params=params, **kwargs)
File ~/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/requests/api.py:59, in request(method, url, **kwargs)
55 # By using the 'with' statement we are sure the session is closed, thus we
56 # avoid leaving sockets open which can trigger a ResourceWarning in some
57 # cases, and look like a memory leak in others.
58 with sessions.Session() as session:
---> 59 return session.request(method=method, url=url, **kwargs)
File ~/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/requests/sessions.py:589, in Session.request(self, method, url, params, data, headers, cookies, files, auth, timeout, allow_redirects, proxies, hooks, stream, verify, cert, json)
584 send_kwargs = {
585 "timeout": timeout,
586 "allow_redirects": allow_redirects,
587 }
588 send_kwargs.update(settings)
--> 589 resp = self.send(prep, **send_kwargs)
591 return resp
File ~/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/requests/sessions.py:703, in Session.send(self, request, **kwargs)
700 start = preferred_clock()
702 # Send the request
--> 703 r = adapter.send(request, **kwargs)
705 # Total elapsed time of the request (approximately)
706 elapsed = preferred_clock() - start
File ~/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/requests/adapters.py:675, in HTTPAdapter.send(self, request, stream, timeout, verify, cert, proxies)
671 raise ProxyError(e, request=request)
673 if isinstance(e.reason, _SSLError):
674 # This branch is for urllib3 v1.22 and later.
--> 675 raise SSLError(e, request=request)
677 raise ConnectionError(e, request=request)
679 except ClosedPoolError as e:
SSLError: HTTPSConnectionPool(host='waterservices.usgs.gov', port=443): Max retries exceeded with url: /nwis/site?bBox=-123.1%2C38.45%2C-122.9%2C38.55&siteOutput=Expanded&format=rdb (Caused by SSLError(SSLEOFError(8, '[SSL: UNEXPECTED_EOF_WHILE_READING] EOF occurred in violation of protocol (_ssl.c:1016)')))
[3]:
info.to_csv(data_path / "site_info.csv", index=False)
info.head()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[3], line 1
----> 1 info.to_csv(data_path / "site_info.csv", index=False)
2 info.head()
NameError: name 'info' is not defined
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
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.
Viewing data in pandas
Pandas has built in methods to inspect DataFrame
objects. We’ll look at a few handy methods:
.head()
: inspect the first few rows of data.tail()
: inspect the last few rows of data.index()
: show the row indexes.columns()
: show the column names.describe()
: statistically describe the data
[4]:
info.head()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[4], line 1
----> 1 info.head()
NameError: name 'info' is not defined
[5]:
info.tail()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[5], line 1
----> 1 info.tail()
NameError: name 'info' is not defined
[6]:
# print the index names
print(info.index)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[6], line 2
1 # print the index names
----> 2 print(info.index)
NameError: name 'info' is not defined
[7]:
# print the column names
print(info.columns)
# print the column names as a list
print(list(info))
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[7], line 2
1 # print the column names
----> 2 print(info.columns)
4 # print the column names as a list
5 print(list(info))
NameError: name 'info' is not defined
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
[8]:
info.describe()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[8], line 1
----> 1 info.describe()
NameError: name 'info' is not defined
Getting data from a pandas dataframe
There are multiple methods to get data out of a pandas dataframe as either a “series”, numpy array, or a list
Let’s start by getting data as a series using a few methods
[9]:
# get a series of site numbers by key
info["site_no"]
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[9], line 2
1 # get a series of site numbers by key
----> 2 info["site_no"]
NameError: name 'info' is not defined
[10]:
# get station names by attribute
info.station_nm
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[10], line 2
1 # get station names by attribute
----> 2 info.station_nm
NameError: name 'info' is not defined
getting data from a dataframe as a numpy array can be accomplished by using .values
[11]:
info.station_nm.values
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[11], line 1
----> 1 info.station_nm.values
NameError: name 'info' is not defined
Selection by position
We can get data by position in the dataframe using the .iloc
attribute
[12]:
info.iloc[0:2, 1:4]
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[12], line 1
----> 1 info.iloc[0:2, 1:4]
NameError: name 'info' is not defined
Selection by label
Pandas allows the user to get data from the dataframe by index and column labels
[13]:
info.loc[0:3, ["site_no", "station_nm", "site_tp_cd"]]
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[13], line 1
----> 1 info.loc[0:3, ["site_no", "station_nm", "site_tp_cd"]]
NameError: name 'info' is not defined
Boolean indexing
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.
Let’s get a dataframe of only groundwater sites from the info
dataframe
[14]:
dfgw = info[info["site_tp_cd"] == "GW"]
dfgw
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[14], line 1
----> 1 dfgw = info[info["site_tp_cd"] == "GW"]
2 dfgw
NameError: name 'info' is not defined
Reading and writing data to .csv
files
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 tutorial documentation
For this part we’ll write a new .csv
file of the groundwater sites that we found in NWIS using to_csv()
.
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
[15]:
csv_file = data_path / "RussianRiverGWsites.csv"
dfgw.to_csv(csv_file, index=False)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[15], line 2
1 csv_file = data_path / "RussianRiverGWsites.csv"
----> 2 dfgw.to_csv(csv_file, index=False)
NameError: name 'dfgw' is not defined
Now we can load the csv file back into a pandas dataframe with the read_csv()
method.
[16]:
dfgw2 = pd.read_csv(csv_file)
dfgw2
[16]:
agency_cd | site_no | station_nm | site_tp_cd | lat_va | long_va | dec_lat_va | dec_long_va | coord_meth_cd | coord_acy_cd | ... | local_time_fg | reliability_cd | gw_file_cd | nat_aqfr_cd | aqfr_cd | aqfr_type_cd | well_depth_va | hole_depth_va | depth_src_cd | project_no | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | USGS | 382713123030401 | 007N011W14NU01M | GW | 382712.5 | 1230303.5 | 38.453472 | -123.050972 | G | 5 | ... | Y | C | YY Y | NaN | NaN | NaN | 81.70 | NaN | S | GGMC200 |
1 | USGS | 382713123030403 | 007N011W14NU03M | GW | 382712.5 | 1230303.5 | 38.453472 | -123.050972 | G | 5 | ... | Y | C | YY Y | NaN | NaN | NaN | 32.56 | NaN | S | GGMC200 |
2 | USGS | 382808122565401 | 007N010W10H001M | GW | 382808.5 | 1225654.2 | 38.469028 | -122.948389 | G | 5 | ... | Y | C | YY Y | NaN | NaN | NaN | 22.00 | 22.0 | D | 00BHM8187 |
3 | USGS | 382819123010001 | 007N010W07D001M | GW | 382819.9 | 1230100.8 | 38.472194 | -123.016889 | G | 5 | ... | Y | C | YY | N100CACSTL | NaN | NaN | 99.00 | 110.0 | D | 9677BHM32 |
4 | USGS | 383001122540701 | 008N009W31C002M | GW | 383001.9 | 1225407.7 | 38.500528 | -122.902139 | G | 5 | ... | Y | C | YY | N100CACSTL | NaN | NaN | 80.00 | 80.0 | D | 9677BHM32 |
5 | USGS | 383003122540401 | 008N009W31C003M | GW | 383003.3 | 1225404.4 | 38.500917 | -122.901222 | G | 1 | ... | Y | C | YY Y | N100CACSTL | NaN | NaN | NaN | NaN | NaN | 967760420 |
6 | USGS | 383003122540402 | 008N009W31C004M | GW | 383003.3 | 1225404.4 | 38.500917 | -122.901222 | G | 1 | ... | Y | C | Y | N100CACSTL | NaN | NaN | NaN | NaN | NaN | 967760420 |
7 | USGS | 383003122540403 | 008N009W31C005M | GW | 383003.3 | 1225404.4 | 38.500917 | -122.901222 | G | 1 | ... | Y | C | YY Y | N100CACSTL | NaN | NaN | 25.00 | NaN | S | 967760420 |
8 | USGS | 383016123041101 | 008N011W27N001M | GW | 383016.3 | 1230411.0 | 38.504528 | -123.069722 | G | 5 | ... | Y | C | YY Y | NaN | NaN | NaN | 30.00 | 30.0 | D | 00BHM8187 |
9 | USGS | 383017122580301 | 008N010W28RU01M | GW | 383017.4 | 1225802.5 | 38.504833 | -122.967361 | G | 5 | ... | Y | C | YY | NaN | NaN | NaN | 98.00 | 107.0 | D | NaN |
10 | USGS | 383034122590701 | 008N010W29H004M | GW | 383034.4 | 1225907.9 | 38.509556 | -122.985528 | G | 5 | ... | Y | C | YY Y | N100CACSTL | NaN | NaN | 99.00 | 110.0 | D | 9677BHM32 |
11 rows × 42 columns
Class exercise 1
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
[17]:
csv_name = data_path / "RussianRiverSWsites.csv"
dfsw = info[info.site_tp_cd == "ST"]
dfsw.to_csv(csv_name, index=False)
dfsw = pd.read_csv(csv_name)
dfsw.head()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[17], line 2
1 csv_name = data_path / "RussianRiverSWsites.csv"
----> 2 dfsw = info[info.site_tp_cd == "ST"]
3 dfsw.to_csv(csv_name, index=False)
5 dfsw = pd.read_csv(csv_name)
NameError: name 'info' is not defined
Get timeseries data from NWIS for surface water sites
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.
(11467002)
[18]:
# filter for only SW sites
dfsw = info[info.site_tp_cd == "ST"]
# create query
sites = dfsw.site_no.to_list()
pcode = "00060"
start_date = "1939-09-22"
end_date = "2023-12-31"
df, metadata = nwis.get_dv(
sites=sites,
start=start_date,
end=end_date,
parameterCd=pcode,
multi_index=False
)
df.head()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[18], line 2
1 # filter for only SW sites
----> 2 dfsw = info[info.site_tp_cd == "ST"]
4 # create query
5 sites = dfsw.site_no.to_list()
NameError: name 'info' is not defined
Awesome! There are two gages in this study area that have daily discharge data!
Inspect the data using .plot()
[19]:
ax = df.plot()
ax.set_ylabel(r"cubic meter per second")
ax.set_yscale("log")
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[19], line 1
----> 1 ax = df.plot()
2 ax.set_ylabel(r"cubic meter per second")
3 ax.set_yscale("log")
NameError: name 'df' is not defined
Updating column names
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.
the rename()
method accepts a dictionary that is formatted {current_col_name: new_col_name}
[20]:
df
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[20], line 1
----> 1 df
NameError: name 'df' is not defined
[21]:
df2 = df[['site_no', '00060_Mean']].pivot(columns=['site_no'])
df2.columns = df2.columns.get_level_values(1)
site_names = {site_no: '_'.join(name.split()[:2]).lower() for site_no, name
in info.groupby('site_no')['station_nm'].first().loc[df2.columns].items()}
df2.rename(columns=site_names, inplace=True)
df2.plot()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[21], line 1
----> 1 df2 = df[['site_no', '00060_Mean']].pivot(columns=['site_no'])
2 df2.columns = df2.columns.get_level_values(1)
4 site_names = {site_no: '_'.join(name.split()[:2]).lower() for site_no, name
5 in info.groupby('site_no')['station_nm'].first().loc[df2.columns].items()}
NameError: name 'df' is not defined
Adding a new column of data to the dataframe
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”
Let’s add a year column (water year) to the dataframe by getting the year from the index and adjusting it using np.where()
[22]:
df2['year'] = np.where(df2.index.month >=10, df2.index.year, df2.index.year - 1)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[22], line 1
----> 1 df2['year'] = np.where(df2.index.month >=10, df2.index.year, df2.index.year - 1)
NameError: name 'df2' is not defined
Manipulating data
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}\).
The conversion for this is:
$ \frac{1 ft^3}{s} \times `:nbsphinx-math:frac{60s}{min}` \times `:nbsphinx-math:frac{60min}{hour}` \times `:nbsphinx-math:frac{24hours}{day}` \rightarrow `:nbsphinx-math:frac{ft^3}{day}` $
and then convert that to acre-feet (1 acre-ft = 43559.9)
[23]:
df2
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[23], line 1
----> 1 df2
NameError: name 'df2' is not defined
[24]:
conv = 60 * 60 * 24
df2["russian_r_cfd"] = df2.russian_r * conv
df2["russian_r_af"] = df2.russian_r_cfd / 43559.9
df2
# now let's plot up discharge from 2018 - 2020
dfplt = df2[(df2["year"] >= 2015) & (df2["year"] <= 2020)]
ax = dfplt["russian_r_af"].plot()
ax.set_yscale("log")
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[24], line 2
1 conv = 60 * 60 * 24
----> 2 df2["russian_r_cfd"] = df2.russian_r * conv
3 df2["russian_r_af"] = df2.russian_r_cfd / 43559.9
4 df2
NameError: name 'df2' is not defined
Class exercise 2:
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"
.
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”.
bonus exercise try to plot both "russian_r_af"
and "austin_c_af"
on the same plot
[25]:
conv = 60 * 60 * 24
df2["austin_c_cfd"] = df2.austin_c * conv
df2["austin_c_af"] = df2.austin_c_cfd / 43559.9
df2
# now let's plot up discharge from 2018 - 2020
dfplt = df2[(df2["year"] >= 2015) & (df2["year"] <= 2020)]
ax = dfplt["austin_c_af"].plot()
ax.set_yscale("log")
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[25], line 2
1 conv = 60 * 60 * 24
----> 2 df2["austin_c_cfd"] = df2.austin_c * conv
3 df2["austin_c_af"] = df2.austin_c_cfd / 43559.9
4 df2
NameError: name 'df2' is not defined
[26]:
# now let's plot up discharge from 2018 - 2020
dfplt = df2[(df2["year"] >= 2015) & (df2["year"] <= 2020)]
ax = dfplt[["russian_r_af","austin_c_af"]].plot()
ax.set_yscale("log")
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[26], line 2
1 # now let's plot up discharge from 2018 - 2020
----> 2 dfplt = df2[(df2["year"] >= 2015) & (df2["year"] <= 2020)]
3 ax = dfplt[["russian_r_af","austin_c_af"]].plot()
4 ax.set_yscale("log")
NameError: name 'df2' is not defined
groupby: grouping data and performing mathematical operations on it
groupby()
is a powerful method that allows for performing statistical operations over a groups of “common” data within a dataframe.
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.
[27]:
df2["day_of_year"] = df2.index.day_of_year
df_mean_day = df2.groupby(by=["day_of_year"], as_index=False)[["russian_r_af", "austin_c_af"]].mean()
ax = df_mean_day[["russian_r_af", "austin_c_af"]].plot()
ax.set_yscale("log")
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[27], line 1
----> 1 df2["day_of_year"] = df2.index.day_of_year
3 df_mean_day = df2.groupby(by=["day_of_year"], as_index=False)[["russian_r_af", "austin_c_af"]].mean()
4 ax = df_mean_day[["russian_r_af", "austin_c_af"]].plot()
NameError: name 'df2' is not defined
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.
[28]:
cimis_file = data_path / "santa_rosa_CIMIS_83.csv"
df_cimis = pd.read_csv(cimis_file)
drop_list = ["qc"] + [f"qc.{i}" for i in range(1, 22)]
df_cimis.drop(columns=drop_list, inplace=True)
df_cimis
[28]:
Stn Id | Stn Name | CIMIS Region | Date | Jul | ETo (in) | Precip (in) | Sol Rad (Ly/day) | Avg Vap Pres (mBars) | Max Air Temp (F) | ... | Rose NNE | Rose ENE | Rose ESE | Rose SSE | Rose SSW | Rose WSW | Rose WNW | Rose NNW | Wind Run (miles) | Avg Soil Temp (F) | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 83 | Santa Rosa | North Coast Valleys | 1/1/1990 | 1 | 0.00 | 0.30 | 42.0 | 7.7 | 48.7 | ... | 0.03 | 0.02 | 0.06 | 0.06 | 0.04 | 0.10 | 0.15 | 0.14 | 32.1 | 46.3 |
1 | 83 | Santa Rosa | North Coast Valleys | 1/2/1990 | 2 | 0.06 | 0.00 | 236.0 | 6.1 | 58.4 | ... | 0.06 | 0.01 | 0.04 | 0.03 | 0.03 | 0.06 | 0.45 | 1.53 | 121.0 | 44.9 |
2 | 83 | Santa Rosa | North Coast Valleys | 1/3/1990 | 3 | 0.04 | 0.00 | 225.0 | 5.8 | 63.0 | ... | 0.08 | 0.04 | 0.10 | 0.04 | 0.06 | 0.05 | 0.21 | 0.23 | 44.5 | 43.3 |
3 | 83 | Santa Rosa | North Coast Valleys | 1/4/1990 | 4 | 0.03 | 0.01 | 190.0 | 6.3 | 56.7 | ... | 0.07 | 0.07 | 0.11 | 0.07 | 0.03 | 0.03 | 0.08 | 0.09 | 30.2 | 42.8 |
4 | 83 | Santa Rosa | North Coast Valleys | 1/5/1990 | 5 | 0.04 | 0.00 | 192.0 | 7.2 | 59.1 | ... | 0.08 | 0.13 | 0.33 | 0.08 | 0.04 | 0.08 | 0.07 | 0.05 | 45.2 | 43.7 |
... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
11948 | 83 | Santa Rosa | North Coast Valleys | 9/18/2022 | 261 | 0.03 | 0.77 | 113.0 | 16.4 | 66.8 | ... | 0.01 | 0.04 | 0.29 | 0.51 | 0.14 | 0.01 | 0.00 | 0.00 | 197.4 | 65.2 |
11949 | 83 | Santa Rosa | North Coast Valleys | 9/19/2022 | 262 | 0.14 | 0.02 | 424.0 | 15.6 | 73.9 | ... | 0.09 | 0.07 | 0.08 | 0.17 | 0.09 | 0.28 | 0.15 | 0.07 | 108.9 | 65.9 |
11950 | 83 | Santa Rosa | North Coast Valleys | 9/20/2022 | 263 | 0.10 | 0.00 | 330.0 | 13.9 | 72.6 | ... | 0.08 | 0.04 | 0.10 | 0.11 | 0.18 | 0.25 | 0.15 | 0.10 | 79.6 | 65.8 |
11951 | 83 | Santa Rosa | North Coast Valleys | 9/21/2022 | 264 | 0.15 | 0.07 | 458.0 | 14.8 | 75.5 | ... | 0.09 | 0.05 | 0.13 | 0.16 | 0.18 | 0.22 | 0.08 | 0.09 | 103.3 | 66.7 |
11952 | 83 | Santa Rosa | North Coast Valleys | 9/22/2022 | 265 | 0.17 | 0.00 | 494.0 | 13.1 | 82.9 | ... | 0.14 | 0.06 | 0.09 | 0.10 | 0.11 | 0.20 | 0.12 | 0.19 | 83.7 | 66.3 |
11953 rows × 27 columns
[29]:
df_mean_clim = df_cimis.groupby(by=["Jul"], as_index=False)[["ETo (in)", "Precip (in)"]].mean()
df_mean_clim = df_mean_clim.rename(columns={"Jul": "day_of_year"})
df_mean_clim.head()
[29]:
day_of_year | ETo (in) | Precip (in) | |
---|---|---|---|
0 | 1 | 0.031212 | 0.335455 |
1 | 2 | 0.031212 | 0.155152 |
2 | 3 | 0.030303 | 0.195484 |
3 | 4 | 0.030909 | 0.190303 |
4 | 5 | 0.034242 | 0.150000 |
Joining two DataFrames
with an inner join
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.
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
[30]:
df_merged = pd.merge(df_mean_day, df_mean_clim, on=["day_of_year",])
df_merged.head()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[30], line 1
----> 1 df_merged = pd.merge(df_mean_day, df_mean_clim, on=["day_of_year",])
2 df_merged.head()
NameError: name 'df_mean_day' is not defined
[31]:
cols = [i for i in list(df_merged) if i != "day_of_year"]
ax = df_merged[cols].plot()
ax.set_yscale("log")
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[31], line 1
----> 1 cols = [i for i in list(df_merged) if i != "day_of_year"]
2 ax = df_merged[cols].plot()
3 ax.set_yscale("log")
NameError: name 'df_merged' is not defined
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.
Stepping back for a minute: exploring Merge and Concat
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.
The merge
method
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.
inner
: operates as a intersection of the data and keeps rows from both dataframes if the row has a matching key in both dataframes.outer
: operates as a union and joins all rows from both dataframes, fills columns with nan where data is not available for a row.left
: keeps all rows from the ‘left’ dataframe and joins columns from the rightright
: keeps all rows from the ‘right’ dataframe and joins columns from the left
[32]:
ex0 = {
"name": [],
"childhood_dream_job": []
}
ex1 = {
"name": [],
"age": []
}
[ ]:
[ ]:
[ ]:
The concat
method
Pandas concat
method allows the user to concatenate dataframes together
[33]:
gw = info[info.site_tp_cd == "GW"]
sw = info[info.site_tp_cd == "ST"]
sw = sw[list(sw)[0:8]]
sw.head()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[33], line 1
----> 1 gw = info[info.site_tp_cd == "GW"]
2 sw = info[info.site_tp_cd == "ST"]
3 sw = sw[list(sw)[0:8]]
NameError: name 'info' is not defined
[34]:
newdf = pd.concat((gw, sw), ignore_index=True)
newdf
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[34], line 1
----> 1 newdf = pd.concat((gw, sw), ignore_index=True)
2 newdf
NameError: name 'gw' is not defined
Back to working with our climate and stream data
Is there anything else we can look at that might give us more insights into this basin?
Let’s try looking at long term, yearly trends to see what’s there.
[35]:
year = []
for dstr in df_cimis.Date.values:
mo, d, yr = [int(i) for i in dstr.split("/")]
if mo < 10:
year.append(yr - 1)
else:
year.append(yr)
df_cimis["year"] = year
df_cimis.head()
[35]:
Stn Id | Stn Name | CIMIS Region | Date | Jul | ETo (in) | Precip (in) | Sol Rad (Ly/day) | Avg Vap Pres (mBars) | Max Air Temp (F) | ... | Rose ENE | Rose ESE | Rose SSE | Rose SSW | Rose WSW | Rose WNW | Rose NNW | Wind Run (miles) | Avg Soil Temp (F) | year | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 83 | Santa Rosa | North Coast Valleys | 1/1/1990 | 1 | 0.00 | 0.30 | 42.0 | 7.7 | 48.7 | ... | 0.02 | 0.06 | 0.06 | 0.04 | 0.10 | 0.15 | 0.14 | 32.1 | 46.3 | 1989 |
1 | 83 | Santa Rosa | North Coast Valleys | 1/2/1990 | 2 | 0.06 | 0.00 | 236.0 | 6.1 | 58.4 | ... | 0.01 | 0.04 | 0.03 | 0.03 | 0.06 | 0.45 | 1.53 | 121.0 | 44.9 | 1989 |
2 | 83 | Santa Rosa | North Coast Valleys | 1/3/1990 | 3 | 0.04 | 0.00 | 225.0 | 5.8 | 63.0 | ... | 0.04 | 0.10 | 0.04 | 0.06 | 0.05 | 0.21 | 0.23 | 44.5 | 43.3 | 1989 |
3 | 83 | Santa Rosa | North Coast Valleys | 1/4/1990 | 4 | 0.03 | 0.01 | 190.0 | 6.3 | 56.7 | ... | 0.07 | 0.11 | 0.07 | 0.03 | 0.03 | 0.08 | 0.09 | 30.2 | 42.8 | 1989 |
4 | 83 | Santa Rosa | North Coast Valleys | 1/5/1990 | 5 | 0.04 | 0.00 | 192.0 | 7.2 | 59.1 | ... | 0.13 | 0.33 | 0.08 | 0.04 | 0.08 | 0.07 | 0.05 | 45.2 | 43.7 | 1989 |
5 rows × 28 columns
[36]:
df2.head()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[36], line 1
----> 1 df2.head()
NameError: name 'df2' is not defined
.aggregate()
The pandas .aggregate()
can be used with .groupy()
to perform multiple statistical operations.
Example usage could be:
agdf = df2.groupby(by=["day_of_year"], as_index=False)["austin_c_af"].aggregate([np.min, np.max, np.mean, np.std])
And would return a dataframe grouped by the day of year and take the min, max, mean, and standard deviation of these data.
Class exercise 3: groupy, aggregate, and merge
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
Hints:
df2
anddf_cimis
are the dataframes these operations should be performed onuse
.groupby()
to group by the year and doprecip
andETo
in separate groupby routinesrename the columns in the ETo aggregated dataframe and the Precip dataframe to
mean_et
,std_et
,mean_precip
,std_precip
np.mean
andnp.std
can be used to calculate mean and std flows
Name your final joined dataframe df_yearly
[37]:
dfg0 = df2.groupby(by=["year"], as_index=False)["russian_r_af"].aggregate(['mean', 'std'])
dfg1 = df_cimis.groupby(by=["year"], as_index=False)["ETo (in)"].aggregate(['mean', 'std'])
dfg2 = df_cimis.groupby(by=["year"], as_index=False)["Precip (in)"].aggregate(['mean', 'std'])
dfm1 = pd.merge(dfg0, dfg1, on=["year"], suffixes=(None, "_et"))
df_yearly = pd.merge(dfm1, dfg2, on=["year"], suffixes=(None, "_precip"))
df_yearly
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[37], line 1
----> 1 dfg0 = df2.groupby(by=["year"], as_index=False)["russian_r_af"].aggregate(['mean', 'std'])
2 dfg1 = df_cimis.groupby(by=["year"], as_index=False)["ETo (in)"].aggregate(['mean', 'std'])
3 dfg2 = df_cimis.groupby(by=["year"], as_index=False)["Precip (in)"].aggregate(['mean', 'std'])
NameError: name 'df2' is not defined
Lets examine the long term flow record for the Russian River
[38]:
df_y_flow = df2.groupby(by=["year"], as_index=False)["russian_r_af"].aggregate(['mean', 'std'])
fig = plt.figure(figsize=(12,4))
lower_ci = df_y_flow["mean"] - 2 * df_y_flow['std']
lower_ci = np.where(lower_ci < 0, 0, lower_ci)
upper_ci = df_y_flow["mean"] + 2 * df_y_flow['std']
ax = df_y_flow["mean"].plot(style="b.-")
ax.fill_between(df_y_flow.index, lower_ci, upper_ci, color="b", alpha=0.5);
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[38], line 1
----> 1 df_y_flow = df2.groupby(by=["year"], as_index=False)["russian_r_af"].aggregate(['mean', 'std'])
3 fig = plt.figure(figsize=(12,4))
5 lower_ci = df_y_flow["mean"] - 2 * df_y_flow['std']
NameError: name 'df2' is not defined
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.
First let’s see if there are any relationships between yearly discharge and climate
[39]:
dfg0 = df2.groupby(by=["year"], as_index=False)["russian_r_af"].aggregate(['mean', 'std', 'sum'])
dfg1 = df_cimis.groupby(by=["year"], as_index=False)["ETo (in)"].aggregate(['mean', 'std', 'sum'])
dfg2 = df_cimis.groupby(by=["year"], as_index=False)["Precip (in)"].aggregate(['mean', 'std', 'sum'])
dfm1 = pd.merge(dfg0, dfg1, on=["year"], suffixes=(None, "_et"))
df_yearly = pd.merge(dfm1, dfg2, on=["year"], suffixes=(None, "_precip"))
fig = plt.figure(figsize=(12,4))
lower_ci = df_yearly["mean"] - 2 * df_yearly['std']
lower_ci = np.where(lower_ci < 0, 0, lower_ci)
upper_ci = df_yearly["mean"] + 2 * df_yearly['std']
ax = df_yearly["mean"].plot(style="b.-", label="flow_af")
ax.fill_between(df_yearly.index, lower_ci, upper_ci, color="b", alpha=0.5)
ax2 = ax.twinx()
ax2.plot(df_yearly.index, df_yearly.sum_precip, "k--", lw=1.5, label="Precip")
ax2.plot(df_yearly.index, df_yearly.sum_et, "r.--", lw=1.5, label="ET")
plt.legend()
plt.show();
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[39], line 1
----> 1 dfg0 = df2.groupby(by=["year"], as_index=False)["russian_r_af"].aggregate(['mean', 'std', 'sum'])
2 dfg1 = df_cimis.groupby(by=["year"], as_index=False)["ETo (in)"].aggregate(['mean', 'std', 'sum'])
3 dfg2 = df_cimis.groupby(by=["year"], as_index=False)["Precip (in)"].aggregate(['mean', 'std', 'sum'])
NameError: name 'df2' is not defined
As expected, precipitation amount is the main driver of the yearly discharge regime here
Baseflow separation
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
[40]:
def _baseflow_low_pass_filter(arr, beta, enforce):
"""
Private method to apply digital baseflow separation filter
(Lyne & Hollick, 1979; Nathan & Mcmahon, 1990;
Boughton, 1993; Chapman & Maxwell, 1996).
This method should not be called by the user!
Method removes "spikes" which would be consistent with storm flow
events from data.
Parameters
----------
arr : np.array
streamflow or municipal pumping time series
beta : float
baseflow filtering parameter that ranges from 0 - 1
values in 0.8 - 0.95 range used in literature for
streamflow baseflow separation
enforce : bool
enforce physical constraint of baseflow less than measured flow
Returns
-------
np.array of filtered data
"""
# prepend 10 records to data for initial spin up
# these records will be dropped before returning data to user
qt = np.zeros((arr.size + 10,), dtype=float)
qt[0:10] = arr[0:10]
qt[10:] = arr[:]
qdt = np.zeros((arr.size + 10,), dtype=float)
qbf = np.zeros((arr.size + 10,), dtype=float)
y = (1. + beta) / 2.
for ix in range(qdt.size):
if ix == 0:
qbf[ix] = qt[ix]
continue
x = beta * qdt[ix - 1]
z = qt[ix] - qt[ix - 1]
qdt[ix] = x + (y * z)
qb = qt[ix] - qdt[ix]
if enforce:
if qb > qt[ix]:
qbf[ix] = qt[ix]
else:
qbf[ix] = qb
else:
qbf[ix] = qb
return qbf[10:]
def baseflow_low_pass_filter(arr, beta=0.9, T=1, enforce=True):
"""
User method to apply digital baseflow separation filter
(Lyne & Hollick, 1979; Nathan & Mcmahon, 1990;
Boughton, 1993; Chapman & Maxwell, 1996).
Method removes "spikes" which would be consistent with storm flow
events from data.
Parameters
----------
arr : np.array
streamflow or municipal pumping time series
beta : float
baseflow filtering parameter that ranges from 0 - 1
values in 0.8 - 0.95 range used in literature for
streamflow baseflow separation
T : int
number of recursive filtering passes to apply to the data
enforce : bool
enforce physical constraint of baseflow less than measured flow
Returns
-------
np.array of baseflow
"""
for _ in range(T):
arr = _baseflow_low_pass_filter(arr, beta, enforce)
return arr
Class exercise 4: baseflow separation
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
.
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.
After performing baseflow separation plot both the baseflow and the total discharge for the years 2015 - 2017.
Hints:
make sure to use the
.values
property when you get discharge data from the pandas dataframefeel free to play with the input parameters
beta
andT
.beta
is commonly between 0.8 - 0.95 for hydrologic problems.T
is the number of filter passes over the data, therefore increasingT
will create a smoother profile. I recommend starting with aT
value of around 5 and adjusting it to see how it changes the baseflow profile.
[41]:
baseflow = baseflow_low_pass_filter(df2["russian_r_af"].values, beta=0.90, T=5)
df2["russian_bf_af"] = baseflow
fig = plt.figure(figsize=(15, 10))
tdf = df2[(df2.year >= 2015) & (df2.year <= 2017)]
ax = tdf[["russian_r_af", "russian_bf_af"]].plot()
ax.set_yscale("log");
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[41], line 1
----> 1 baseflow = baseflow_low_pass_filter(df2["russian_r_af"].values, beta=0.90, T=5)
2 df2["russian_bf_af"] = baseflow
4 fig = plt.figure(figsize=(15, 10))
NameError: name 'df2' is not defined
Looking at trends and seasonal cycles in the data with statsmodels
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
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.
[42]:
df2["russian_bf_af"] = baseflow_low_pass_filter(df2["russian_r_af"].values, beta=0.9, T=5)
# drop nan values prior to running decomposition
dfsm = df2[df2['russian_bf_af'].notna()]
dfsm.head()
# decompose
decomposition = sm.tsa.seasonal_decompose(dfsm["russian_bf_af"], period=365, model="additive")
decomposition.plot()
mpl.rcParams['figure.figsize'] = [14, 14];
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[42], line 1
----> 1 df2["russian_bf_af"] = baseflow_low_pass_filter(df2["russian_r_af"].values, beta=0.9, T=5)
2 # drop nan values prior to running decomposition
3 dfsm = df2[df2['russian_bf_af'].notna()]
NameError: name 'df2' is not defined
Now let’s plot our baseflow trend against the baseflow we calculated
[43]:
fig, ax = plt.subplots(figsize=(8, 8))
ax = df2["russian_r_af"].plot(ax=ax)
ax.plot(decomposition.trend.index, decomposition.trend, "r-", lw=1.5)
ax.set_yscale('log')
plt.savefig("russian_river.png")
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[43], line 3
1 fig, ax = plt.subplots(figsize=(8, 8))
----> 3 ax = df2["russian_r_af"].plot(ax=ax)
4 ax.plot(decomposition.trend.index, decomposition.trend, "r-", lw=1.5)
5 ax.set_yscale('log')
NameError: name 'df2' is not defined

What’s going on at the end of the 1970’s? Why is the baseflow trend so low and peak flows are also low?
From “California’s Most Significant Droughts”, Chapter 3
The setting for the 1976-77 drought differed significantly from the dry times of the 1920s-30s. Although only a two-year event, its hydrology was severe. Based on 114 years of computed statewide runoff, 1977 occupies rank 114 (driest year) and 1976 is in rank 104. The drought was notable for the impacts experienced by water agencies that were unprepared for such conditions. One reason for the lack of preparedness was the perception of relatively ample water supplies in most areas of the state.
We can also see the analysis that there is an overall reduction in baseflow from 1940’s to the 2020’s
Writing yearly plots of data to a PDF document
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.
[44]:
# add the decomposition trend to the dataframe to make plotting easy
df2["baseflow_trend"] = decomposition.trend
pdf_file = data_path / "all_years_russian.pdf"
with PdfPages(pdf_file) as outpdf:
for year in sorted(df2.year.unique()):
plt.figure()
tdf = df2[df2.year == year]
tdf[["russian_r_af", "russian_bf_af", "baseflow_trend"]].plot()
plt.xlabel("date")
plt.ylabel("acre-feet")
plt.ylim([1, 300000])
plt.yscale("log")
plt.tight_layout()
outpdf.savefig()
plt.close('all')
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[44], line 2
1 # add the decomposition trend to the dataframe to make plotting easy
----> 2 df2["baseflow_trend"] = decomposition.trend
4 pdf_file = data_path / "all_years_russian.pdf"
5 with PdfPages(pdf_file) as outpdf:
NameError: name 'decomposition' is not defined
Go into the “data/pandas” and open up the file “all_years_russian.pdf” to check out the plots.
Back to the basics:
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
Adjusting a whole column of data values
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
[45]:
df3 = df2.copy()
# Adding values
df3["russian_r"] += 10000
# Subtracting values
df3["russian_r"] -= 10000
# multiplication
df3["russian_r"] *= 5
# division
df3["russian_r"] /= 5
# more complex operations
df3["russian_r"] = np.log10(df3["russian_r"].values)
df3.head()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[45], line 1
----> 1 df3 = df2.copy()
2 # Adding values
3 df3["russian_r"] += 10000
NameError: name 'df2' is not defined
updating values by position
We can update single or multiple values by their position in the dataframe using .iloc
[46]:
df3.iloc[0, 1] = 999
df3.head()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[46], line 1
----> 1 df3.iloc[0, 1] = 999
2 df3.head()
NameError: name 'df3' is not defined
updating values based on location
We can update values in the dataframe based on their index and column headers too using .loc
[47]:
df3.loc[df3.index[0], 'russian_r'] *= 100
df3.head()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[47], line 1
----> 1 df3.loc[df3.index[0], 'russian_r'] *= 100
2 df3.head()
NameError: name 'df3' is not defined
Class exercise 5: setting data
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.
Then use either a location based method to change the value in austin_c on 10-4-1939 to 0.
Finally display the first 5 records in the notebook.
[48]:
df3["russian_r_cfd"] = np.log10(df3["russian_r_cfd"])
df3.loc[df3.index[3], "austin_c"] = 0
df3.head()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[48], line 1
----> 1 df3["russian_r_cfd"] = np.log10(df3["russian_r_cfd"])
2 df3.loc[df3.index[3], "austin_c"] = 0
3 df3.head()
NameError: name 'df3' is not defined
Now for some powerful analysis using fourier analysis to extract the data’s signal.
Some background for your evening viewing: https://youtu.be/spUNpyF58BY
We have to detrend the data for fast fourier transforms to work properly. Here’s a discussion on why:
https://groups.google.com/forum/#!topic/comp.dsp/kYDZqetr_TU
Fortunately we can easily do this in python using scipy!
https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.detrend.html
Let’s create a new dataframe with only the data we’re interested in analysing.
[49]:
df_data = df2[["russian_r_af", "year"]]
df_data = df_data.rename(columns={"russian_r_af": "Q"})
df_data.describe()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[49], line 1
----> 1 df_data = df2[["russian_r_af", "year"]]
2 df_data = df_data.rename(columns={"russian_r_af": "Q"})
3 df_data.describe()
NameError: name 'df2' is not defined
But first let’s also set the water year by date on this dataframe and drop the single nan value from the dataframe
[50]:
df_data['water_year'] = df_data.index.shift(30+31+31,freq='d').year
df_data.dropna(inplace=True)
df_data["detrended"] = detrend(df_data.Q.values)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[50], line 1
----> 1 df_data['water_year'] = df_data.index.shift(30+31+31,freq='d').year
2 df_data.dropna(inplace=True)
3 df_data["detrended"] = detrend(df_data.Q.values)
NameError: name 'df_data' is not defined
[51]:
fig = plt.figure(figsize=(12,6))
ax = plt.subplot(3,1,1)
plt.plot(df_data.detrended)
plt.title('detrended signal')
plt.subplot(3,1,2)
plt.plot(df_data.index, df_data.Q)
plt.title('Raw Signal')
plt.subplot(3,1,3)
plt.plot(df_data.index, df_data.Q - df_data.detrended)
plt.title('Difference')
plt.tight_layout()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[51], line 3
1 fig = plt.figure(figsize=(12,6))
2 ax = plt.subplot(3,1,1)
----> 3 plt.plot(df_data.detrended)
4 plt.title('detrended signal')
6 plt.subplot(3,1,2)
NameError: name 'df_data' is not defined

Evaluate and plot the Period Spectrum to see timing of recurrence
Here we’ve created a function that performs fast fourier transforms and then plot’s the spectrum for signals of various lengths.
[52]:
def fft_and_plot(df, plot_dominant_periods=4):
N = len(df)
yf = np.fft.fft(df.detrended)
yf = np.abs(yf[:int(N / 2)])
# get the right frequency
# https://docs.scipy.org/doc/numpy/reference/generated/numpy.fft.fftfreq.html#numpy.fft.fftfreq
d = 1. # day
f = np.fft.fftfreq(N, d)[:int(N/2)]
f[0] = .00001
per = 1./f # days
fig = plt.figure(figsize=(12,6))
ax = plt.subplot(2,1,1)
plt.plot(per, yf)
plt.xscale('log')
top=np.argsort(yf)[-plot_dominant_periods:]
j=(10-plot_dominant_periods)/10
for i in top:
plt.plot([per[i], per[i]], [0,np.max(yf)], 'r:')
plt.text(per[i], j*np.max(yf), f"{per[i] :.2f}")
j+=0.1
plt.title('Period Spectrum')
plt.grid()
ax.set_xlabel('Period (days)')
plt.xlim([1, 1e4])
plt.subplot(2,1,2)
plt.plot(df.index,df.Q)
plt.title('Raw Signal')
plt.tight_layout()
Now let’s look at the whole signal
[53]:
fft_and_plot(df_data)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[53], line 1
----> 1 fft_and_plot(df_data)
NameError: name 'df_data' is not defined
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
Okay, what about years before dams were installed on the Russian River.
Let’s get in the wayback machine and look at only the years prior to 1954!
[54]:
fft_and_plot(df_data[df_data.index.year < 1954])
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[54], line 1
----> 1 fft_and_plot(df_data[df_data.index.year < 1954])
NameError: name 'df_data' is not defined
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.
What if we focus in on after 1954?
[55]:
fft_and_plot(df_data[(df_data.index.year > 1954)])
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[55], line 1
----> 1 fft_and_plot(df_data[(df_data.index.year > 1954)])
NameError: name 'df_data' is not defined
It looks like flows are more controlled with a tighter periods.
So post dam construction, flows are much more controlled with regular release schedules.