Plot QVA Probabilistic¶
Plots a probabilistic Quantitative Volcanic Ash Concentration forecast with the agreed colormaps per threshold:
| Threshold band (mg/m3) | Colormap |
|---|---|
| [0.2 - 2[ | Blues |
| [2 - 5[ | YlOrBr |
| [5 - 10[ | Reds |
| [10 - Inf[ | RbPu |
Intervals are lower-bound inclusive, upper-bound exclusive.
The netcdf file used in this example (qva_prob_Ruapehu_202501081200.nc), as well as this notebook are is available in https://github.com/rosatrancoso/notebook-share/blob/master/docs/qva/.
In [1]:
Copied!
# Contents of the netcdf file
! ncdump -h qva_prob_Ruapehu_202501081200.nc
# Contents of the netcdf file
! ncdump -h qva_prob_Ruapehu_202501081200.nc
netcdf qva_prob_Ruapehu_202501081200 {
dimensions:
time = UNLIMITED ; // (24 currently)
latitude = 361 ;
longitude = 441 ;
levels = 12 ;
thresholds = 4 ;
variables:
float latitude(latitude) ;
latitude:_FillValue = NaNf ;
latitude:long_name = "latitude degrees north from the equator" ;
latitude:units = "degrees_north" ;
latitude:point_spacing = "even" ;
float longitude(longitude) ;
longitude:_FillValue = NaNf ;
longitude:long_name = "longitude degrees east from the greenwich meridian" ;
longitude:units = "degrees_east" ;
longitude:point_spacing = "even" ;
double levels(levels) ;
levels:_FillValue = NaN ;
levels:long_name = "Top of flight level layer" ;
levels:units = "feet" ;
levels:axis = "Z" ;
levels:positive = "up" ;
int64 time(time) ;
time:axis = "T" ;
time:long_name = "time" ;
time:units = "hours since 2025-01-08T13:00:00" ;
time:calendar = "proleptic_gregorian" ;
double thresholds(thresholds) ;
thresholds:_FillValue = NaN ;
float frequency_of_excedance(thresholds, time, levels, latitude, longitude) ;
frequency_of_excedance:_FillValue = NaNf ;
frequency_of_excedance:long_name = "Frequency of excedance based on 30 members of GEFS" ;
frequency_of_excedance:units = "fraction" ;
frequency_of_excedance:thresholds = 0.2, 2., 5., 10. ;
// global attributes:
:title = "HYSPLIT Model Concentration Output" ;
:Conventions = "CF-1.5" ;
:volcano_name = "Ruapehu" ;
:eruption_lon_degrees = 175.57f ;
:eruption_lat_degrees = -39.28f ;
:eruption_vent_meters_msl = 2797.f ;
:eruption_height_meters_msl = 12500.f ;
:eruption_mass_eruption_rate_kgs = 87500. ;
:eruption_start_time = "2025-01-08T12:00:00" ;
:eruption_duration_s = 3600. ;
}
In [2]:
Copied!
import numpy as np
import pandas as pd
import xarray as xr
import matplotlib.pyplot as plt
import matplotlib as mpl
import cartopy.crs as ccrs
import cartopy.feature as cf
from math import floor,ceil
# x-positions of gridlines (need to specify when over 180 deg)
xlocs = [100,110, 120, 130, 140, 150, 160, 170, 179.99999, -170, -160, -150, -140, -130, -120, -110, -100]
# central longitude for projection
clon = 180
# mask values below "numerical zero"
almost_zero = 1e-16
import numpy as np
import pandas as pd
import xarray as xr
import matplotlib.pyplot as plt
import matplotlib as mpl
import cartopy.crs as ccrs
import cartopy.feature as cf
from math import floor,ceil
# x-positions of gridlines (need to specify when over 180 deg)
xlocs = [100,110, 120, 130, 140, 150, 160, 170, 179.99999, -170, -160, -150, -140, -130, -120, -110, -100]
# central longitude for projection
clon = 180
# mask values below "numerical zero"
almost_zero = 1e-16
Colormap settings¶
In [3]:
Copied!
thresholds = np.array([0.2,2,5,10])
print(f'concentration thresholds (mg/m3)= {thresholds}')
# Colormaps per threshold band
colormaps = {
0.2: 'Blues',
2: 'YlOrBr',
5: 'Reds',
10: 'RdPu',
}
# contour levels for colorbar
contour_levels = [0.01,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1]
norm = mpl.colors.BoundaryNorm(boundaries=contour_levels, ncolors=256)
for i in range(len(thresholds)):
print(i, thresholds[i], colormaps[thresholds[i]])
fig, ax = plt.subplots(figsize=(6, 1), layout='constrained')
fig.colorbar(mpl.cm.ScalarMappable(cmap=colormaps[thresholds[i]], norm=norm),
cax=ax, orientation='horizontal', label=colormaps[thresholds[i]])
thresholds = np.array([0.2,2,5,10])
print(f'concentration thresholds (mg/m3)= {thresholds}')
# Colormaps per threshold band
colormaps = {
0.2: 'Blues',
2: 'YlOrBr',
5: 'Reds',
10: 'RdPu',
}
# contour levels for colorbar
contour_levels = [0.01,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1]
norm = mpl.colors.BoundaryNorm(boundaries=contour_levels, ncolors=256)
for i in range(len(thresholds)):
print(i, thresholds[i], colormaps[thresholds[i]])
fig, ax = plt.subplots(figsize=(6, 1), layout='constrained')
fig.colorbar(mpl.cm.ScalarMappable(cmap=colormaps[thresholds[i]], norm=norm),
cax=ax, orientation='horizontal', label=colormaps[thresholds[i]])
concentration thresholds (mg/m3)= [ 0.2 2. 5. 10. ] 0 0.2 Blues 1 2.0 YlOrBr 2 5.0 Reds 3 10.0 RdPu
Read input file¶
In [4]:
Copied!
ds = xr.open_dataset('qva_prob_Ruapehu_202501081200.nc')
ds
ds = xr.open_dataset('qva_prob_Ruapehu_202501081200.nc')
ds
Out[4]:
<xarray.Dataset> Size: 734MB
Dimensions: (latitude: 361, longitude: 441, levels: 12,
time: 24, thresholds: 4)
Coordinates:
* latitude (latitude) float32 1kB -84.28 -84.03 ... 5.47 5.72
* longitude (longitude) float32 2kB 120.6 120.8 ... 230.3 230.6
* levels (levels) float64 96B 5e+03 1e+04 ... 5.5e+04 6e+04
* time (time) datetime64[ns] 192B 2025-01-08T13:00:00 .....
* thresholds (thresholds) float64 32B 0.2 2.0 5.0 10.0
Data variables:
frequency_of_excedance (thresholds, time, levels, latitude, longitude) float32 734MB ...
Attributes:
title: HYSPLIT Model Concentration Output
Conventions: CF-1.5
volcano_name: Ruapehu
eruption_lon_degrees: 175.57
eruption_lat_degrees: -39.28
eruption_vent_meters_msl: 2797.0
eruption_height_meters_msl: 12500.0
eruption_mass_eruption_rate_kgs: 87500.0
eruption_start_time: 2025-01-08T12:00:00
eruption_duration_s: 3600.0Get full extent of the plume at all levels and times¶
This helps to have the same zoom extent in all plots that fits the data.
In [5]:
Copied!
# determine full extent of the plume at all levels and times
da = ds['frequency_of_excedance']
da_m = da.where(da > almost_zero, drop=True)
print(f'before mask = {da.shape}')
print(f'after mask = {da_m.shape}')
bnd = [floor(da_m.longitude.values.min()), ceil(da_m.longitude.values.max()),
floor(da_m.latitude.values.min()), ceil(da_m.latitude.values.max())]
print(f'full extent of plume = {bnd}')
# determine full extent of the plume at all levels and times
da = ds['frequency_of_excedance']
da_m = da.where(da > almost_zero, drop=True)
print(f'before mask = {da.shape}')
print(f'after mask = {da_m.shape}')
bnd = [floor(da_m.longitude.values.min()), ceil(da_m.longitude.values.max()),
floor(da_m.latitude.values.min()), ceil(da_m.latitude.values.max())]
print(f'full extent of plume = {bnd}')
before mask = (4, 24, 12, 361, 441) after mask = (4, 24, 12, 42, 124) full extent of plume = [175, 207, -42, -31]
Plot¶
Choose time and level to plot
In [6]:
Copied!
# time index
it = 2
# level index of top of the layer
ilev = 7
print(ds.time[it].values)
print(f'FL{ds.levels[ilev].values/100:.0f}')
# time index
it = 2
# level index of top of the layer
ilev = 7
print(ds.time[it].values)
print(f'FL{ds.levels[ilev].values/100:.0f}')
2025-01-08T15:00:00.000000000 FL400
Prepare plot title
In [10]:
Copied!
volcano_name = ds.attrs['volcano_name']
olon = ds.attrs['eruption_lon_degrees']
olat = ds.attrs['eruption_lat_degrees']
h_vent = ds.attrs['eruption_vent_meters_msl']
h_top = ds.attrs['eruption_height_meters_msl']
mer = ds.attrs['eruption_mass_eruption_rate_kgs']
dur = ds.attrs['eruption_duration_s']
time = pd.to_datetime(ds.time[it].values)
eruption_start_time = pd.to_datetime(ds.attrs['eruption_start_time'])
lead_time = (time - eruption_start_time).total_seconds()/3600.
level = ds.levels[ilev].values
level_below = ds.levels[ilev-1].values if ilev > 0 else 0
level = (level/100).astype(int)
level_below = (level_below/100).astype(int)
tit_str1 = f'{volcano_name} ({olon:.4f} E, {olat:.4f} N, {h_vent:.1f} m asl)'
tit_str2 = f'H = {h_top/1000.:.1f} km asl; MER = {mer:.1e} kg/s, D = {dur/3600.:.0f} h'
tit_str3 = f'Valid at {time} (+{lead_time:.0f} h after eruption) for FL{level_below}-FL{level}'
tit_str = f'{tit_str1}\n{tit_str2}\n{tit_str3}'
print(tit_str)
volcano_name = ds.attrs['volcano_name']
olon = ds.attrs['eruption_lon_degrees']
olat = ds.attrs['eruption_lat_degrees']
h_vent = ds.attrs['eruption_vent_meters_msl']
h_top = ds.attrs['eruption_height_meters_msl']
mer = ds.attrs['eruption_mass_eruption_rate_kgs']
dur = ds.attrs['eruption_duration_s']
time = pd.to_datetime(ds.time[it].values)
eruption_start_time = pd.to_datetime(ds.attrs['eruption_start_time'])
lead_time = (time - eruption_start_time).total_seconds()/3600.
level = ds.levels[ilev].values
level_below = ds.levels[ilev-1].values if ilev > 0 else 0
level = (level/100).astype(int)
level_below = (level_below/100).astype(int)
tit_str1 = f'{volcano_name} ({olon:.4f} E, {olat:.4f} N, {h_vent:.1f} m asl)'
tit_str2 = f'H = {h_top/1000.:.1f} km asl; MER = {mer:.1e} kg/s, D = {dur/3600.:.0f} h'
tit_str3 = f'Valid at {time} (+{lead_time:.0f} h after eruption) for FL{level_below}-FL{level}'
tit_str = f'{tit_str1}\n{tit_str2}\n{tit_str3}'
print(tit_str)
Ruapehu (175.5700 E, -39.2800 N, 2797.0 m asl) H = 12.5 km asl; MER = 8.8e+04 kg/s, D = 1 h Valid at 2025-01-08 15:00:00 (+3 h after eruption) for FL350-FL400
Plots for each threshold
In [19]:
Copied!
fig, axs = plt.subplots(len(thresholds), 1, figsize=(10,12), subplot_kw={"projection": ccrs.PlateCarree(clon)})
for ithreshold in range(len(thresholds)):
ax = axs.flatten()[ithreshold]
da = ds['frequency_of_excedance'].isel(time=it, levels=ilev, thresholds=ithreshold)
# need to mask to lower countour level so that the lower limit of colorbar is correct (otherwise everything is lower-level color)
da_m = da.where(da > contour_levels[0])
cmap = colormaps[thresholds[ithreshold]]
norm = mpl.colors.BoundaryNorm(boundaries=contour_levels, ncolors=256)
da_m.plot(ax=ax, norm=norm, cmap=cmap, transform=ccrs.PlateCarree(), cbar_kwargs={'label': "Frequency [-]"})
# add volcano location
ax.plot(olon, olat, marker='^' , color='k', transform=ccrs.PlateCarree())
# zoom in to see the data
ax.set_extent(bnd)
# ax.set_extent([175, 180, -40, -35])
# add decorations
# ax.add_feature(cf.OCEAN)
ax.add_feature(cf.COASTLINE)
ax.add_feature(cf.BORDERS, lw=0.5)
ax.add_feature(cf.LAND, edgecolor='black', facecolor='wheat', linewidth=0.5, alpha=0.5)
ax.add_feature(cf.LAKES, edgecolor='black', facecolor=cf.COLORS['water']) #'none')
if clon == 180:
gl = ax.gridlines(xlocs=xlocs, draw_labels=True, ls=':', color='gray')
else:
gl = ax.gridlines(draw_labels=True, ls=':', color='gray')
gl.top_labels = False
gl.right_labels = False
ax.set_title(f'Frequency of excedance of {thresholds[ithreshold]} mg/m3')
fig.suptitle(tit_str)#, fontsize='small')
fig.tight_layout()
fig, axs = plt.subplots(len(thresholds), 1, figsize=(10,12), subplot_kw={"projection": ccrs.PlateCarree(clon)})
for ithreshold in range(len(thresholds)):
ax = axs.flatten()[ithreshold]
da = ds['frequency_of_excedance'].isel(time=it, levels=ilev, thresholds=ithreshold)
# need to mask to lower countour level so that the lower limit of colorbar is correct (otherwise everything is lower-level color)
da_m = da.where(da > contour_levels[0])
cmap = colormaps[thresholds[ithreshold]]
norm = mpl.colors.BoundaryNorm(boundaries=contour_levels, ncolors=256)
da_m.plot(ax=ax, norm=norm, cmap=cmap, transform=ccrs.PlateCarree(), cbar_kwargs={'label': "Frequency [-]"})
# add volcano location
ax.plot(olon, olat, marker='^' , color='k', transform=ccrs.PlateCarree())
# zoom in to see the data
ax.set_extent(bnd)
# ax.set_extent([175, 180, -40, -35])
# add decorations
# ax.add_feature(cf.OCEAN)
ax.add_feature(cf.COASTLINE)
ax.add_feature(cf.BORDERS, lw=0.5)
ax.add_feature(cf.LAND, edgecolor='black', facecolor='wheat', linewidth=0.5, alpha=0.5)
ax.add_feature(cf.LAKES, edgecolor='black', facecolor=cf.COLORS['water']) #'none')
if clon == 180:
gl = ax.gridlines(xlocs=xlocs, draw_labels=True, ls=':', color='gray')
else:
gl = ax.gridlines(draw_labels=True, ls=':', color='gray')
gl.top_labels = False
gl.right_labels = False
ax.set_title(f'Frequency of excedance of {thresholds[ithreshold]} mg/m3')
fig.suptitle(tit_str)#, fontsize='small')
fig.tight_layout()
In [ ]:
Copied!