This post discuss Hovmöller diagrams drawing in python. Hovmöller diagrams is used to plot the time evolution of vertical profiles of scalar quantities such as temperature, density, or concentrations of constituents in the atmosphere or ocean. The axes of a Hovmöller diagram are typically longitude or latitude and time with the value of some field represented through color or shading.
For this exercice Monthly variation of temperature over West-Africa following a latiude and longitude gradient will be ploted.
Load packages¶
In [1]:
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.gridspec as gridspec
import matplotlib.pyplot as plt
import metpy.calc as mpcalc
import numpy as np
import xarray as xr
import pandas as pd
In [3]:
#Setting Working directory
import os
os.chdir("D:/Weather_Processed_data" )
os.getcwd()
Out[3]:
Import data¶
Near-surface temperature from CRU TS4.04 daily datasets is used. It can be downloaded from this :http://dap.ceda.ac.uk/thredds/fileServer/badc/cru/data/cru_ts/cru_ts_4.04/data/tmp/cru_ts4.04.1901.2019.tmp.dat.nc.gz
In [4]:
# import the dataset
df = xr.open_dataset('cru_ts4.04.1901.2019.tmp.dat.nc')
In [5]:
df #print content
Out[5]:
Study area - West-Africa¶
In [6]:
# Latitude and Longitude extent of the study area
lat_slice = slice(5,30)
lon_slice = slice(-15,30)
# Get data, selecting time, lat/lon slice
data=df['tmp'].sel(lat=lat_slice,lon=lon_slice)
In [7]:
# Start figure
fig = plt.figure(figsize=(10,50))
# Use gridspec to help size elements of plot; small top plot and big bottom plot
gs = gridspec.GridSpec(nrows=2, ncols=1, hspace=0.03)
# Tick labels
# Define x, y axis labels
x_tick_labels = [u'20\N{DEGREE SIGN}W',u'15 \N{DEGREE SIGN}W', u'10\N{DEGREE SIGN}W', u'5\N{DEGREE SIGN}W',
u'0\N{DEGREE SIGN}',u'5\N{DEGREE SIGN}E',u'10\N{DEGREE SIGN}E',u'15\N{DEGREE SIGN}E',
u'20\N{DEGREE SIGN}E',
u'25\N{DEGREE SIGN}E',u'30\N{DEGREE SIGN}E']
y_tick_labels = [u'5\N{DEGREE SIGN}N', u'10 \N{DEGREE SIGN}N',
u'15\N{DEGREE SIGN}N',u'20\N{DEGREE SIGN}N',u'25\N{DEGREE SIGN}N',u'30\N{DEGREE SIGN}N']
# Set geographic reference
ax1 = fig.add_subplot(gs[0, 0], projection=ccrs.PlateCarree(central_longitude=0))
ax1.set_extent([-15,30,0,30], ccrs.PlateCarree(central_longitude=0)) #lon min, max, lat min, max
ax1.set_xticks([-20,-15,-10, -5, 0, 5,10,20,25,30])
ax1.set_xticklabels(x_tick_labels)
ax1.set_yticks([5,10,15,20,25,30])
ax1.set_yticklabels(y_tick_labels)
ax1.grid(linestyle='dotted', linewidth=2)
# Add geopolitical boundaries for map reference
ax1.add_feature(cfeature.COASTLINE.with_scale('50m'))
ax1.add_feature(cfeature.LAKES.with_scale('50m'), color='black', linewidths=0.5)
ax1.add_feature(cfeature.BORDERS, linewidths=1)
# Set some titles
plt.title('Hovmoller Diagram over West-Africa', loc='left',fontsize=15)
plt.show()
Compute monthly mean around the datasets¶
In [8]:
monmean=data.groupby('time.month').mean()
monmean
Out[8]:
Compute weights over latitude dimension¶
In [9]:
# Get times and make array of datetime objects
month = monmean.month.values
# Compute weights and take weighted average over latitude dimension
weights = np.cos(np.deg2rad(monmean.lat.values))
avg_data = (monmean * weights[None, :, None]).sum(dim='lat') / np.sum(weights)
#Specify longitude values for chosen domain
lons = monmean.lon.values
Monthly Longitudunal over time series¶
In [10]:
# Define x, y axis labels
x_tick_labels = [u'15 \N{DEGREE SIGN}W', u'10\N{DEGREE SIGN}W', u'5\N{DEGREE SIGN}W',
u'0\N{DEGREE SIGN}',
u'5\N{DEGREE SIGN}E',u'10\N{DEGREE SIGN}E',u'15\N{DEGREE SIGN}E',u'20\N{DEGREE SIGN}E',
u'25\N{DEGREE SIGN}E',u'30\N{DEGREE SIGN}E']
y_tick_labels=['JAN','FEB','MAR','APR','MAY','JUN','JUL','AUG','SEP','OCT','NOV','DEV']
fig = plt.figure(figsize=(10,7))
#gs = gridspec.GridSpec(nrows=1, ncols=1, hspace=0.03)
#Bottom plot for Hovmoller diagram
ax= plt.axes()
ax.invert_yaxis() # Reverse the time order to do oldest first
# Plot of chosen variable averaged over latitude and slightly smoothed
z=mpcalc.smooth_n_point(avg_data, 9, 2)
clevs= np.arange(z.min(),z.max(),1) # color levels
#cf = ax.contourf(lons,month, mpcalc.smooth_n_point(avg_data, 9, 2), clevs, cmap=plt.cm.bwr, extend='both')
cf = ax.contourf(lons,month, mpcalc.smooth_n_point(avg_data, 9, 2), clevs, cmap=plt.cm.RdBu_r, extend='both')
cs = ax.contour(lons,month, mpcalc.smooth_n_point(avg_data, 9, 2), clevs, colors='k', linewidths=False)
cbar = plt.colorbar(cf, orientation='horizontal', pad=0.1, aspect=50, extend=True)
cbar.set_label('Deg [C]')
# Make some ticks and tick labels
ax.set_xticks([-15,-10, -5, 0, 5, 10,15,20,25,30])
ax.set_xticklabels(x_tick_labels)
ax.set_yticks(month)
ax.set_yticklabels(y_tick_labels)
#import time series
time=df.time.values.astype('datetime64[ms]').astype('O')
# Set some titles
plt.title('Mean Temperature', loc='left', fontsize=10)
plt.title('Time Range: {0:%Y/%m/%d} - {1:%Y/%m/%d}'.format(time[0], time[-1]),
loc='right', fontsize=10)
plt.show()
In [11]:
#Get more Colormaps by doing a tab completion on the plt.cm module: plt.cm.<TAB>
Monthly Latidunal over time series¶
In [12]:
# Compute weights and take weighted average over latitude dimension
weights = np.cos(np.deg2rad(monmean.lon.values))
avg_data = (monmean * weights[None,None,:]).sum(dim='lon') / np.sum(weights)
# Specify latitude values for chosen domain
lats = data.lat.values
In [13]:
#Reordering dimensions
avg_data=avg_data.transpose('lat', 'month')
# set figure dimensions and axis
fig = plt.figure(figsize=(15,8))
ax2 = plt.axes()
# Plot of chosen variable averaged over longitude and slightly smoothed
z=mpcalc.smooth_n_point(avg_data, 9, 2)
clevs= np.arange(z.min(),z.max(),1) # color levels
cf = ax2.contourf(month,lats,z, clevs, cmap=plt.cm.Spectral_r, extend='both')
cs = ax2.contour(month,lats,z, clevs, colors='k', linewidths=0)
cbar = plt.colorbar(cf, orientation='vertical', pad=0.04, aspect=50, extendrect=True)
cbar.set_label('Deg [C]')
#Define x and y labels plots
x_tick_labels=['JAN','FEB','MAR','APR','MAY','JUN','JUL','AUG','SEP','OCT','NOV','DEV']
ax2.set_xticks(month)
ax2.set_xticklabels(x_tick_labels)
# Set some titles
plt.title('Mean Temperature', loc='left', fontsize=10)
plt.title('Time Range: {0:%Y/%m/%d} - {1:%Y/%m/%d}'.format(time[0], time[-1]),
loc='right', fontsize=10)
plt.ylabel("Latitude \N{DEGREE SIGN} N", Fontsize=10)
plt.xlabel("Monthly Mean", Fontsize=10)
plt.show()
In [14]:
#part of this code is inspired from this post: https://unidata.github.io/python-gallery/examples/Hovmoller_Diagram.html
Thank you very much. This is very nice !
ReplyDeleteNicely done bro.
ReplyDeleteWell done!!
ReplyDeleteGood right
ReplyDeleteVery nice and useful Thierry, thank you. Keep it up
ReplyDelete