Accessing Gridded Climate Data

• Sign up to the DEA Sandbox to run this notebook interactively from a browser

• Compatibility: Notebook currently compatible with the DEA Sandbox & NCI environments

• Products used: ERA5 and ANUClim (these datasets are external to the Digital Earth Australia platform)

Background & Description

This notebook demonstrates how to access and use Australian-specific climate data stored on the National Computing Infrastructure (NCI), and global climate reanalysis data (ERA5) hosted on Google’s cloud infrastructure. The content is structured into two parts.

1. The first section shows how to access an Australian climate dataset through the THREDDS platform. NCI’s THREDDS provides a wide range of climate, satellite, and other gridded environmental datasets, which can be accessed programmatically via xarray.

2. The second section focuses on the use of a pre-processed version of the ECMWF Reanalysis v5 (ERA5) Climate Fields. These datasets are also accessed using xarray.

Note: Additional details about the datasets are provided further in the notebook.

Load packages

Import Python packages that are used for the analysis.

%matplotlib inline

import os
import gcsfs
import requests
import datacube
import xarray as xr
import pandas as pd
import numpy as np
import contextily as ctx
import geopandas as gpd
from bs4 import BeautifulSoup
import matplotlib.pyplot as plt
from odc.geo.xr import assign_crs
from odc.geo.geom import Geometry, BoundingBox

import sys
sys.path.insert(1, '../Tools/')
from dea_tools.dask import create_local_dask_cluster

Analysis parameters

The following code cell is used to define variable needed later.

• vector_path: a file path to the polygon that defines our analysis region, the default example uses a boundary of Victoria.

# Boundary of victoria
vector_path = '../Supplementary_data/External_data_Climate/victoria_boundary.geojson'

Section 1: ANUClim climate data

In the first section of this notebook, we demonstrates retrieval of the Australian National University’s Climate data product (‘ANUClim’). These datasets are produced using topographically conditional spatial interpolation of Australia’s extensive network of weather stations.

The available variables include: * tavg: Average surface air temperature * rain: Total precipitation * vpd : Vapour pressure deficit * srad: Incoming solar radiation * evap: Pan evaporation

The entire archive can be viewed here

Data is available at a daily and monthly cadence, and has a 1 x 1 km spatial resolution, making it one of the highest resolution climate data products available for Australia. The time range of the data depends on the climate variable. For example, rainfall goes back to 1900, while air temperature begins in 1960. Data is updated annually.

Note: ANUClim is stored in the 'EPSG:4283' projection

---

Analysis parameters for ANUClim

• var: the climate variable to load.

• cadence: specifies whether to load monthly or daily data. Must be either ‘month’ or ‘day’.

• year_start & year_end: the year range of data you want to load, use integer values here.

# Which variable should we load?
var = 'rain' # other options: vpd, srad, rain, evap,  tavg

# Define the cadence of the data to load
cadence= 'month' # other options: day

# Define the time window (inclusive of end year) (use integers only)
year_start, year_end = 2020, 2020

Load our region of interest

# open the vector file
gdf = gpd.read_file(vector_path)

# convert to GDA94 to match ANUClim's CRS
gdf = gdf.to_crs('epsg:4283')

# define bounding box of the geometry
# we'll use this to select data before masking
bb = gdf.total_bounds

Load the climate data through THREDDS

Extract the links for loading the target year data from THREDDS using xarray.

# the top level folder on the NCI
thredds_data_link = 'https://thredds.nci.org.au/thredds/dodsC/gh70/ANUClimate/v2-0/stable/'
thredds_catalog_link=f'https://thredds.nci.org.au/thredds/catalog/gh70/ANUClimate/v2-0/stable/{cadence}/'

# list of years to load
years = [str(i) for i in range(year_start, year_end+1)]

i=0
pp = []

# loop through years and load
for y in years:
    print(" {:02}/{:02}\r".format(i + 1, len(years)), end="")

    # get list of all the file links on the THREDDS server
    ff = f'{thredds_catalog_link}{var}/{y}/catalog.html'
    soup = BeautifulSoup(requests.get(ff).content, "html.parser")

    # extract just the href links
    file_names = []
    for link in soup.select('a[href*=".html"]'): # loop over all <a> elements where the href attribute contains ".html"
        href = link["href"]
        if "dataset" in href: # add only file names of datasets to list
            file_names.append(href)

    # create a list of links for xarray
    list_of_links = [thredds_data_link+f.replace('catalog.html?dataset=gh70/', '') for f in file_names]

    # open with xarray and tidy up
    ds = xr.open_mfdataset(list_of_links)
    ds = assign_crs(ds, crs='epsg:4283') #GDA94
    ds = ds.drop_vars('crs')[var]
    ds.attrs['nodata'] = np.nan

    # select data from within the bounding box of the geometry
    # this is faster than masking, which we can do later
    ds = ds.sel(lat=slice(bb[3], bb[1]), lon=slice(bb[0],bb[2]))

    pp.append(ds)
    i += 1

# join the years together
climate = xr.concat(pp, dim='time').sortby('time')
climate.load() #bring into memory

# mask out all pixels outside of our ROI
# convert to a Geometry object for masking
geometry = gdf.iloc[0].geometry
geom = Geometry(geom=geometry, crs=gdf.crs) # use same CRS of the vector file
climate = climate.odc.mask(poly=geom)

 01/01

Plot annual data

Below we resample the climate data into an annual series. As we have loaded rainfall in the default example, we calculate an annual sum. This statistic may not make sense if you are loading a different variable, e.g. temperature.

climate_annual = climate.resample(time='1YE').sum()
climate_annual = climate_annual.odc.mask(poly=geom) # need to remask after a sum to remove zeros

/env/lib/python3.10/site-packages/odc/geo/_xr_interop.py:500: UserWarning: grid_mapping=crs is not pointing to valid coordinate
  warnings.warn(

fig, ax = plt.subplots(figsize=(6, 4), layout='constrained')

im = climate_annual.isel(time=0).plot(cmap='viridis', add_labels=False, robust=True, add_colorbar=False, ax=ax)
ax.set_title(climate_annual.isel(time=0).time.dt.year.item())

cb = fig.colorbar(im, ax=ax, shrink=0.75, orientation='vertical')
cb.ax.set_title(var + ' \n mm yr\u207B\u00B9', fontsize=8);

Section 2: ECMWF Reanalysis v5 (ERA5)

We can also load ERA5 reanalysis datasets using Google’s collection of Analysis-Ready, Cloud Optimized ERA5 data. ERA5 has a coarse spatial resolution (~25 km) compared with ANUClim, but provides hourly reconstructions of past climate and has many more climate fields at both surface and pressure levels.

A list of variables to load is available in the data summary table. Some commonly used long and short variable names are the following: * 2 metre temperature : 2t * total_precipitation : tp * 2m_dewpoint_temperature : 2d * evaporation : e * 10m_u_component_of_wind : 10u * 10m_v_component_of_wind : 10v * Surface short-wave (solar) radiation downwards : ssrd

Analysis parameters for ERA5

• era5_var: Provide a variable name to load

• time_range: length of time to load

era5_var = 'total_precipitation'
time_range = ('2020-01-01', '2020-12-31')

Open a dask client

This helps with ERA5 data as it has a larger volume of time steps to summarise.

create_local_dask_cluster()

Client

Client-421a5cbe-6110-11f0-8194-a257d99bdb0e

Connection method: Cluster object	Cluster type: distributed.LocalCluster
Dashboard: /user/chad.burton@ga.gov.au/proxy/8787/status

Launch dashboard in JupyterLab

Cluster Info

LocalCluster

ca9e844b

Dashboard: /user/chad.burton@ga.gov.au/proxy/8787/status	Workers: 1
Total threads: 2	Total memory: 15.00 GiB
Status: running	Using processes: True

Scheduler Info

Scheduler

Scheduler-05038069-2097-41e3-92b4-66463e14f10b

Comm: tcp://127.0.0.1:38757	Workers: 1
Dashboard: /user/chad.burton@ga.gov.au/proxy/8787/status	Total threads: 2
Started: Just now	Total memory: 15.00 GiB

Workers

Worker: 0

Comm: tcp://127.0.0.1:41071	Total threads: 2
Dashboard: /user/chad.burton@ga.gov.au/proxy/44255/status	Memory: 15.00 GiB
Nanny: tcp://127.0.0.1:44099
Local directory: /tmp/dask-scratch-space/worker-ptuthvpu

Define a function for loading ERA5

def load_era5(
    x=None,
    y=None,
    crs="EPSG:4326",
    time=None,
    bands=None,
    chunks=None,
    path="gs://gcp-public-data-arco-era5/ar/full_37-1h-0p25deg-chunk-1.zarr-v3",
):
    # Lazily load Zarr from Google Cloud Platform
    era5_ds = xr.open_zarr(path, chunks=chunks, storage_options=dict(token="anon"))

    # Select bands
    era5_ds = era5_ds[bands] if bands is not None else era5_ds

    # Clip to extent
    if x is not None:
        bbox = BoundingBox.from_xy(x, y, crs=crs).to_crs(crs)
        era5_ds = era5_ds.sel(
            longitude=slice(bbox.left, bbox.right),
            latitude=slice(bbox.top, bbox.bottom),
        )

    # Select time
    if time is not None:
        era5_ds = era5_ds.sel(time=slice(time[0], time[-1]))

    return era5_ds.odc.assign_crs(crs)

Convert the polygon back to EPSG:4326

# convert our vector back to EPSG 4326 for ERA5
gdf = gdf.to_crs('epsg:4326')

# define bounding box of the geometry
# we'll use this to select data before masking
bb = gdf.total_bounds

Load total precipitation

era5_climate = load_era5(
    x=(bb[0],bb[2]),
    y=(bb[1],bb[3]),
    crs="EPSG:4326",
    time=time_range,
    bands=[era5_var],
    chunks=dict(x=10, y=10, time=2000), # really small chunks
    path="gs://gcp-public-data-arco-era5/ar/full_37-1h-0p25deg-chunk-1.zarr-v3",
)

era5_climate = era5_climate[era5_var]
era5_climate

<xarray.DataArray 'total_precipitation' (time: 8784, latitude: 21, longitude: 36)> Size: 27MB
dask.array<getitem, shape=(8784, 21, 36), dtype=float32, chunksize=(2000, 21, 36), chunktype=numpy.ndarray>
Coordinates:
  * latitude     (latitude) float32 84B -34.0 -34.25 -34.5 ... -38.75 -39.0
  * longitude    (longitude) float32 144B 141.0 141.2 141.5 ... 149.5 149.8
  * time         (time) datetime64[ns] 70kB 2020-01-01 ... 2020-12-31T23:00:00
    spatial_ref  int32 4B 4326
Attributes:
    long_name:   Total precipitation
    short_name:  tp
    units:       m

xarray.DataArray

'total_precipitation'

• time: 8784
• latitude: 21
• longitude: 36

• dask.array<chunksize=(104, 21, 36), meta=np.ndarray>

  Array Chunk Bytes 25.33 MiB 5.77 MiB Shape (8784, 21, 36) (2000, 21, 36) Dask graph 6 chunks in 4 graph layers Data type float32 numpy.ndarray

• Coordinates: (4)
  • latitude
    (latitude)
    float32
    -34.0 -34.25 -34.5 ... -38.75 -39.0

    long_name :

    latitude

    units :

    degrees_north

    array([-34.  , -34.25, -34.5 , -34.75, -35.  , -35.25, -35.5 , -35.75, -36.  ,
           -36.25, -36.5 , -36.75, -37.  , -37.25, -37.5 , -37.75, -38.  , -38.25,
           -38.5 , -38.75, -39.  ], dtype=float32)

  • longitude
    (longitude)
    float32
    141.0 141.2 141.5 ... 149.5 149.8

    long_name :

    longitude

    units :

    degrees_east

    array([141.  , 141.25, 141.5 , 141.75, 142.  , 142.25, 142.5 , 142.75, 143.  ,
           143.25, 143.5 , 143.75, 144.  , 144.25, 144.5 , 144.75, 145.  , 145.25,
           145.5 , 145.75, 146.  , 146.25, 146.5 , 146.75, 147.  , 147.25, 147.5 ,
           147.75, 148.  , 148.25, 148.5 , 148.75, 149.  , 149.25, 149.5 , 149.75],
          dtype=float32)

  • time
    (time)
    datetime64[ns]
    2020-01-01 ... 2020-12-31T23:00:00

    array(['2020-01-01T00:00:00.000000000', '2020-01-01T01:00:00.000000000',
           '2020-01-01T02:00:00.000000000', ..., '2020-12-31T21:00:00.000000000',
           '2020-12-31T22:00:00.000000000', '2020-12-31T23:00:00.000000000'],
          dtype='datetime64[ns]')

  • spatial_ref
    ()
    int32
    4326

    spatial_ref :

    GEOGCRS["WGS 84",ENSEMBLE["World Geodetic System 1984 ensemble",MEMBER["World Geodetic System 1984 (Transit)"],MEMBER["World Geodetic System 1984 (G730)"],MEMBER["World Geodetic System 1984 (G873)"],MEMBER["World Geodetic System 1984 (G1150)"],MEMBER["World Geodetic System 1984 (G1674)"],MEMBER["World Geodetic System 1984 (G1762)"],MEMBER["World Geodetic System 1984 (G2139)"],MEMBER["World Geodetic System 1984 (G2296)"],ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ENSEMBLEACCURACY[2.0]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],CS[ellipsoidal,2],AXIS["geodetic latitude (Lat)",north,ORDER[1],ANGLEUNIT["degree",0.0174532925199433]],AXIS["geodetic longitude (Lon)",east,ORDER[2],ANGLEUNIT["degree",0.0174532925199433]],USAGE[SCOPE["Horizontal component of 3D system."],AREA["World."],BBOX[-90,-180,90,180]],ID["EPSG",4326]]

    crs_wkt :

    GEOGCRS["WGS 84",ENSEMBLE["World Geodetic System 1984 ensemble",MEMBER["World Geodetic System 1984 (Transit)"],MEMBER["World Geodetic System 1984 (G730)"],MEMBER["World Geodetic System 1984 (G873)"],MEMBER["World Geodetic System 1984 (G1150)"],MEMBER["World Geodetic System 1984 (G1674)"],MEMBER["World Geodetic System 1984 (G1762)"],MEMBER["World Geodetic System 1984 (G2139)"],MEMBER["World Geodetic System 1984 (G2296)"],ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ENSEMBLEACCURACY[2.0]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],CS[ellipsoidal,2],AXIS["geodetic latitude (Lat)",north,ORDER[1],ANGLEUNIT["degree",0.0174532925199433]],AXIS["geodetic longitude (Lon)",east,ORDER[2],ANGLEUNIT["degree",0.0174532925199433]],USAGE[SCOPE["Horizontal component of 3D system."],AREA["World."],BBOX[-90,-180,90,180]],ID["EPSG",4326]]

    semi_major_axis :

    6378137.0

    semi_minor_axis :

    6356752.314245179

    inverse_flattening :

    298.257223563

    reference_ellipsoid_name :

    WGS 84

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    WGS 84

    horizontal_datum_name :

    World Geodetic System 1984 ensemble

    grid_mapping_name :

    latitude_longitude

    array(4326, dtype=int32)

• Indexes: (3)
  • latitude
    PandasIndex

    PandasIndex(Index([ -34.0, -34.25,  -34.5, -34.75,  -35.0, -35.25,  -35.5, -35.75,  -36.0,
           -36.25,  -36.5, -36.75,  -37.0, -37.25,  -37.5, -37.75,  -38.0, -38.25,
            -38.5, -38.75,  -39.0],
          dtype='float32', name='latitude'))

  • longitude
    PandasIndex

    PandasIndex(Index([ 141.0, 141.25,  141.5, 141.75,  142.0, 142.25,  142.5, 142.75,  143.0,
           143.25,  143.5, 143.75,  144.0, 144.25,  144.5, 144.75,  145.0, 145.25,
            145.5, 145.75,  146.0, 146.25,  146.5, 146.75,  147.0, 147.25,  147.5,
           147.75,  148.0, 148.25,  148.5, 148.75,  149.0, 149.25,  149.5, 149.75],
          dtype='float32', name='longitude'))

  • time
    PandasIndex

    PandasIndex(DatetimeIndex(['2020-01-01 00:00:00', '2020-01-01 01:00:00',
                   '2020-01-01 02:00:00', '2020-01-01 03:00:00',
                   '2020-01-01 04:00:00', '2020-01-01 05:00:00',
                   '2020-01-01 06:00:00', '2020-01-01 07:00:00',
                   '2020-01-01 08:00:00', '2020-01-01 09:00:00',
                   ...
                   '2020-12-31 14:00:00', '2020-12-31 15:00:00',
                   '2020-12-31 16:00:00', '2020-12-31 17:00:00',
                   '2020-12-31 18:00:00', '2020-12-31 19:00:00',
                   '2020-12-31 20:00:00', '2020-12-31 21:00:00',
                   '2020-12-31 22:00:00', '2020-12-31 23:00:00'],
                  dtype='datetime64[ns]', name='time', length=8784, freq=None))

• Attributes: (3)

  long_name :

  Total precipitation

  short_name :

  tp

  units :

  m

Resample to annual

As ERA5 has an hourly time step, we need to use dask to resample into annual time fields as there is a lot of data even though the spatial extent of the subset of data is small. This can take a few minutes on the default sandbox.

era5_climate_annual = era5_climate.resample(time='1YE').sum().load()

Mask with the polygon

geometry = gdf.iloc[0].geometry
geom = Geometry(geom=geometry, crs=gdf.crs)
era5_climate_annual = era5_climate_annual.odc.mask(poly=geom) # need to remask after a sum to remove zeros

Plot an annual summary data

fig, ax = plt.subplots(figsize=(5, 4), layout='constrained')

im = (era5_climate_annual*1000).isel(time=0).plot(cmap='viridis', add_labels=False, robust=True, add_colorbar=False, ax=ax);  # multiply by 1000 to convert into mm
ax.set_title(era5_climate_annual.isel(time=0).time.dt.year.item())

cb = fig.colorbar(im, ax=ax, shrink=0.75, orientation='vertical')
cb.ax.set_title(var + ' \n mm yr\u207B\u00B9', fontsize=8);

---

Additional information

License: The code in this notebook is licensed under the Apache License, Version 2.0. Digital Earth Australia data is licensed under the Creative Commons by Attribution 4.0 license.

Contact: If you need assistance, please post a question on the Open Data Cube Discord chat or on the GIS Stack Exchange using the open-data-cube tag (you can view previously asked questions here). If you would like to report an issue with this notebook, you can file one on GitHub.

Last modified: April 2025

Compatible datacube version:

print(datacube.__version__)

1.8.19

Tags

Tags: NCI compatible, sandbox compatible, :index:`precipitation, air temperature, external data, climate data, downloading data, no_testing