tobac example: Tracking of deep convection based on OLR from convection permitting model simulations
This example notebook demonstrates the use of tobac to track deep convection based on the outgoing longwave radiation (OLR) from convection permitting simulations.
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. Simulations were performed with a horizontal grid spacing of 4.5 km.
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:
[9]:
# Import a range of python libraries used in this notebook:
import xarray
import numpy as np
import pandas as pd
import os,sys
import shutil
import datetime
from six.moves import urllib
from glob import glob
import matplotlib.pyplot as plt
%matplotlib inline
[2]:
# Import tobac itself:
import tobac
[4]:
# Disable a few 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 example data:
The actual download is only necessary once for all example notebooks.
[5]:
data_out='../'
[7]:
# # Download the data: This only has to be done once for all tobac examples and can take a while
# file_path='https://zenodo.org/record/3195910/files/climate-processes/tobac_example_data-v1.0.1.zip'
# #file_path='http://zenodo..'
# tempfile='temp.zip'
# print('start downloading data')
# request=urllib.request.urlretrieve(file_path,tempfile)
# print('start extracting data')
# shutil.unpack_archive(tempfile,data_out)
# #zf = zipfile.ZipFile(tempfile)
# #zf.extractall(data_out)
# os.remove(tempfile)
# print('data extracted')
start downloading data
start extracting data
data extracted
[6]:
data_file=os.path.join(data_out,'*','data','Example_input_OLR_model.nc')
data_file = glob(data_file)[0]
[10]:
#Load Data from downloaded file:
OLR=xarray.open_dataset(data_file)['OLR']
[11]:
#Set up directory to save output and plots:
savedir='Save'
if not os.path.exists(savedir):
os.makedirs(savedir)
plot_dir="Plot"
if not os.path.exists(plot_dir):
os.makedirs(plot_dir)
Feature detection:
Feature detection is performed based on OLR field and a set of thresholds.
[12]:
# Determine temporal and spatial sampling:
dxy,dt=tobac.utils.get_spacings(OLR)
(<xarray.DataArray 'OLR' (time: 96, south_north: 110, west_east: 132)>
[1393920 values with dtype=float32]
Coordinates:
* time (time) datetime64[ns] 2013-06-19T19:05:00 ... 2013-06-20T03:00:00
* south_north (south_north) int64 156 157 158 159 160 ... 261 262 263 264 265
* west_east (west_east) int64 201 202 203 204 205 ... 328 329 330 331 332
latitude (south_north, west_east) float32 ...
longitude (south_north, west_east) float32 ...
x (west_east) float64 ...
y (south_north) float64 ...
x_0 (west_east) int64 ...
y_0 (south_north) int64 ...
Attributes:
units: W m-2,)
{}
converting xarray to iris and back
(<iris 'Cube' of OLR / (W m-2) (time: 96; south_north: 110; west_east: 132)>,)
{}
(4500.0, 300)
[14]:
# Dictionary containing keyword arguments for feature detection step (Keywords could also be given directly in the function call).
parameters_features={}
parameters_features['position_threshold']='weighted_diff'
parameters_features['sigma_threshold']=0.5
parameters_features['min_num']=4
parameters_features['target']='minimum'
parameters_features['threshold']=[250,225,200,175,150]
[15]:
# Perform feature detection:
print('starting feature detection')
Features=tobac.themes.tobac_v1.feature_detection_multithreshold(OLR,dxy,**parameters_features)
Features.to_netcdf(os.path.join(savedir,'Features.netcdf'))
print('feature detection performed and saved')
starting feature detection
( frame idx hdim_1 hdim_2 num threshold_value feature
0 0 1 48.000000 104.999952 3 250 1
1 0 2 49.000000 100.000000 1 250 2
2 0 4 59.462142 84.000000 2 250 3
3 0 7 65.000000 67.000000 1 250 4
4 0 9 67.000000 62.346732 2 250 5
... ... ... ... ... ... ... ...
2716 95 28 67.000000 110.000000 1 175 2717
2717 95 29 71.666555 122.473032 62 175 2718
2718 95 30 69.802409 15.037454 10 175 2719
2719 95 31 75.380948 109.634011 19 175 2720
2720 95 34 88.139481 113.625911 68 150 2721
[2721 rows x 7 columns], <iris 'Cube' of OLR / (W m-2) (time: 96; south_north: 110; west_east: 132)>)
{}
frame idx hdim_1 hdim_2 num threshold_value feature \
0 0 1 48.000000 104.999952 3 250 1
1 0 2 49.000000 100.000000 1 250 2
2 0 4 59.462142 84.000000 2 250 3
3 0 7 65.000000 67.000000 1 250 4
4 0 9 67.000000 62.346732 2 250 5
... ... ... ... ... ... ... ...
2716 95 28 67.000000 110.000000 1 175 2717
2717 95 29 71.666555 122.473032 62 175 2718
2718 95 30 69.802409 15.037454 10 175 2719
2719 95 31 75.380948 109.634011 19 175 2720
2720 95 34 88.139481 113.625911 68 150 2721
time timestr south_north west_east \
0 2013-06-19 19:05:00 2013-06-19 19:05:00 204.000000 305.999952
1 2013-06-19 19:05:00 2013-06-19 19:05:00 205.000000 301.000000
2 2013-06-19 19:05:00 2013-06-19 19:05:00 215.462142 285.000000
3 2013-06-19 19:05:00 2013-06-19 19:05:00 221.000000 268.000000
4 2013-06-19 19:05:00 2013-06-19 19:05:00 223.000000 263.346732
... ... ... ... ...
2716 2013-06-20 03:00:00 2013-06-20 03:00:00 223.000000 311.000000
2717 2013-06-20 03:00:00 2013-06-20 03:00:00 227.666555 323.473032
2718 2013-06-20 03:00:00 2013-06-20 03:00:00 225.802409 216.037454
2719 2013-06-20 03:00:00 2013-06-20 03:00:00 231.380948 310.634011
2720 2013-06-20 03:00:00 2013-06-20 03:00:00 244.139481 314.625911
projection_y_coordinate y latitude \
0 9.202500e+05 204.000000 [29.56323250795596]
1 9.247500e+05 205.000000 [29.61300277709961]
2 9.718296e+05 215.462142 [30.06779582764317]
3 9.967500e+05 221.000000 [30.317230224609375]
4 1.005750e+06 223.000000 [30.40456474862914]
... ... ... ...
2716 1.005750e+06 223.000000 [30.334190368652344]
2717 1.026749e+06 227.666555 [30.500885027226275]
2718 1.018361e+06 225.802409 [30.5510583635925]
2719 1.043464e+06 231.380948 [30.678966958276725]
2720 1.100878e+06 244.139481 [31.19504789023921]
longitude projection_x_coordinate x
0 [-90.0455955414817] 1.379250e+06 305.999952
1 [-90.2796630859375] 1.356750e+06 301.000000
2 [-91.01834241734284] 1.284750e+06 285.000000
3 [-91.81805419921875] 1.208250e+06 268.000000
4 [-92.03704073277446] 1.187310e+06 263.346732
... ... ... ...
2716 [-89.76840209960938] 1.401750e+06 311.000000
2717 [-89.16383151140309] 1.457879e+06 323.473032
2718 [-94.29081414933923] 9.744185e+05 216.037454
2719 [-89.76749015468663] 1.400103e+06 310.634011
2720 [-89.54780638125565] 1.418067e+06 314.625911
[2721 rows x 17 columns]
frame idx hdim_1 hdim_2 num threshold_value feature \
0 0 1 48.000000 104.999952 3 250 1
1 0 2 49.000000 100.000000 1 250 2
2 0 4 59.462142 84.000000 2 250 3
3 0 7 65.000000 67.000000 1 250 4
4 0 9 67.000000 62.346732 2 250 5
... ... ... ... ... ... ... ...
2716 95 28 67.000000 110.000000 1 175 2717
2717 95 29 71.666555 122.473032 62 175 2718
2718 95 30 69.802409 15.037454 10 175 2719
2719 95 31 75.380948 109.634011 19 175 2720
2720 95 34 88.139481 113.625911 68 150 2721
time timestr south_north west_east \
0 2013-06-19 19:05:00 2013-06-19 19:05:00 204.000000 305.999952
1 2013-06-19 19:05:00 2013-06-19 19:05:00 205.000000 301.000000
2 2013-06-19 19:05:00 2013-06-19 19:05:00 215.462142 285.000000
3 2013-06-19 19:05:00 2013-06-19 19:05:00 221.000000 268.000000
4 2013-06-19 19:05:00 2013-06-19 19:05:00 223.000000 263.346732
... ... ... ... ...
2716 2013-06-20 03:00:00 2013-06-20 03:00:00 223.000000 311.000000
2717 2013-06-20 03:00:00 2013-06-20 03:00:00 227.666555 323.473032
2718 2013-06-20 03:00:00 2013-06-20 03:00:00 225.802409 216.037454
2719 2013-06-20 03:00:00 2013-06-20 03:00:00 231.380948 310.634011
2720 2013-06-20 03:00:00 2013-06-20 03:00:00 244.139481 314.625911
projection_y_coordinate y latitude \
0 9.202500e+05 204.000000 [29.56323250795596]
1 9.247500e+05 205.000000 [29.61300277709961]
2 9.718296e+05 215.462142 [30.06779582764317]
3 9.967500e+05 221.000000 [30.317230224609375]
4 1.005750e+06 223.000000 [30.40456474862914]
... ... ... ...
2716 1.005750e+06 223.000000 [30.334190368652344]
2717 1.026749e+06 227.666555 [30.500885027226275]
2718 1.018361e+06 225.802409 [30.5510583635925]
2719 1.043464e+06 231.380948 [30.678966958276725]
2720 1.100878e+06 244.139481 [31.19504789023921]
longitude projection_x_coordinate x
0 [-90.0455955414817] 1.379250e+06 305.999952
1 [-90.2796630859375] 1.356750e+06 301.000000
2 [-91.01834241734284] 1.284750e+06 285.000000
3 [-91.81805419921875] 1.208250e+06 268.000000
4 [-92.03704073277446] 1.187310e+06 263.346732
... ... ... ...
2716 [-89.76840209960938] 1.401750e+06 311.000000
2717 [-89.16383151140309] 1.457879e+06 323.473032
2718 [-94.29081414933923] 9.744185e+05 216.037454
2719 [-89.76749015468663] 1.400103e+06 310.634011
2720 [-89.54780638125565] 1.418067e+06 314.625911
[2721 rows x 17 columns]
['time', 'south_north', 'west_east', 'latitude', 'longitude', 'x', 'y', 'x_0', 'y_0']
feature detection performed and saved
Segmentation:
Segmentation is performed with watershedding based on the detected features and a single threshold value.
[16]:
# Dictionary containing keyword options for the segmentation step:
parameters_segmentation={}
parameters_segmentation['target']='minimum'
parameters_segmentation['method']='watershed'
parameters_segmentation['threshold']=250
[18]:
# Perform segmentation and save results:
print('Starting segmentation based on OLR.')
Mask_OLR,Features_OLR=tobac.themes.tobac_v1.segmentation(Features,OLR,dxy,**parameters_segmentation)
print('segmentation OLR performed, start saving results to files')
Mask_OLR.to_netcdf(os.path.join(savedir,'Mask_Segmentation_OLR.nc'))
Features_OLR.to_netcdf(os.path.join(savedir,'Features_OLR.nc'))
print('segmentation OLR performed and saved')
Starting segmentation based on OLR.
<xarray.DataArray 'OLR' (time: 96, south_north: 110, west_east: 132)>
array([[[289.05252, 288.72162, ..., 280.26904, 280.21326],
[289.42636, 289.08936, ..., 280.36798, 280.3152 ],
...,
[289.93964, 295.0356 , ..., 283.53815, 274.44553],
[296.27737, 296.89362, ..., 279.30847, 274.90115]],
[[289.3221 , 288.9656 , ..., 280.25012, 280.19217],
[289.7054 , 289.34625, ..., 280.35687, 280.3002 ],
...,
[295.6821 , 286.90628, ..., 286.26834, 275.19684],
[296.15982, 296.80145, ..., 281.44797, 275.65378]],
...,
[[300.36115, 300.2792 , ..., 278.89587, 278.8872 ],
[300.26105, 300.16174, ..., 278.63638, 278.67062],
...,
[296.05872, 296.319 , ..., 288.2268 , 288.3182 ],
[295.82617, 296.04742, ..., 287.981 , 287.89703]],
[[300.3252 , 300.24365, ..., 278.6296 , 278.71307],
[300.2927 , 300.2258 , ..., 278.47296, 278.5284 ],
...,
[295.60107, 296.1242 , ..., 288.0909 , 288.09134],
[295.24384, 295.66608, ..., 287.80396, 287.7296 ]]], dtype=float32)
Coordinates:
* time (time) datetime64[ns] 2013-06-19T19:05:00 ... 2013-06-20T03:00:00
* south_north (south_north) int64 156 157 158 159 160 ... 261 262 263 264 265
* west_east (west_east) int64 201 202 203 204 205 ... 328 329 330 331 332
latitude (south_north, west_east) float32 27.684788 ... 32.014114
longitude (south_north, west_east) float32 -95.00824 ... -88.65869
x (west_east) float64 9.068e+05 9.112e+05 ... 1.492e+06 1.496e+06
y (south_north) float64 7.042e+05 7.088e+05 ... 1.195e+06
x_0 (west_east) int64 201 202 203 204 205 ... 328 329 330 331 332
y_0 (south_north) int64 156 157 158 159 160 ... 261 262 263 264 265
Attributes:
units: W m-2
segmentation OLR performed, start saving results to files
segmentation OLR performed and saved
Trajectory linking:
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.**
[19]:
# Arguments for trajectory linking:
parameters_linking={}
parameters_linking['v_max']=20
parameters_linking['stubs']=2
parameters_linking['order']=1
parameters_linking['extrapolate']=1
parameters_linking['memory']=0
parameters_linking['adaptive_stop']=0.2
parameters_linking['adaptive_step']=0.95
parameters_linking['subnetwork_size']=100
parameters_linking['method_linking']= 'predict'
[20]:
# Perform linking and save results to file:
Track=tobac.themes.tobac_v1.linking_trackpy(Features,OLR,dt=dt,dxy=dxy,**parameters_linking)
Track.to_netcdf(os.path.join(savedir,'Track.nc'))
Frame 95: 18 trajectories present.
Visualisation:
[17]:
# Set extent of maps created in the following cells:
axis_extent=[-95,-89,28,32]
[18]:
# 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.plot.map_tracks(Track,axis_extent=axis_extent,axes=ax_map)
[19]:
# Create animation of tracked clouds and outlines with OLR as a background field
animation_test_tobac=tobac.plot.animation_mask_field(Track,Features,OLR,Mask_OLR,
axis_extent=axis_extent,#figsize=figsize,orientation_colorbar='horizontal',pad_colorbar=0.2,
vmin=80,vmax=330,
plot_outline=True,plot_marker=True,marker_track='x',plot_number=True,plot_features=True)
[20]:
# Display animation:
from IPython.display import HTML, Image, display
HTML(animation_test_tobac.to_html5_video())
[20]:
[ ]:
# # Save animation to file:
# savefile_animation=os.path.join(plot_dir,'Animation.mp4')
# animation_test_tobac.save(savefile_animation,dpi=200)
# print(f'animation saved to {savefile_animation}')
[21]:
# Lifetimes of tracked clouds:
fig_lifetime,ax_lifetime=plt.subplots()
tobac.plot.plot_lifetime_histogram_bar(Track,axes=ax_lifetime,bin_edges=np.arange(0,200,20),density=False,width_bar=10)
ax_lifetime.set_xlabel('lifetime (min)')
ax_lifetime.set_ylabel('counts')
[21]:
Text(0, 0.5, 'counts')
[ ]: