tobac example: Tracking isolated convection based on updraft velocity and total condensate

This example notebook demonstrates the use of tobac to track isolated deep convective clouds in cloud-resolving model simulation output based on vertical velocity and total condensate mixing ratio.

The simulation results used in this example were performed as part of the ACPC deep convection intercomparison case study (http://acpcinitiative.org/Docs/ACPC_DCC_Roadmap_171019.pdf) with WRF using the Morrison microphysics scheme.

The data used in this example is downloaded from “zenodo link” automatically as part of the notebooks (This only has to be done once for all the tobac example notebooks).

Import libraries:

import iris
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import iris.plot as iplt
import iris.quickplot as qplt
import shutil
from six.moves import urllib
from pathlib import Path
%matplotlib inline

# Import tobac itself:
import tobac
print('using tobac version', str(tobac.__version__))

using tobac version 1.5.3

#Disable a couple of warnings:
import warnings
warnings.filterwarnings('ignore', category=UserWarning, append=True)
warnings.filterwarnings('ignore', category=RuntimeWarning, append=True)
warnings.filterwarnings('ignore', category=FutureWarning, append=True)
warnings.filterwarnings('ignore',category=pd.io.pytables.PerformanceWarning)

Download and load example data:

The actual dowloading is only necessary once for all example notebooks.

data_out=Path('../')

# Download the data: This only has to be done once for all tobac examples and can take a while
data_file = list(data_out.rglob('data/Example_input_midlevelUpdraft.nc'))
if len(data_file) == 0:
    file_path='https://zenodo.org/records/3195910/files/climate-processes/tobac_example_data-v1.0.1.zip'
    #file_path='http://zenodo..'
    tempfile=Path('temp.zip')
    print('start downloading data')
    request=urllib.request.urlretrieve(file_path, tempfile)
    print('start extracting data')
    shutil.unpack_archive(tempfile, data_out)
    tempfile.unlink()
    print('data extracted')
    data_file = list(data_out.rglob('data/Example_input_midlevelUpdraft.nc'))

Load Data from downloaded file:

data_file_W_mid_max = list(data_out.rglob('data/Example_input_midlevelUpdraft.nc'))[0]
data_file_TWC = list(data_out.rglob('data/Example_input_Condensate.nc'))[0]

W_mid_max=iris.load_cube(str(data_file_W_mid_max),'w')
TWC=iris.load_cube(str(data_file_TWC),'TWC')

# Display information about the two cubes for vertical velocity and total condensate mixing ratio:
display(W_mid_max)
display(TWC)

W (m s-1)	time	south_north	west_east
Shape	47	198	198
Dimension coordinates
time	x	-	-
south_north	-	x	-
west_east	-	-	x
Auxiliary coordinates
projection_y_coordinate	-	x	-
y	-	x	-
latitude	-	x	x
longitude	-	x	x
projection_x_coordinate	-	-	x
x	-	-	x
Scalar coordinates
bottom_top_stag	39, bound=(31, 47)
model_level_number	39, bound=(31, 47)
Cell methods
0	model_level_number: maximum
Attributes
Conventions	'CF-1.5'
FieldType	104
MemoryOrder	'XYZ'
description	'z-wind component'
stagger	'Z'

Twc (kg kg-1)	time	bottom_top	south_north	west_east
Shape	47	94	198	198
Dimension coordinates
time	x	-	-	-
bottom_top	-	x	-	-
south_north	-	-	x	-
west_east	-	-	-	x
Auxiliary coordinates
model_level_number	-	x	-	-
projection_y_coordinate	-	-	x	-
y	-	-	x	-
latitude	-	-	x	x
longitude	-	-	x	x
projection_x_coordinate	-	-	-	x
x	-	-	-	x
Attributes
Conventions	'CF-1.5'

#Set up directory to save output and plots:
savedir=Path("Save")
if not savedir.is_dir():
    savedir.mkdir()
plot_dir=Path("Plot")
if not plot_dir.is_dir():
    plot_dir.mkdir()

Feature detection:

Perform feature detection based on midlevel maximum vertical velocity and a range of threshold values.

# Determine temporal and spatial sampling of the input data:
dxy,dt=tobac.get_spacings(W_mid_max)

# Keyword arguments for feature detection step:
parameters_features={}
parameters_features['position_threshold']='weighted_diff'
parameters_features['sigma_threshold']=0.5
parameters_features['min_distance']=0
parameters_features['sigma_threshold']=1
parameters_features['threshold']=[3,5,10] #m/s
parameters_features['n_erosion_threshold']=0
parameters_features['n_min_threshold']=3

# Perform feature detection and save results:
print('start feature detection based on midlevel column maximum vertical velocity')
dxy,dt=tobac.get_spacings(W_mid_max)
Features=tobac.feature_detection_multithreshold(W_mid_max,dxy,**parameters_features)
print('feature detection performed start saving features')
Features.to_hdf(savedir / 'Features.h5', 'table')
print('features saved')

start feature detection based on midlevel column maximum vertical velocity

feature detection performed start saving features

features saved

Segmentation:

Perform segmentation based on 3D total condensate field to determine cloud volumes associated to identified features:

parameters_segmentation_TWC={}
parameters_segmentation_TWC['method']='watershed'
parameters_segmentation_TWC['threshold']=0.1e-3  # kg/kg mixing ratio

print('Start segmentation based on total water content')
Mask_TWC,Features_TWC=tobac.segmentation_3D(Features,TWC,dxy,**parameters_segmentation_TWC)
print('segmentation TWC performed, start saving results to files')
iris.save([Mask_TWC], savedir / 'Mask_Segmentation_TWC.nc',zlib=True,complevel=4)
Features_TWC.to_hdf(savedir / 'Features_TWC.h5','table')
print('segmentation TWC performed and saved')

Start segmentation based on total water content

segmentation TWC performed, start saving results to files

segmentation TWC performed and saved

Trajectory linking:

Detected features are linked into cloud trajectories using the trackpy library (http://soft-matter.github.io/trackpy). This takes the feature positions determined in the feature detection step into account but does not include information on the shape of the identified objects.

# Keyword arguments for linking step:
parameters_linking={}
parameters_linking['method_linking']='predict'
parameters_linking['adaptive_stop']=0.2
parameters_linking['adaptive_step']=0.95
parameters_linking['extrapolate']=0
parameters_linking['order']=1
parameters_linking['subnetwork_size']=100
parameters_linking['memory']=0
parameters_linking['time_cell_min']=5*60
parameters_linking['method_linking']='predict'
parameters_linking['v_max']=10

# Perform linking and save results:
Track=tobac.linking_trackpy(Features,W_mid_max,dt=dt,dxy=dxy,**parameters_linking)
Track.to_hdf(savedir / 'Track.h5', 'table')

Frame 46: 18 trajectories present.

Visualisation:

# Set extent for maps plotted in the following cells ( in the form [lon_min,lon_max,lat_min,lat_max])
axis_extent=[-95,-93.8,29.5,30.6]

# Plot map with all individual tracks:
import cartopy.crs as ccrs
fig_map,ax_map=plt.subplots(figsize=(10,10),subplot_kw={'projection': ccrs.PlateCarree()})
ax_map=tobac.map_tracks(Track,axis_extent=axis_extent,axes=ax_map)

# Create animation showing tracked cells with outline of cloud volumes and the midlevel vertical velocity as a background field:
animation_tobac=tobac.animation_mask_field(track=Track,features=Features,field=W_mid_max,mask=Mask_TWC,
                                          axis_extent=axis_extent,#figsize=figsize,orientation_colorbar='horizontal',pad_colorbar=0.2,
                                          vmin=0,vmax=20,extend='both',cmap='Blues',
                                          interval=500,figsize=(10,7),
                                          plot_outline=True,plot_marker=True,marker_track='x',plot_number=True,plot_features=True)

# Display animation:
from IPython.display import HTML, Image, display
HTML(animation_tobac.to_html5_video())

<Figure size 640x480 with 0 Axes>

# # Save animation to file
# savefile_animation=plot_dir / 'Animation.mp4')
# animation_tobac.save(savefile_animation,dpi=200)
# print(f'animation saved to {savefile_animation}')

# Updraft lifetimes of tracked cells:
fig_lifetime,ax_lifetime=plt.subplots()
tobac.plot_lifetime_histogram_bar(Track,axes=ax_lifetime,bin_edges=np.arange(0,120,10),density=False,width_bar=8)
ax_lifetime.set_xlabel('lifetime (min)')
ax_lifetime.set_ylabel('counts')

Text(0, 0.5, 'counts')