Calculating Climatology and Departures from Time Series Data#

Author: Tom Vo

Date: 05/27/22

Last Updated: 2/27/23

Related APIs:

The data used in this example can be found through the Earth System Grid Federation (ESGF) search portal.

Overview#

Suppose we have two netCDF4 files for air temperature data (tas).

File 1: Monthly frequency from 1850-01-16 to 2014-12-16
- We want to calculate the annual and seasonal cycle climatologies and departures using this file.
File 2: Hourly frequency from 2010-01-01 to 2015-01-01 (subset).
- We want to calculate the daily cycle climatologies and departures using this file.

Notebook Setup#

Create an Anaconda environment for this notebook using the command below, then select the kernel in Jupyter.

conda create -n xcdat_notebook -c conda-forge python xarray netcdf4 xcdat xesmf matplotlib nc-time-axis jupyter

xesmf is required for horizontal regridding with xESMF
matplotlib is an optional dependency required for plotting with xarray
nc-time-axis is an optional dependency required for matplotlib to plot cftime coordinates

[1]:

%matplotlib inline

import matplotlib.pyplot as plt
import pandas as pd
import xcdat

1. Open Sample Datasets#

We are using xarray’s OPeNDAP support to read a netCDF4 dataset files directly from their source. The data is not loaded over the network until we perform operations on it (e.g., temperature unit adjustment).

More information on the xarray’s OPeNDAP support can be found here.

File 2: Hourly Frequency#

The size of this file is approximately 1.17 GB, so we will be chunking our request using Dask to avoid hitting the OPeNDAP file size request limit for this ESGF node.

[3]:

filepath2 = "http://esgf.nci.org.au/thredds/dodsC/master/CMIP6/CMIP/CSIRO/ACCESS-ESM1-5/historical/r10i1p1f1/3hr/tas/gn/v20200605/tas_3hr_ACCESS-ESM1-5_historical_r10i1p1f1_gn_201001010300-201501010000.nc"
ds_hourly = xcdat.open_dataset(filepath2, chunks={"time": "auto"}, add_bounds=["T"])

# Unit adjust (-273.15, K to C)
ds_hourly["tas"] = ds_hourly.tas - 273.15
ds_hourly

[3]:

<xarray.Dataset> Size: 2GB
Dimensions:    (lat: 145, bnds: 2, lon: 192, time: 14608)
Coordinates:
  * lat        (lat) float64 1kB -90.0 -88.75 -87.5 -86.25 ... 87.5 88.75 90.0
  * lon        (lon) float64 2kB 0.0 1.875 3.75 5.625 ... 354.4 356.2 358.1
    height     float64 8B ...
  * time       (time) object 117kB 2010-01-01 03:00:00 ... 2015-01-01 00:00:00
Dimensions without coordinates: bnds
Data variables:
    lat_bnds   (lat, bnds) float64 2kB dask.array<chunksize=(145, 2), meta=np.ndarray>
    lon_bnds   (lon, bnds) float64 3kB dask.array<chunksize=(192, 2), meta=np.ndarray>
    tas        (time, lat, lon) float32 2GB dask.array<chunksize=(1205, 145, 192), meta=np.ndarray>
    time_bnds  (time, bnds) object 234kB 2010-01-01 03:00:00 ... 2015-01-01 0...
Attributes: (12/48)
    Conventions:                     CF-1.7 CMIP-6.2
    activity_id:                     CMIP
    branch_method:                   standard
    branch_time_in_child:            0.0
    branch_time_in_parent:           87658.0
    creation_date:                   2020-06-05T04:54:56Z
    ...                              ...
    variant_label:                   r10i1p1f1
    version:                         v20200605
    license:                         CMIP6 model data produced by CSIRO is li...
    cmor_version:                    3.4.0
    tracking_id:                     hdl:21.14100/b79e6a05-c482-46cf-b3b8-83b...
    DODS_EXTRA.Unlimited_Dimension:  time

2. Calculate Climatology#

In this example, we will be calculating the weighted climatology of the tas variable for its seasonal, annual, and daily cycles.

Helpful knowledge:

Masked (missing) data is automatically handled.
- The weight of masked (missing) data is excluded when averages are calculated. This is the same as giving them a weight of 0.
If desired, use the reference_period argument to calculate a climatology based on a climatological reference period (a subset of the entire time series). If no value is provided, the climatological reference period will be the full period covered by the dataset.

Daily Climatology#

Groups time coordinates by month and day.

Leap days (if present) are dropped if the CF calendar type is "gregorian", "proleptic_gregorian", or "standard".

[13]:

daily_climo = ds_hourly.temporal.climatology("tas", freq="day", weighted=True)

[14]:

daily_climo.tas

[14]:

<xarray.DataArray 'tas' (time: 365, lat: 145, lon: 192)> Size: 81MB
dask.array<truediv, shape=(365, 145, 192), dtype=float64, chunksize=(1, 145, 192), chunktype=numpy.ndarray>
Coordinates:
  * lat      (lat) float64 1kB -90.0 -88.75 -87.5 -86.25 ... 87.5 88.75 90.0
  * lon      (lon) float64 2kB 0.0 1.875 3.75 5.625 ... 352.5 354.4 356.2 358.1
    height   float64 8B ...
  * time     (time) object 3kB 0001-01-01 00:00:00 ... 0001-12-31 00:00:00
Attributes:
    operation:  temporal_avg
    mode:       climatology
    freq:       day
    weighted:   True

3. Calculate Departures (Anomalies)#

In this example, we will be calculating the weighted departures of the tas variable for its seasonal, annual, and daily cycles.

Helpful knowledge:

What are anomalies?
- In climatology, “anomalies” refer to the difference between the value during a given time interval (e.g., the January average surface air temperature) and the long-term average value for that time interval (e.g., the average surface temperature over the last 30 Januaries).
How is the climatology calculated?
- In the departures API, the reference climatology is calculated internally so there is no need to pass one to this method.
- You can still calculate the reference climatology using the climatology API.
- If desired, use the reference_period argument to calculate anomalies relative to a climatological reference period (a subset of the entire time series). If no value is provided, the climatological reference period will be the full period covered by the dataset.
Masked (missing) data is automatically handled.
- The weight of masked (missing) data is excluded when averages are calculated. This is the same as giving them a weight of 0.

Daily Anomalies#

Leap days (if present) are dropped if the CF calendar type is "gregorian", "proleptic_gregorian", or "standard".

[21]:

daily_departures = ds_hourly.temporal.departures("tas", freq="day", weighted=True)

[22]:

daily_departures.tas

[22]:

<xarray.DataArray 'tas' (time: 1826, lat: 145, lon: 192)> Size: 407MB
dask.array<sub, shape=(1826, 145, 192), dtype=float64, chunksize=(1, 145, 192), chunktype=numpy.ndarray>
Coordinates:
  * lat      (lat) float64 1kB -90.0 -88.75 -87.5 -86.25 ... 87.5 88.75 90.0
  * lon      (lon) float64 2kB 0.0 1.875 3.75 5.625 ... 352.5 354.4 356.2 358.1
    height   float64 8B ...
  * time     (time) object 15kB 2010-01-01 00:00:00 ... 2015-01-01 00:00:00
Attributes:
    operation:  temporal_avg
    mode:       departures
    freq:       day
    weighted:   True

Calculating Climatology and Departures from Time Series Data

Contents

Calculating Climatology and Departures from Time Series Data#

Overview#

Notebook Setup#

1. Open Sample Datasets#

File 1: Monthly Frequency#

File 2: Hourly Frequency#

2. Calculate Climatology#

Seasonal Climatology#

Custom Seasonal Climatology#

Annual Climatology#

Daily Climatology#

3. Calculate Departures (Anomalies)#

Seasonal Anomalies#

Custom Seasonal Anomalies#

Annual Anomalies#

Daily Anomalies#