Tutorial 6: Implementing the Analysis

Good Research Practices

Content creators: Marguerite Brown, Zane Mitrevica, Natalie Steinemann, Yuxin Zhou

Content reviewers: Katrina Dobson, Sloane Garelick, Maria Gonzalez, Paul Heubel, Nahid Hasan, Sherry Mi, Beatriz Cosenza Muralles, Cheng Zhang

Content editors: Jenna Pearson, Chi Zhang, Ohad Zivan

Production editors: Wesley Banfield, Paul Heubel, Jenna Pearson, Konstantine Tsafatinos, Chi Zhang, Ohad Zivan

Our 2024 Sponsors: CMIP, NFDI4Earth

Tutorials Objectives

In Tutorials 5-8, you will learn about the research process. This includes how to

5. Draft analyses of data to test a hypothesis

6. Implement analysis of data

7. Interpret results in the context of existing knowledge

8. Communicate your results and conclusions

By the end of these tutorials you will be able to:

• Understand the principles of good research practices

• Learn to view a scientific data set or question through the lens of equity: Who is represented by this data and who is not? Who has access to this information? Who is in a position to use it?

Activity: Implement the Analysis

In this tutorial, you will be implementing a linear regression model as outlined in Step 5 on real-world CO₂ and temperature records.

The CO₂ and temperature records we will be analyzing are both examples of paleoclimate data (for more information, refer back to Step 3). The CO₂ record (Bereiter et al., 2015) was generated by measuring the CO₂ concentration in ancient air bubbles trapped inside ice from multiple ice cores retrieved from Antarctica. The temperature record (Shakun et al., 2015) is based on chemical analysis done on the shells of planktic foraminifera. The foraminifera shells were identified and picked from deep-sea sediments, and the temperature record combined multiple sea-surface temperature records from a range of sites globally.

Why are we focusing on these two records specifically? The CO₂ record from Antarctic ice core is the gold standard of air CO₂ variability on glacial-interglacial time scales, and it has a temporal resolution unmatched by any other reconstruction methods. The temperature record comes from sediment cores all over the global ocean, and therefore is likely representative of the global sea surface temperature (SST) variability. All SST records were shifted to a mean of zero and combined as unweighted global averages. Polar air temperature records are also available from ice core studies, but such records may represent an exaggerated view of the global temperature because of polar amplification.

If you would like to learn more, the data sources are listed at the bottom of the page.

# imports

import matplotlib.pyplot as plt
import pandas as pd
import seaborn as sns
import numpy as np
from scipy import interpolate
from scipy import stats
import os
import pooch
import tempfile

Figure Settings

Show code cell source Hide code cell source

# @title Figure Settings
import ipywidgets as widgets  # interactive display

%config InlineBackend.figure_format = 'retina'
plt.style.use(
    "https://raw.githubusercontent.com/neuromatch/climate-course-content/main/cma.mplstyle"
)

Helper functions

Show code cell source Hide code cell source

# @title Helper functions


def pooch_load(filelocation=None, filename=None, processor=None):
    shared_location = "/home/jovyan/shared/Data/tutorials/W2D1_FutureClimate-IPCCIPhysicalBasis"  # this is different for each day
    user_temp_cache = tempfile.gettempdir()

    if os.path.exists(os.path.join(shared_location, filename)):
        file = os.path.join(shared_location, filename)
    else:
        file = pooch.retrieve(
            filelocation,
            known_hash=None,
            fname=os.path.join(user_temp_cache, filename),
            processor=processor,
        )

    return file

# time series
# read SST data "Shakun2015_SST.txt"
filename_Shakun2015_SST = "Shakun2015_SST.txt"
url_Shakun2015_SST = "https://osf.io/kmy5w/download"
SST = pd.read_table(pooch_load(url_Shakun2015_SST, filename_Shakun2015_SST))
SST.set_index("Age", inplace=True)
SST

Downloading data from 'https://osf.io/kmy5w/download' to file '/tmp/Shakun2015_SST.txt'.

---------------------------------------------------------------------------
TimeoutError                              Traceback (most recent call last)
File ~/micromamba/envs/climatematch/lib/python3.11/site-packages/urllib3/connectionpool.py:534, in HTTPConnectionPool._make_request(self, conn, method, url, body, headers, retries, timeout, chunked, response_conn, preload_content, decode_content, enforce_content_length)
    533 try:
--> 534     response = conn.getresponse()
    535 except (BaseSSLError, OSError) as e:

File ~/micromamba/envs/climatematch/lib/python3.11/site-packages/urllib3/connection.py:571, in HTTPConnection.getresponse(self)
    570 # Get the response from http.client.HTTPConnection
--> 571 httplib_response = super().getresponse()
    573 try:

File ~/micromamba/envs/climatematch/lib/python3.11/http/client.py:1415, in HTTPConnection.getresponse(self)
   1414 try:
-> 1415     response.begin()
   1416 except ConnectionError:

File ~/micromamba/envs/climatematch/lib/python3.11/http/client.py:330, in HTTPResponse.begin(self)
    329 while True:
--> 330     version, status, reason = self._read_status()
    331     if status != CONTINUE:

File ~/micromamba/envs/climatematch/lib/python3.11/http/client.py:291, in HTTPResponse._read_status(self)
    290 def _read_status(self):
--> 291     line = str(self.fp.readline(_MAXLINE + 1), "iso-8859-1")
    292     if len(line) > _MAXLINE:

File ~/micromamba/envs/climatematch/lib/python3.11/socket.py:718, in SocketIO.readinto(self, b)
    717 try:
--> 718     return self._sock.recv_into(b)
    719 except timeout:

File ~/micromamba/envs/climatematch/lib/python3.11/ssl.py:1314, in SSLSocket.recv_into(self, buffer, nbytes, flags)
   1311         raise ValueError(
   1312           "non-zero flags not allowed in calls to recv_into() on %s" %
   1313           self.__class__)
-> 1314     return self.read(nbytes, buffer)
   1315 else:

File ~/micromamba/envs/climatematch/lib/python3.11/ssl.py:1166, in SSLSocket.read(self, len, buffer)
   1165 if buffer is not None:
-> 1166     return self._sslobj.read(len, buffer)
   1167 else:

TimeoutError: The read operation timed out

The above exception was the direct cause of the following exception:

ReadTimeoutError                          Traceback (most recent call last)
File ~/micromamba/envs/climatematch/lib/python3.11/site-packages/requests/adapters.py:696, in HTTPAdapter.send(self, request, stream, timeout, verify, cert, proxies)
    695 try:
--> 696     resp = conn.urlopen(
    697         method=request.method,
    698         url=url,
    699         body=request.body,  # type: ignore[arg-type]  # urllib3 stubs don't accept Iterable[bytes | str]
    700         headers=request.headers,  # type: ignore[arg-type]  # urllib3#3072
    701         redirect=False,
    702         assert_same_host=False,
    703         preload_content=False,
    704         decode_content=False,
    705         retries=self.max_retries,
    706         timeout=resolved_timeout,
    707         chunked=chunked,
    708     )
    710 except (ProtocolError, OSError) as err:

File ~/micromamba/envs/climatematch/lib/python3.11/site-packages/urllib3/connectionpool.py:842, 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)
    840     new_e = ProtocolError("Connection aborted.", new_e)
--> 842 retries = retries.increment(
    843     method, url, error=new_e, _pool=self, _stacktrace=sys.exc_info()[2]
    844 )
    845 retries.sleep()

File ~/micromamba/envs/climatematch/lib/python3.11/site-packages/urllib3/util/retry.py:498, in Retry.increment(self, method, url, response, error, _pool, _stacktrace)
    497 if read is False or method is None or not self._is_method_retryable(method):
--> 498     raise reraise(type(error), error, _stacktrace)
    499 elif read is not None:

File ~/micromamba/envs/climatematch/lib/python3.11/site-packages/urllib3/util/util.py:39, in reraise(tp, value, tb)
     38         raise value.with_traceback(tb)
---> 39     raise value
     40 finally:

File ~/micromamba/envs/climatematch/lib/python3.11/site-packages/urllib3/connectionpool.py:788, 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)
    787 # Make the request on the HTTPConnection object
--> 788 response = self._make_request(
    789     conn,
    790     method,
    791     url,
    792     timeout=timeout_obj,
    793     body=body,
    794     headers=headers,
    795     chunked=chunked,
    796     retries=retries,
    797     response_conn=response_conn,
    798     preload_content=preload_content,
    799     decode_content=decode_content,
    800     **response_kw,
    801 )
    803 # Everything went great!

File ~/micromamba/envs/climatematch/lib/python3.11/site-packages/urllib3/connectionpool.py:536, in HTTPConnectionPool._make_request(self, conn, method, url, body, headers, retries, timeout, chunked, response_conn, preload_content, decode_content, enforce_content_length)
    535 except (BaseSSLError, OSError) as e:
--> 536     self._raise_timeout(err=e, url=url, timeout_value=read_timeout)
    537     raise

File ~/micromamba/envs/climatematch/lib/python3.11/site-packages/urllib3/connectionpool.py:367, in HTTPConnectionPool._raise_timeout(self, err, url, timeout_value)
    366 if isinstance(err, SocketTimeout):
--> 367     raise ReadTimeoutError(
    368         self, url, f"Read timed out. (read timeout={timeout_value})"
    369     ) from err
    371 # See the above comment about EAGAIN in Python 3.

ReadTimeoutError: HTTPSConnectionPool(host='files.ca-1.osf.io', port=443): Read timed out. (read timeout=30)

During handling of the above exception, another exception occurred:

ReadTimeout                               Traceback (most recent call last)
Cell In[4], line 5
      1 # time series
      2 # read SST data "Shakun2015_SST.txt"
      3 filename_Shakun2015_SST = "Shakun2015_SST.txt"
      4 url_Shakun2015_SST = "https://osf.io/kmy5w/download"
----> 5 SST = pd.read_table(pooch_load(url_Shakun2015_SST, filename_Shakun2015_SST))
      6 SST.set_index("Age", inplace=True)
      7 SST

Cell In[3], line 11, in pooch_load(filelocation, filename, processor)
      7 
      8     if os.path.exists(os.path.join(shared_location, filename)):
      9         file = os.path.join(shared_location, filename)
     10     else:
---> 11         file = pooch.retrieve(
     12             filelocation,
     13             known_hash=None,
     14             fname=os.path.join(user_temp_cache, filename),

File ~/micromamba/envs/climatematch/lib/python3.11/site-packages/pooch/core.py:242, in retrieve(url, known_hash, fname, path, processor, downloader, progressbar)
    239 if downloader is None:
    240     downloader = choose_downloader(url, progressbar=progressbar)
--> 242 stream_download(url, full_path, known_hash, downloader, pooch=None)
    244 if known_hash is None:
    245     get_logger().info(
    246         "SHA256 hash of downloaded file: %s\n"
    247         "Use this value as the 'known_hash' argument of 'pooch.retrieve'"
   (...)    250         file_hash(str(full_path)),
    251     )

File ~/micromamba/envs/climatematch/lib/python3.11/site-packages/pooch/core.py:823, in stream_download(url, fname, known_hash, downloader, pooch, retry_if_failed)
    819 try:
    820     # Stream the file to a temporary so that we can safely check its
    821     # hash before overwriting the original.
    822     with temporary_file(path=str(fname.parent)) as tmp:
--> 823         downloader(url, tmp, pooch)
    824         hash_matches(tmp, known_hash, strict=True, source=str(fname.name))
    825         shutil.move(tmp, str(fname))

File ~/micromamba/envs/climatematch/lib/python3.11/site-packages/pooch/downloaders.py:230, in HTTPDownloader.__call__(self, url, output_file, pooch, check_only)
    228     # pylint: enable=consider-using-with
    229 try:
--> 230     response = requests.get(url, timeout=timeout, **kwargs)
    231     response.raise_for_status()
    232     content = response.iter_content(chunk_size=self.chunk_size)

File ~/micromamba/envs/climatematch/lib/python3.11/site-packages/requests/api.py:87, in get(url, params, **kwargs)
     74 def get(
     75     url: _t.UriType, params: _t.ParamsType = None, **kwargs: Unpack[_t.GetKwargs]
     76 ) -> Response:
     77     r"""Sends a GET request.
     78 
     79     :param url: URL for the new :class:`Request` object.
   (...)     84     :rtype: requests.Response
     85     """
---> 87     return request("get", url, params=params, **kwargs)

File ~/micromamba/envs/climatematch/lib/python3.11/site-packages/requests/api.py:71, in request(method, url, **kwargs)
     67 # By using the 'with' statement we are sure the session is closed, thus we
     68 # avoid leaving sockets open which can trigger a ResourceWarning in some
     69 # cases, and look like a memory leak in others.
     70 with sessions.Session() as session:
---> 71     return session.request(method=method, url=url, **kwargs)

File ~/micromamba/envs/climatematch/lib/python3.11/site-packages/requests/sessions.py:651, in Session.request(self, method, url, params, data, headers, cookies, files, auth, timeout, allow_redirects, proxies, hooks, stream, verify, cert, json)
    646 send_kwargs = {
    647     "timeout": timeout,
    648     "allow_redirects": allow_redirects,
    649 }
    650 send_kwargs.update(settings)
--> 651 resp = self.send(prep, **send_kwargs)
    653 return resp

File ~/micromamba/envs/climatematch/lib/python3.11/site-packages/requests/sessions.py:805, in Session.send(self, request, **kwargs)
    802 if allow_redirects:
    803     # Redirect resolving generator.
    804     gen = self.resolve_redirects(r, request, **kwargs)
--> 805     history = [resp for resp in gen]
    806 else:
    807     history = []

File ~/micromamba/envs/climatematch/lib/python3.11/site-packages/requests/sessions.py:805, in <listcomp>(.0)
    802 if allow_redirects:
    803     # Redirect resolving generator.
    804     gen = self.resolve_redirects(r, request, **kwargs)
--> 805     history = [resp for resp in gen]
    806 else:
    807     history = []

File ~/micromamba/envs/climatematch/lib/python3.11/site-packages/requests/sessions.py:292, in SessionRedirectMixin.resolve_redirects(self, resp, req, stream, timeout, verify, cert, proxies, yield_requests, **adapter_kwargs)
    290     yield req  # type: ignore[misc]  # Internal use only, returns PreparedRequest
    291 else:
--> 292     resp = self.send(
    293         req,
    294         stream=stream,
    295         timeout=timeout,
    296         verify=verify,
    297         cert=cert,
    298         proxies=proxies,
    299         allow_redirects=False,
    300         **adapter_kwargs,
    301     )
    303     extract_cookies_to_jar(self.cookies, prepared_request, resp.raw)
    305     # extract redirect url, if any, for the next loop

File ~/micromamba/envs/climatematch/lib/python3.11/site-packages/requests/sessions.py:784, in Session.send(self, request, **kwargs)
    781 start = preferred_clock()
    783 # Send the request
--> 784 r = adapter.send(request, **kwargs)
    786 # Total elapsed time of the request (approximately)
    787 elapsed = preferred_clock() - start

File ~/micromamba/envs/climatematch/lib/python3.11/site-packages/requests/adapters.py:742, in HTTPAdapter.send(self, request, stream, timeout, verify, cert, proxies)
    740     raise SSLError(e, request=request)
    741 elif isinstance(e, ReadTimeoutError):
--> 742     raise ReadTimeout(e, request=request)
    743 elif isinstance(e, _InvalidHeader):
    744     raise InvalidHeader(e, request=request)

ReadTimeout: HTTPSConnectionPool(host='files.ca-1.osf.io', port=443): Read timed out. (read timeout=30)

# read CO2 dataantarctica2015co2composite_cleaned.txt
filename_antarctica2015co2composite_cleaned = "antarctica2015co2composite_cleaned.txt"
url_antarctica2015co2composite_cleaned = "https://osf.io/45fev/download"
CO2 = pd.read_table(
    pooch_load(
        url_antarctica2015co2composite_cleaned,
        filename_antarctica2015co2composite_cleaned,
    )
)
CO2.set_index("age_gas_calBP", inplace=True)
CO2

Downloading data from 'https://osf.io/45fev/download' to file '/tmp/antarctica2015co2composite_cleaned.txt'.

SHA256 hash of downloaded file: e86ec8dba4ca9a1e8404117bdd370fd830352d07164fe8093fff12330e371aa8
Use this value as the 'known_hash' argument of 'pooch.retrieve' to ensure that the file hasn't changed if it is downloaded again in the future.

	co2_ppm	co2_1s_ppm
age_gas_calBP
-51.03	368.02	0.06
-48.00	361.78	0.37
-46.28	359.65	0.10
-44.41	357.11	0.16
-43.08	353.95	0.04
...	...	...
803925.28	202.92	2.06
804009.87	207.50	0.92
804522.67	204.86	1.64
805132.44	202.23	0.69
805668.87	207.29	2.20

1901 rows × 2 columns

# plot
# set up two subplots in a grid of 2 rows and 1 column
# also make sure the two plots share the same x(time) axis
fig, axes = plt.subplots(2, 1, sharex=True)
# move the two subplots closer to each other
fig.subplots_adjust(hspace=-0.5)
axes[0].plot(SST.index, SST["SST stack"], color="C4")
axes[1].plot(CO2.index / 1000, CO2["co2_ppm"], color="C1")

# beautification
# since sharex=True in plt.subplots(), this sets the x axis limit for both panels
axes[1].set_xlim((0, 805))
# axis labels
axes[1].set_xlabel("Age (ka BP)")
axes[0].set_ylabel(r"Sea Surface Temperature" "\n" "Anomaly (°C)", color="C4")
axes[1].set_ylabel(r"CO${}_\mathrm{2}$ (ppm)", color="C1")

# despine makes the plots look cleaner
sns.despine(ax=axes[0], top=True, right=False, bottom=True, left=True)
sns.despine(ax=axes[1], top=True, right=True, bottom=False, left=False)
# clean up top panel x axis ticks
axes[0].xaxis.set_ticks_position("none")
# move top panel xlabel to the right side
axes[0].yaxis.set_label_position("right")
# the following code ensures the subplots don't overlap
for ax in axes:
    ax.set_zorder(10)
    ax.set_facecolor("none")
# color the axis
axes[0].spines["right"].set_color("C4")
axes[1].spines["left"].set_color("C1")
axes[0].tick_params(axis="y", colors="C4")
axes[1].tick_params(axis="y", colors="C1")

Now that we’ve taken a look at the two time series, let’s make a scatter plot between them and fit a linear regression model through the data.

# in this code block, we will make a scatter plot of CO2 and temperature
# and fit a linear regression model through the data


def age_model_interp(CO2_age, CO2, SST_age):
    """
    This helper function linearly interpolates CO2 data, which
    have a very high temporal resolution, to temperature data,
    which have a relatively low resolution
    """
    f = interpolate.interp1d(CO2_age, CO2)
    all_ages = f(SST_age)
    return all_ages


# interpolate CO2 data to SST age
CO2_interpolated = age_model_interp(CO2.index / 1000, CO2["co2_ppm"], SST.index)

# plot
# set up two subplots in a grid of 2 rows and 1 column
# also make sure the two plots share the same x(time) axis
fig, ax = plt.subplots(1, 1, sharex=True)

ax.scatter(CO2_interpolated, SST["SST stack"], color="gray")

# regression
X = CO2_interpolated
y = SST["SST stack"]
res = stats.linregress(X, y)  # ordinary least sqaure

x_fit = np.arange(180, 280)
# intercept
y_fit = x_fit * res.slope + res.intercept
ax.plot(x_fit, y_fit, color="k")

# beautification
# axis labels
ax.set_xlabel(r"CO${}_\mathrm{2}$ (ppm)")
ax.set_ylabel(r"Sea Surface Temperature" "\n" "Anomaly (°C)")
print(
    "Pearson (r^2) value: "
    + "{:.2f}".format(res.rvalue**2)
    + " \nwith a p-value of: "
    + "{:.2e}".format(res.pvalue)
)

Figure Making Through the Equity Lense

Click here for some information

Are the colors in your figure distinguishable for people with color-vision deficiencies?

More readings on this topic:

Contrast checker: https://www.color-blindness.com/coblis-color-blindness-simulator/

Coloring for color blindness: https://davidmathlogic.com/colorblind

Python-specific color palettes that are friendly to those with color-vision deficiency: https://seaborn.pydata.org/tutorial/color_palettes.html

Resources

Data from the following sources are used in this tutorial:

CO₂: Bereiter, B., Eggleston, S., Schmitt, J., Nehrbass-Ahles, C., Stocker, T.F., Fischer, H., Kipfstuhl, S., Chappellaz, J., 2015. Revision of the EPICA Dome C CO2 record from 800 to 600 kyr before present. Geophysical Research Letters 42, 542–549. https://doi.org/10.1002/2014GL061957

Temperature: Shakun, J.D., Lea, D.W., Lisiecki, L.E., Raymo, M.E., 2015. An 800-kyr record of global surface ocean δ18O and implications for ice volume-temperature coupling. Earth and Planetary Science Letters 426, 58–68. https://doi.org/10.1016/j.epsl.2015.05.042 (not Open Access)