Open In Colab   Open in Kaggle

Tutorial 2: Selection, Interpolation and Slicing#

Week 1, Day 1, Climate System Overview

Content creators: Sloane Garelick, Julia Kent

Content reviewers: Katrina Dobson, Younkap Nina Duplex, Danika Gupta, Maria Gonzalez, Will Gregory, Nahid Hasan, Paul Heubel, Sherry Mi, Beatriz Cosenza Muralles, Jenna Pearson, Agustina Pesce, Chi Zhang, Ohad Zivan

Content editors: Paul Heubel, 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

project pythia#

Pythia credit: Rose, B. E. J., Kent, J., Tyle, K., Clyne, J., Banihirwe, A., Camron, D., May, R., Grover, M., Ford, R. R., Paul, K., Morley, J., Eroglu, O., Kailyn, L., & Zacharias, A. (2023). Pythia Foundations (Version v2023.05.01) https://zenodo.org/record/8065851

CMIP.png#

Tutorial Objectives#

Estimated timing of tutorial: 25 minutes

In the previous tutorial, we learned how to use Xarray to create DataArray and Dataset objects. Global climate datasets can be very large with multiple variables, and DataArrays and Datasets are very useful tools for organizing, comparing, and interpreting such data. However, sometimes we are not interested in examining a global dataset but wish to examine a specific time or location. For example, we might want to look at climate variables in a particular region of Earth, and potentially compare that to another region. In order to carry out such analyses, it’s useful to be able to extract and compare subsets of data from a global dataset.

In this tutorial, you will explore multiple computational tools in Xarray that allow you to select data from a specific spatial and temporal range. In particular, you will practice using:

  • .sel(): select data based on coordinate values or date

  • .interp(): interpolate to any latitude/longitude location to extract data

  • slice(): to select a range (or slice) along one or more coordinates, we can pass a Python slice object to .sel()

Setup#

# imports
from datetime import timedelta
import numpy as np
import pandas as pd
import xarray as xr
import matplotlib.pyplot as plt

Install and import feedback gadget#

Hide code cell source 
# @title Install and import feedback gadget

!pip3 install vibecheck datatops --quiet

from vibecheck import DatatopsContentReviewContainer
def content_review(notebook_section: str):
    return DatatopsContentReviewContainer(
        "",  # No text prompt
        notebook_section,
        {
            "url": "https://pmyvdlilci.execute-api.us-east-1.amazonaws.com/klab",
            "name": "comptools_4clim",
            "user_key": "l5jpxuee",
        },
    ).render()


feedback_prefix = "W1D1_T2"

[notice] A new release of pip is available: 24.2 -> 24.3.1
[notice] To update, run: pip install --upgrade pip

Figure Settings#

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")

Video 1: Solar Radiation and Earth’s Energy Budget#

Submit your feedback#

Hide code cell source 
# @title Submit your feedback
content_review(f"{feedback_prefix}_Solar_Radiation_Video")

If you want to download the slides: https://osf.io/download/gh5us/

Submit your feedback#

Hide code cell source 
# @title Submit your feedback
content_review(f"{feedback_prefix}_Solar_Radiation_Slides")

Section 1: Subsetting and Selection by Coordinate Values#

Since Xarray allows us to label coordinates, you can select data based on coordinate names and values, rather than array indices. We’ll explore this briefly here. To explore these Xarray tools, first recreate the synthetic temperature and pressure DataArrays you generated in the previous Tutorial 1, and combine these two DataArrays into a Dataset.

# temperature data
rand_data = 283 + 5 * np.random.randn(5, 3, 4)
times_index = pd.date_range("2018-01-01", periods=5)
lons = np.linspace(-120, -60, 4)
lats = np.linspace(25, 55, 3)
temperature = xr.DataArray(
    rand_data, coords=[times_index, lats, lons], dims=["time", "lat", "lon"]
)
temperature.attrs["units"] = "Kelvin"
temperature.attrs["standard_name"] = "air_temperature"

# pressure data
pressure_data = 1000.0 + 5 * np.random.randn(5, 3, 4)
pressure = xr.DataArray(
    pressure_data, coords=[times_index, lats, lons], dims=["time", "lat", "lon"]
)
pressure.attrs["units"] = "hPa"
pressure.attrs["standard_name"] = "air_pressure"

# combinate temperature and pressure DataArrays into a Dataset called 'ds'
ds = xr.Dataset(data_vars={"Temperature": temperature, "Pressure": pressure})
ds

<xarray.Dataset> Size: 1kB
Dimensions:      (time: 5, lat: 3, lon: 4)
Coordinates:
  * time         (time) datetime64[ns] 40B 2018-01-01 2018-01-02 ... 2018-01-05
  * lat          (lat) float64 24B 25.0 40.0 55.0
  * lon          (lon) float64 32B -120.0 -100.0 -80.0 -60.0
Data variables:
    Temperature  (time, lat, lon) float64 480B 275.9 281.2 278.3 ... 286.4 278.6
    Pressure     (time, lat, lon) float64 480B 996.1 995.0 ... 1.009e+03 997.6

To refresh your memory from the previous tutorial, take a look at the DataArrays you created for temperature and pressure by clicking on those variables in the dataset above.

Section 1.1: NumPy-like Selection#

Suppose you want to extract all the spatial data for one single date: January 2, 2018. It’s possible to achieve that with NumPy-like index selection:

indexed_selection = temperature[
    1, :, :
]  # index 1 along axis 0 is the time slice we want...
indexed_selection

<xarray.DataArray (lat: 3, lon: 4)> Size: 96B
array([[282.05710807, 278.58385409, 280.71299931, 295.5668969 ],
       [286.86994484, 284.72952072, 285.53542154, 281.07689852],
       [273.93279516, 292.27008303, 286.84415105, 278.79981171]])
Coordinates:
    time     datetime64[ns] 8B 2018-01-02
  * lat      (lat) float64 24B 25.0 40.0 55.0
  * lon      (lon) float64 32B -120.0 -100.0 -80.0 -60.0
Attributes:
    units:          Kelvin
    standard_name:  air_temperature

However, notice that this requires us (the user) to have detailed knowledge of the order of the axes and the meaning of the indices along those axes. By having named coordinates in Xarray, we can avoid this issue.

Section 1.2: .sel()#

Rather than using a NumPy-like index selection, in Xarray, we can instead select data based on coordinate values using the .sel() method, which takes one or more named coordinate(s) as a keyword argument:

named_selection = temperature.sel(time="2018-01-02")
named_selection

<xarray.DataArray (lat: 3, lon: 4)> Size: 96B
array([[282.05710807, 278.58385409, 280.71299931, 295.5668969 ],
       [286.86994484, 284.72952072, 285.53542154, 281.07689852],
       [273.93279516, 292.27008303, 286.84415105, 278.79981171]])
Coordinates:
    time     datetime64[ns] 8B 2018-01-02
  * lat      (lat) float64 24B 25.0 40.0 55.0
  * lon      (lon) float64 32B -120.0 -100.0 -80.0 -60.0
Attributes:
    units:          Kelvin
    standard_name:  air_temperature

We got the same result as when we used the NumPy-like index selection, but

  • we didn’t have to know anything about how the array was created or stored

  • our code is agnostic about how many dimensions we are dealing with

  • the intended meaning of our code is much clearer!

By using the .sel() method in Xarray, we can easily isolate data from a specific time. You can also isolate data from a specific coordinate.

Coding Exercises 1.2#

  1. Write a line of code to select the temperature data from the grid point with the coordinates (25, -120).

coordinate_selection = ...
coordinate_selection

Ellipsis

Click for solution

Submit your feedback#

Hide code cell source 
# @title Submit your feedback
content_review(f"{feedback_prefix}_Coding_Exercise_1_2")

Section 1.3: Approximate Selection and Interpolation#

The spatial and temporal resolution of climate data often differs between datasets or a dataset may be incomplete. Therefore, with time and space data, we frequently want to sample “near” the coordinate points in our dataset. For example, we may want to analyze data from a specific coordinate or a specific time, but may not have a value from that specific location or date. In that case, we would want to use the data from the closest coordinate or time step. Here are a few simple ways to achieve that.

Section 1.3.1: Nearest-neighbor Sampling#

Suppose we want to know the temperature from 2018-01-07. However, the last day on our time axis is 2018-01-05. We can therefore sample within two days of our desired date of 2018-01-07. We can do this using the .sel() method we used earlier but with the added flexibility of performing nearest neighbor sampling and specifying an optional tolerance. This is called an inexact lookup because we are not searching for a perfect match, although there may be one. Here the tolerance is the maximum distance away from our desired point Xarray will search for the nearest neighbor.

temperature.sel(time="2018-01-07", method="nearest", tolerance=timedelta(days=2))

<xarray.DataArray (lat: 3, lon: 4)> Size: 96B
array([[279.61714087, 282.09886047, 274.87710805, 281.00120504],
       [275.91970126, 278.80432915, 287.39513783, 289.83979555],
       [281.8850557 , 278.43863234, 286.35977048, 278.59065652]])
Coordinates:
    time     datetime64[ns] 8B 2018-01-05
  * lat      (lat) float64 24B 25.0 40.0 55.0
  * lon      (lon) float64 32B -120.0 -100.0 -80.0 -60.0
Attributes:
    units:          Kelvin
    standard_name:  air_temperature

Notice that the resulting data is from the date 2018-01-05.

Section 1.3.2: Interpolation#

The latitude values of our dataset are 25ºN, 40ºN, 55ºN, and the longitude values are 120ºW, 100ºW, 80ºW, 60ºW. But suppose we want to extract a timeseries for Boulder, Colorado, USA (40°N, 105°W). Since lon=-105 is not a point on our longitude axis, this requires interpolation between data points.

We can do this using the .interp() method (see the docs here), which works similarly to .sel(). Using .interp(), we can interpolate to any latitude/longitude location using an interpolation method of our choice. In the example below, you will linearly interpolate between known points.

temperature.interp(lon=-105, lat=40, method="linear")

<xarray.DataArray (time: 5)> Size: 40B
array([283.8375919 , 285.26462675, 279.6042351 , 281.53411262,
       278.08317218])
Coordinates:
  * time     (time) datetime64[ns] 40B 2018-01-01 2018-01-02 ... 2018-01-05
    lon      int64 8B -105
    lat      int64 8B 40
Attributes:
    units:          Kelvin
    standard_name:  air_temperature

In this case, we specified a linear interpolation method, yet one can choose other methods as well (e.g., nearest, cubic, quadratic). Note that the temperature values we extracted in the code cell above are not actual values in the dataset, but are instead calculated based on linear interpolations between values that are in the dataset.

Section 1.4: Slicing Along Coordinates#

Frequently we want to select a range (or slice) along one or more coordinate(s). For example, you may wish to only assess average annual temperatures in equatorial regions. We can achieve this by passing a Python slice object to .sel(). The calling sequence for slice always looks like slice(start, stop[, step]), where step is optional. In this case, let’s only look at values between 110ºW-70ºW and 25ºN-40ºN:

temperature.sel(
    time=slice("2018-01-01", "2018-01-03"), lon=slice(-110, -70), lat=slice(25, 45)
)

<xarray.DataArray (time: 3, lat: 2, lon: 2)> Size: 96B
array([[[281.24503913, 278.33921522],
        [282.86476396, 277.8111101 ]],

       [[278.58385409, 280.71299931],
        [284.72952072, 285.53542154]],

       [[290.84626352, 286.80501845],
        [278.22233254, 286.28767412]]])
Coordinates:
  * time     (time) datetime64[ns] 24B 2018-01-01 2018-01-02 2018-01-03
  * lat      (lat) float64 16B 25.0 40.0
  * lon      (lon) float64 16B -100.0 -80.0
Attributes:
    units:          Kelvin
    standard_name:  air_temperature

Section 1.5: One More Selection Method: .loc#

All of these operations can also be done within square brackets on the .loc attribute of the DataArray:

temperature.loc['2018-01-02']

<xarray.DataArray (lat: 3, lon: 4)> Size: 96B
array([[282.05710807, 278.58385409, 280.71299931, 295.5668969 ],
       [286.86994484, 284.72952072, 285.53542154, 281.07689852],
       [273.93279516, 292.27008303, 286.84415105, 278.79981171]])
Coordinates:
    time     datetime64[ns] 8B 2018-01-02
  * lat      (lat) float64 24B 25.0 40.0 55.0
  * lon      (lon) float64 32B -120.0 -100.0 -80.0 -60.0
Attributes:
    units:          Kelvin
    standard_name:  air_temperature

This is sort of in between the NumPy-style selection

temp[1,:,:]

and the fully label-based selection using .sel()

With .loc, we make use of the coordinate values, but lose the ability to specify the names of the various dimensions. Instead, the slicing must be done in the correct order:

temperature.loc['2018-01-01':'2018-01-03', 25:45, -110:-70]

<xarray.DataArray (time: 3, lat: 2, lon: 2)> Size: 96B
array([[[281.24503913, 278.33921522],
        [282.86476396, 277.8111101 ]],

       [[278.58385409, 280.71299931],
        [284.72952072, 285.53542154]],

       [[290.84626352, 286.80501845],
        [278.22233254, 286.28767412]]])
Coordinates:
  * time     (time) datetime64[ns] 24B 2018-01-01 2018-01-02 2018-01-03
  * lat      (lat) float64 16B 25.0 40.0
  * lon      (lon) float64 16B -100.0 -80.0
Attributes:
    units:          Kelvin
    standard_name:  air_temperature

One advantage of using .loc is that we can use NumPy-style slice notation like 25:45, rather than the more verbose slice(25,45). But of course that also works:

temperature.loc["2018-01-01":"2018-01-03", slice(25, 45), -110:-70]

<xarray.DataArray (time: 3, lat: 2, lon: 2)> Size: 96B
array([[[281.24503913, 278.33921522],
        [282.86476396, 277.8111101 ]],

       [[278.58385409, 280.71299931],
        [284.72952072, 285.53542154]],

       [[290.84626352, 286.80501845],
        [278.22233254, 286.28767412]]])
Coordinates:
  * time     (time) datetime64[ns] 24B 2018-01-01 2018-01-02 2018-01-03
  * lat      (lat) float64 16B 25.0 40.0
  * lon      (lon) float64 16B -100.0 -80.0
Attributes:
    units:          Kelvin
    standard_name:  air_temperature

What doesn’t work is passing the slices in a different order to the dimensions of the dataset:

# This will generate an error
# temperature.loc[-110:-70, 25:45,'2018-01-01':'2018-01-03']

Summary#

In this tutorial, we have explored the practical use of .sel(), .interp(), .loc, and slicing techniques to extract data from specific spatial and temporal ranges. These methods are valuable when we are interested in only certain pieces of large datasets.

Resources#

Code and data for this tutorial is based on existing content from Project Pythia.