Calculate the monthly climatological precipitation from ERA5

This notebook demonstrates how to retrieve, process, and plot the data displayed in the ERA Explorer

==> For further training materials, consider the Copernicus Climate Change Service (C3S) data tutorials

In this example we will be using earthkit to retrieve the data and standard Python packages to process and plot it.

import pandas as pd
import matplotlib.pyplot as plt
import calendar
import earthkit.data

If you want to run this code yourself, make sure to get a CDS API key from the CDS website.

Following the guidance there, you’ll then need to set up your local .cdsapirc to be able to authenticate with the CDS and download data. Then, you will be able to customise the date ranges, variables, and much more!

This is the latitude/longitude that we will extract the nearest grid cell from. You can change these to plot different locations. You can see the currently used ones in the ERA Explorer URL. In this example we are plotting data at the gridpoint closest to Brussels, Belgium.

lat = 50.86
lng = 4.35

This is the variable we are getting from the ERA dataset, and the time period

variable = "total_precipitation"
date_range = ["1991-01-01", "2020-12-31"]

For convenience, we make a function to handle the data retrieval:

def retrieve_data(variable, date_range, lat, lng):
    # Define the dataset and request parameters
    dataset = "reanalysis-era5-single-levels-timeseries"
    request = {
        "variable": [
        variable,  # Variable to retrieve
        ],
        "date": date_range,  # Date range for the data
        "location": {"longitude": lng, "latitude": lat},  # Location coordinates
        "data_format": "netcdf"  # Format of the retrieved data
    }

    # Use "earthkit" to retrieve the data
    ekds = earthkit.data.from_source(
        "cds", dataset, request
    ).to_xarray()

    return ekds

Get the data. This will download a NetCDF file

data = retrieve_data(variable, date_range, lat, lng)

And let’s make some bespoke functions to process the data. Each one is described with a “””docstring”””

# Make a function to compute the monthly precipitation climatology
def precipMonthlyClimatology():
    """
    Calculate the monthly climatology of precipitation.

    This function reads precipitation data from a NetCDF file,
    converts the time coordinate to a pandas datetime index,
    and then resamples the data to calculate the monthly 
    climatology. The resulting climatology is returned in millimeters.

    Returns:
        pandas.DataFrame: A DataFrame containing the monthly climatology
        of precipitation in millimeters, indexed by month.
    """


    data_tp_pt = data.tp

    # Convert the time coordinate to a pandas datetime index
    time_index = pd.to_datetime(data_tp_pt.valid_time.values)

    # Create a DataFrame for easier manipulation
    df = pd.DataFrame(data_tp_pt.values, index=time_index, columns=['tp'])

    df_monthly = df.resample('ME').sum()
    df_monthly['month'] = df_monthly.index.month
    monthly_climatology = df_monthly.groupby('month').mean() * 1000

    # Get the actual lat/lon used
    nearest_lat = data_tp_pt.latitude.values
    nearest_lng = data_tp_pt.longitude.values

    return monthly_climatology, nearest_lat, nearest_lng


# Call our function
clim, nearest_lat, nearest_lng = precipMonthlyClimatology()

And finally, let’s set up the plot nicely. This can easily be customised, but for now we do something similar to ERA Explorer

# Determine suffix for latitude (N/S) and longitude (E/W) based on their sign
latSuffix = 'N' if nearest_lat > 0 else 'S'  # 'N' for +ve latitude, 'S' for -ve latitude
lngSuffix = 'E' if nearest_lng > 0 else 'W'  # 'E' for +ve longitude, 'W' for -ve longitude

# Create a new figure with specified size
plt.figure(figsize=(8, 6))

# Plot the data as a bar chart
bars = plt.bar(
    clim.index,               # X-axis values (e.g., months)
    clim['tp'],               # Y-axis values (e.g., precipitation)
    width=0.8,                # Width of each bar
    color='#596CAA',          # Bar fill color
    edgecolor='white',        # Color of the bar edges
    label=f'{abs(nearest_lat):.2f} °{latSuffix:s}, {abs(nearest_lng):.2f} °{lngSuffix:s}'
)

# Add labels on top of each bar
for bar in bars:
    height = bar.get_height()
    plt.text(
        bar.get_x() + bar.get_width() / 2,  # Center of the bar
        height + 0.5,                       # Slightly above the bar
        f'{height:.0f}',                    # Rounded bar height
        ha='center',                        # Horizontal alignment
        va='bottom',                        # Vertical alignment
        fontsize=10,                        # Font size
        color='black'                       # Text color
    )

# Add a legend with a transparent background
plt.legend(framealpha=0)

# Customize axis labels
plt.xlabel('Month', fontsize=12)
plt.ylabel('Precipitation [mm]', fontsize=12)

# Set custom x-axis tick labels to month abbreviations
plt.xticks(
    ticks=clim.index,
    labels=[calendar.month_abbr[i] for i in clim.index]
)

# Add a title to the plot
plt.title(
    f'Monthly Precipitation Climatologies from ERA5 '
    f'({date_range[0][:4]}-{date_range[1][:4]})',
    fontsize=14
)

# Adjust layout to ensure no elements overlap
plt.tight_layout()

# Display the final plot
plt.show()