Calculate percentage of days when Tmean > 90th percentile (TG90p)

Example notebook that runs icclim.

The example calculates the percentage of days when tas exceeds the 90th percentile (TG90p) for the dataset chosen by the user on C4I.

We assume to have the tas variable in netCDF files in a ./data/latest folder for model CMCC and for one member r1i1p1f1.

The data can be downloaded using the metalink provided with this notebook. The data described in a .metalink file can be downloaded with tools such as aria2 or a browser plugin such as DownThemAll! If you wish to use a different dataset, you can use the climate 4 impact portal to search and select the data you wish to use and a metalink file to the ESGF data will be provided.

The data is read using xarray and a plot of the time series over a specific region is generated, as well as an average spatial map. Several output types examples are shown.

The dataset expected for this notebook is the tas variable needed to calculate TG90p for one climate model, one experiment and one member. The time period should be continuous.

The following time period is considered: 2081-01-01 to 2100-12-31 using the period 1981-01-01 to 2000-12-31 as a reference. Plots are shown over European region.

Install packages

import datetime
from pathlib import Path

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import numpy as np
import xarray as xr

import icclim

%pip install icclim matplotlib nc_time_axis

Specification of index parameters

# studied period
dt1 = datetime.datetime(2081, 1, 1, tzinfo=datetime.timezone.utc)
dt2 = datetime.datetime(2100, 12, 31, tzinfo=datetime.timezone.utc)

# reference period
dt1r = datetime.datetime(1981, 1, 1, tzinfo=datetime.timezone.utc)
dt2r = datetime.datetime(2000, 12, 31, tzinfo=datetime.timezone.utc)

out_f = "tg90p_icclim.nc"
data_dir = Path("data/latest")
filenames = [str(f) for f in data_dir.glob("tas_day_CMCC-ESM2*.nc")]
filenames

icclim.index(
    index_name="TG90p",
    in_files=filenames,
    slice_mode="JJA",
    base_period_time_range=[dt1r, dt2r],
    time_range=[dt1, dt2],
    out_unit="%",
    out_file=out_f,
    logs_verbosity="HIGH",
)

Plot preparation

with xr.open_dataset(out_f, decode_times=False) as ds:
    tg90_xr = ds
    ds["time"] = xr.decode_cf(ds).time

# Select a single x,y combination from the data
longitude = tg90_xr["TG90p"]["lon"].sel(lon=3.5, method="nearest")
latitude = tg90_xr["TG90p"]["lat"].sel(lat=44.2, method="nearest")

Subset and Plot

# Slice the data spatially using a single lat/lon point
one_point = tg90_xr["TG90p"].sel(lat=latitude, lon=longitude)
# Use xarray to create a quick time series plot
one_point.plot.line()
plt.show()

# You can clean up your plot as you wish using standard matplotlib approaches
f, ax = plt.subplots(figsize=(12, 6))
one_point.plot.line(
    hue="lat",
    marker="o",
    ax=ax,
    color="grey",
    markerfacecolor="purple",
    markeredgecolor="purple",
)
ax.set(title="Time Series For a Single Lat / Lon Location")

plt.show()

# Convert to dataframe -- then this can easily be exported to a csv
one_point_df = one_point.to_dataframe()
# View just the first 5 rows of the data
one_point_df.head()

# Export data to .csv file

# Time subsetting: this is just an example on how to do it
start_date = "2081-01-01"
end_date = "2082-12-31"

tg90_two = tg90_xr["TG90p"].sel(time=slice(start_date, end_date))

# Quickly plot the data using xarray.plot()
tg90_two.plot(x="lon", y="lat", col="time", col_wrap=1)

plt.suptitle("Two Time Steps of TG90P", y=1.03)
plt.show()

# Set spatial extent and centre
central_lat = 47.0
central_lon = 1.0
extent = [-30, 30, 30, 56]  # Western Europe

tg90 = tg90_xr["TG90p"]

# Calculate time average
tg90_avg = tg90.mean(dim="time", keep_attrs=True)

# Set plot projection
map_proj = ccrs.AlbersEqualArea(
    central_longitude=central_lon, central_latitude=central_lat
)

# Define plot
f, ax = plt.subplots(figsize=(14, 6), subplot_kw={"projection": map_proj})

# Plot data with proper colormap scale range
levels = np.arange(0, 1.0, 0.1)
p = tg90_avg.plot(levels=levels, cmap="RdBu_r", transform=ccrs.PlateCarree())

# Plot information
plt.suptitle(
    "Percentage of days when Tas > 90th percentil Period 2081-2100 Reference 1981-2000 TG90P",
    y=1,
)

# Add the coastlines to axis and set extent
ax.coastlines()
ax.gridlines()
ax.set_extent(extent)

# Save plot as png
plt.savefig("c4i_tg90p_icclim.png")

# Re-order longitude so that there is no blank line at 0 deg because 0 deg is within our spatial selection
tg90_avg.coords["lon"] = (tg90_avg.coords["lon"] + 180) % 360 - 180
tg90_avg = tg90_avg.sortby(tg90_avg.lon)

# Define plot
f, ax = plt.subplots(figsize=(14, 6), subplot_kw={"projection": map_proj})

# Define colorscale
levels = np.arange(0, 1.0, 0.1)

# Contours lines
p = tg90_avg.plot.contour(
    levels=levels, colors="k", linewidths=0.5, transform=ccrs.PlateCarree()
)

# Contour filled colors
p = tg90_avg.plot.contourf(
    levels=levels, cmap="RdBu_r", extend="both", transform=ccrs.PlateCarree()
)

# Plot information
plt.suptitle(
    "Percentage of days when tas > 90th percentile - Period 2081-2100 - Reference 1981-2000 - TG90P",
    y=1,
)

# Add the coastlines to axis and set extent
ax.coastlines()
ax.gridlines()
ax.set_extent(extent)

# Save plot as png
plt.savefig("c4i_tg90p_contours_icclim.png")