Python matplotlib example Hovmoeller diagram#

Software requirements:

  • Python 3

  • numpy

  • xarray

  • matplotlib

  • cartopy

Run the Hovmoeller diagram example script:

python Hovmoeller_diagram_with_additional_map.py

Script Hovmoeller_diagram_with_additional_map.py:

#!/usr/bin/env python
# coding: utf-8

# DKRZ tutorial
#
# Hovmöller diagram with additional map
#
# A Hovmöller diagram is often used for meterological data. The diagram
# shows the contours of a spatial variable in time and space:
#
#     x-axis: longitude or latitude
#     y-axis: time
#
# Here for demonstration purposes, we use the variable northward wind, get
# the data slice along a specified range of latitudes (40° - 55°) and plot
# it against time.
#
# Learning content
#
# - read netCDF file
# - extract slices
# - convert pressure levels to height
# - create a Hovmöller diagram
# - add a map on top of the plot
#
# 06.12.2022  copyright DKRZ, kmf

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

def main():
    # Read the CMIP5 atmos historical dataset of the northward wind.

    infile2 = os.environ['HOME']+\
              '/data/CMIP5/atmos/va_Amon_MPI-ESM-LR_historical_r1i1p1_198001-198912.nc'

    dsh = xr.open_dataset(infile2)

    # The grid of the data contains longitudes from 0° to 360°, but we want to
    # have them from -180° to 180°.

    dsh = dsh.assign_coords(lon=(((dsh.lon + 180) % 360) - 180)).sortby('lon')

    # Select variable **va** from dataset and convert the pressure levels from
    # Pa to hPa, and overwrite plev in the dataset.

    vnwind = dsh.va

    levels = dsh.plev / 100
    levels.attrs['units'] = 'hPa'

    vnwind['plev'] = levels

    # Select the latitude range to be used and select one pressure level.

    lat_slice = [40., 55.]
    sel_lev = 300

    # Next, get the slices of the data.

    data = vnwind.sel(plev=sel_lev,
                      lat=lat_slice,
                      method='nearest')

    # Compute the weights and the weighted mean over latitude dimension.

    weights = np.cos(np.deg2rad(data.lat))
    weights.name = "weights"

    data_weighted = data.weighted(weights)
    data_weighted_mean = data_weighted.mean('lat')
    data_weighted_mean.name = 'va'

    # Create datetime objects

    vtimes = vnwind.time.values.astype('datetime64[ms]').astype('O')

    # Create the Hovmöller plot and add a small map of the region above it.

    projection = ccrs.PlateCarree()

    fig = plt.figure(figsize=(12, 10), constrained_layout=False)

    # define space for the two plots
    gs = fig.add_gridspec(nrows=4, ncols=1, hspace=0)

    #-- create the small upper map
    ax1 = plt.subplot(gs[0, 0], projection=projection)
    ax1.set_extent([-180., 180., lat_slice[0], lat_slice[1]],
                   crs=ccrs.PlateCarree())
    ax1.add_feature(cfeature.LAND)
    ax1.add_feature(cfeature.OCEAN)
    ax1.add_feature(cfeature.BORDERS)
    ax1.add_feature(cfeature.COASTLINE)
    ax1.gridlines(draw_labels=True, lw=0.01,crs=ccrs.PlateCarree())

    plt.title(u'Hovmöller Diagram', loc='left')
    plt.title('MPI-ESM-LR RCP-4.5', loc='right')

    #-- create lower Hovmoeller diagram
    clevs = np.arange(-30, 31, 5)

    ax2 = fig.add_subplot(gs[1:, 0], sharex=ax1)

    cf = ax2.contourf(vnwind.lon.values, vtimes, data_weighted_mean,
                      clevs,
                      cmap='viridis',
                      extend='both')

    cs = ax2.contour(vnwind.lon.values, vtimes, data_weighted_mean,
                     clevs,
                     colors='k',
                     linewidths=0.25)

    cbar = plt.colorbar(cf,
                        orientation='horizontal',
                        pad=0.05,
                        shrink=0.7,
                        aspect=50,
                        extendrect=True)
    cbar.set_label('m $s^{-1}$')

    ax2.set_xticks([-180, -120, -60, 0, 60, 120, 180])
    ax2.set_xticklabels([u'180\N{DEGREE SIGN}', u'120\N{DEGREE SIGN}W',
                         u'60\N{DEGREE SIGN}W', u'0\N{DEGREE SIGN}',
                         u'60\N{DEGREE SIGN}E', u'120\N{DEGREE SIGN}E',
                         u'180\N{DEGREE SIGN}'])
    ax2.set_yticks(vtimes[4::8])
    ax2.set_yticklabels(vtimes[4::8])

    plt.title(str(sel_lev)+'-hPa V-wind   latitudes: '+str(lat_slice),
              loc='left',
              fontsize=10)
    plt.title('Time Range: {0:%Y%m%d %HZ} - {1:%Y%m%d %HZ}'.format(vtimes[0], vtimes[-1]),
              loc='right',
              fontsize=10)

    plt.savefig('plot_hovmoeller_2.png',
                 bbox_inches='tight',
                 facecolor='white',
                 dpi=100)

if __name__ == '__main__':
    main()

Plot result:

image0