Calculate the monthly climatological temperature thresholds 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 = "2m_temperature"
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 temperature climatology
def temperatureMonthlyClimatology():
"""
Processes temperature data to calculate monthly climatology and thresholds.
This function reads temperature data from a NetCDF file, processes it to calculate
daily minimum and maximum temperatures, and then computes monthly climatology
statistics. It also calculates the percentage of frost days, tropical nights, and
summer days per month over a 30-year period.
Returns:
tuple: A tuple containing the following elements:
- monthly_climatology (pd.DataFrame): Climatological averages for daily
minimum and maximum temperatures per month, converted to Celsius.
- frost_day_climatology (pd.Series): Average percentage of frost days
per month.
- tropical_night_climatology (pd.Series): Average percentage of tropical
nights per month.
- summer_day_climatology (pd.Series): Average percentage of summer days
per month.
"""
data_t2m_pt = data.t2m
# Convert the time coordinate to a pandas datetime index
time_index = pd.to_datetime(data_t2m_pt.valid_time.values)
# Create a DataFrame for easier manipulation
df = pd.DataFrame(data_t2m_pt.values, index=time_index, columns=['t2m'])
# Resample to find daily minimum and maximum
daily_min = df.resample('D').min()
daily_max = df.resample('D').max()
# Combine the daily min and max into a single DataFrame
daily_stats = pd.DataFrame({
'daily_min': daily_min['t2m'],
'daily_max': daily_max['t2m']
})
# Define the frost night threshold (e.g., 0°C in Kelvin)
frostDayThreshold = 0 + 273.15
tropicalNightThreshold = 20 + 273.15
summerDayThreshold = 25 + 273.15
# Check if the daily minimum is below the threshold
daily_min['frost_day'] = daily_min['t2m'] < frostDayThreshold
daily_min['tropical_night'] = daily_min['t2m'] > tropicalNightThreshold
daily_max['summer_day'] = daily_max['t2m'] > summerDayThreshold
# Count the number of frost days (days below the threshold) per month
frost_day_counts = daily_min.resample('ME').sum()['frost_day']
tropical_night_counts = daily_min.resample('ME').sum()['tropical_night']
summer_day_counts = daily_max.resample('ME').sum()['summer_day']
# Count the total number of days per month
total_days_per_month = daily_min.resample('ME').size()
# Calculate the percentage of frost nights per month
frost_day_pc = (frost_day_counts / total_days_per_month) * 100
tropical_night_pc = (tropical_night_counts / total_days_per_month) * 100
summer_day_pc = (summer_day_counts / total_days_per_month) * 100
# Convert the index to month
frost_day_pc.index = frost_day_pc.index.month
tropical_night_pc.index = tropical_night_pc.index.month
summer_day_pc.index = summer_day_pc.index.month
# Group by month and calculate the average percentage over the 30-year period
frost_day_climatology = frost_day_pc.groupby(frost_day_pc.index).mean()
tropical_night_climatology = tropical_night_pc.groupby(tropical_night_pc.index).mean()
summer_day_climatology = summer_day_pc.groupby(summer_day_pc.index).mean()
# Extract the month from DateTimeIndex
daily_stats['month'] = daily_stats.index.month
# Group by month and calculate climatological averages for daily min and max
grouped_by_month = daily_stats.groupby('month')
monthly_climatology = grouped_by_month.mean() - 273.15 # Convert to Celsius
# Get the actual lat/lon used
nearest_lat = data_t2m_pt.latitude.values
nearest_lng = data_t2m_pt.longitude.values
return monthly_climatology, frost_day_climatology, \
tropical_night_climatology, summer_day_climatology, \
nearest_lat, nearest_lng
# Call our function
clim, frost_clim, tropical_clim, summer_clim, \
nearest_lat, nearest_lng = temperatureMonthlyClimatology()
And finally, let’s set up the plot nicely. This can easily be customised, but for now we do something similar to ERA Explorer
# Set line width, marker properties, and edge properties for the plot
lw = 4 # Line width for the plot lines
marker = 'o' # Marker style for the data points
markersize = 12 # Size of the markers
markeredgecolor = 'white' # Color of the marker edges
markeredgewidth = 2 # Width of the marker edges
# Determine the suffix for latitude (N/S) and longitude (E/W) based on their signs
latSuffix = 'N' if nearest_lat > 0 else 'S'
lngSuffix = 'E' if nearest_lng > 0 else 'W'
# Create a new figure with specified size
plt.figure(figsize=(7, 6))
# Plot frost days climatology with customized markers and lines
plt.plot(
clim.index,
frost_clim,
marker=marker,
markersize=markersize,
markeredgecolor=markeredgecolor,
markeredgewidth=markeredgewidth,
label=(
f'Frost days ({abs(nearest_lat):.2f} °{latSuffix}, '
f'{abs(nearest_lng):.2f} °{lngSuffix})'
),
color='#25408F',
lw=lw,
clip_on=False
)
# Plot summer days climatology with the same styling as frost days
plt.plot(
clim.index,
summer_clim,
marker=marker,
markersize=markersize,
markeredgecolor=markeredgecolor,
markeredgewidth=markeredgewidth,
label=(
f'Summer days ({abs(nearest_lat):.2f} °{latSuffix}, '
f'{abs(nearest_lng):.2f} °{lngSuffix})'
),
color='#FFC900',
lw=lw,
clip_on=False
)
# Plot tropical nights climatology with the same styling as other plots
plt.plot(
clim.index,
tropical_clim,
marker=marker,
markersize=markersize,
markeredgecolor=markeredgecolor,
markeredgewidth=markeredgewidth,
label=(
f'Tropical nights ({abs(nearest_lat):.2f} °{latSuffix}, '
f'{abs(nearest_lng):.2f} °{lngSuffix})'
),
color='#CF7E92',
lw=lw,
clip_on=False
)
# Add legend to the plot with no background frame
plt.legend(framealpha=0)
# Customize x and y axis labels
plt.xlabel('Month', fontsize=12)
plt.ylabel('Percentage of days per month [%]', fontsize=12)
# Set x-axis ticks to correspond to month indices and display month abbreviations
plt.xticks(
ticks=clim.index,
labels=[calendar.month_abbr[i] for i in clim.index]
)
# Set the limits for the y-axis to range from 0 to 100
plt.ylim(0, 100)
# Add a title with the specified date range and custom font size
plt.title(
f'Monthly Temperature Threshold Climatologies from ERA5 '
f'({date_range[0][:4]}-{date_range[1][:4]})',
fontsize=14
)
# Adjust layout to ensure all elements are properly displayed and do not overlap
plt.tight_layout()
# Display the final plot
plt.show()
