tobac example: Tracking deep convection based on OLR from geostationary satellite retrievals
This example notebook demonstrates the use of tobac to track isolated deep convective clouds based on outgoing longwave radiation (OLR) calculated based on a combination of two different channels of the GOES-13 imaging instrument.
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).
[3]:
# Import libraries:
import xarray
import numpy as np
import pandas as pd
import os
from six.moves import urllib
from glob import glob
import matplotlib.pyplot as plt
%matplotlib inline
[1]:
# 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:
This has to be done only once for all tobac examples.
[5]:
data_out='../'
[5]:
# # 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'
# tempfile='temp.zip'
# print('start downloading data')
# request=urllib.request.urlretrieve(file_path,tempfile)
# print('start extracting data')
# zf = zipfile.ZipFile(tempfile)
# zf.extractall(data_out)
# print('example data saved in')
Load data:
[6]:
data_file=os.path.join(data_out,'*','data','Example_input_OLR_satellite.nc')
data_file = glob(data_file)[0]
[10]:
# Load Data from downloaded file:
OLR=xarray.open_dataset(data_file)['olr']
[9]:
# Display information about the input data cube:
display(OLR)
| Olr (W m^-2) | time | latitude | longitude |
|---|---|---|---|
| Shape | 54 | 131 | 184 |
| Dimension coordinates | |||
| time | x | - | - |
| latitude | - | x | - |
| longitude | - | - | x |
| Attributes | |||
| Conventions | CF-1.5 | ||
[12]:
#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 identification:
Identify features based on OLR field and a set of threshold values
[13]:
# Determine temporal and spatial sampling of the input data:
dxy,dt=tobac.utils.get_spacings(OLR,grid_spacing=4000)
(<xarray.DataArray 'olr' (time: 54, lat: 131, lon: 184)>
[1301616 values with dtype=float64]
Coordinates:
* time (time) datetime64[ns] 2013-06-19T19:02:22 ... 2013-06-20T02:45:18
* lat (lat) float64 28.03 28.07 28.11 28.15 ... 32.88 32.92 32.96 32.99
* lon (lon) float64 -94.99 -94.95 -94.91 -94.87 ... -88.08 -88.04 -88.01
Attributes:
long_name: OLR
units: W m^-2,)
{'grid_spacing': 4000}
converting xarray to iris and back
(<iris 'Cube' of OLR / (W m^-2) (time: 54; latitude: 131; longitude: 184)>,)
{'grid_spacing': 4000}
(4000, 476)
[14]:
# Keyword arguments for the feature detection step
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]:
# Feature detection and save results to file:
print('starting feature detection')
Features=tobac.themes.tobac_v1.feature_detection_multithreshold(OLR,dxy,**parameters_features)
Features.to_netcdf(os.path.join(savedir,'Features.nc'))
print('feature detection performed and saved')
starting feature detection
( frame idx hdim_1 hdim_2 num threshold_value feature
0 0 2 0.737149 177.294423 4 250 1
1 0 3 3.354292 162.980569 9 250 2
2 0 5 32.000000 138.000000 1 250 3
3 0 8 37.479003 146.971379 6 250 4
4 0 13 65.790337 2.067482 3 250 5
... ... ... ... ... ... ... ...
2733 53 24 64.980999 91.292535 16 225 2734
2734 53 25 83.297251 157.068621 89 225 2735
2735 53 26 103.000000 1.000000 1 225 2736
2736 53 27 129.853687 11.461736 4 225 2737
2737 53 29 46.832556 29.433805 24 200 2738
[2738 rows x 7 columns], <iris 'Cube' of OLR / (W m^-2) (time: 54; latitude: 131; longitude: 184)>)
{}
frame idx hdim_1 hdim_2 num threshold_value feature \
0 0 2 0.737149 177.294423 4 250 1
1 0 3 3.354292 162.980569 9 250 2
2 0 5 32.000000 138.000000 1 250 3
3 0 8 37.479003 146.971379 6 250 4
4 0 13 65.790337 2.067482 3 250 5
... ... ... ... ... ... ... ...
2733 53 24 64.980999 91.292535 16 225 2734
2734 53 25 83.297251 157.068621 89 225 2735
2735 53 26 103.000000 1.000000 1 225 2736
2736 53 27 129.853687 11.461736 4 225 2737
2737 53 29 46.832556 29.433805 24 200 2738
time timestr latitude longitude
0 2013-06-19 19:02:22 2013-06-19 19:02:22 28.060909 -88.223109
1 2013-06-19 19:02:22 2013-06-19 19:02:22 28.160780 -88.769332
2 2013-06-19 19:02:22 2013-06-19 19:02:22 29.253914 -89.722603
3 2013-06-19 19:02:22 2013-06-19 19:02:22 29.462995 -89.380251
4 2013-06-19 19:02:22 2013-06-19 19:02:22 30.543369 -94.909851
... ... ... ... ...
2733 2013-06-20 02:45:18 2013-06-20 02:45:18 30.512484 -91.504982
2734 2013-06-20 02:45:18 2013-06-20 02:45:18 31.211441 -88.994935
2735 2013-06-20 02:45:18 2013-06-20 02:45:18 31.963307 -94.950587
2736 2013-06-20 02:45:18 2013-06-20 02:45:18 32.988057 -94.551362
2737 2013-06-20 02:45:18 2013-06-20 02:45:18 29.819931 -93.865540
[2738 rows x 11 columns]
frame idx hdim_1 hdim_2 num threshold_value feature \
0 0 2 0.737149 177.294423 4 250 1
1 0 3 3.354292 162.980569 9 250 2
2 0 5 32.000000 138.000000 1 250 3
3 0 8 37.479003 146.971379 6 250 4
4 0 13 65.790337 2.067482 3 250 5
... ... ... ... ... ... ... ...
2733 53 24 64.980999 91.292535 16 225 2734
2734 53 25 83.297251 157.068621 89 225 2735
2735 53 26 103.000000 1.000000 1 225 2736
2736 53 27 129.853687 11.461736 4 225 2737
2737 53 29 46.832556 29.433805 24 200 2738
time timestr latitude longitude
0 2013-06-19 19:02:22 2013-06-19 19:02:22 28.060909 -88.223109
1 2013-06-19 19:02:22 2013-06-19 19:02:22 28.160780 -88.769332
2 2013-06-19 19:02:22 2013-06-19 19:02:22 29.253914 -89.722603
3 2013-06-19 19:02:22 2013-06-19 19:02:22 29.462995 -89.380251
4 2013-06-19 19:02:22 2013-06-19 19:02:22 30.543369 -94.909851
... ... ... ... ...
2733 2013-06-20 02:45:18 2013-06-20 02:45:18 30.512484 -91.504982
2734 2013-06-20 02:45:18 2013-06-20 02:45:18 31.211441 -88.994935
2735 2013-06-20 02:45:18 2013-06-20 02:45:18 31.963307 -94.950587
2736 2013-06-20 02:45:18 2013-06-20 02:45:18 32.988057 -94.551362
2737 2013-06-20 02:45:18 2013-06-20 02:45:18 29.819931 -93.865540
[2738 rows x 11 columns]
['time', 'lat', 'lon']
feature detection performed and saved
Segmentation:
Segmentation is performed based on the OLR field and a threshold value to determine the cloud areas.
[16]:
# Keyword arguments for the segmentation step:
parameters_segmentation={}
parameters_segmentation['target']='minimum'
parameters_segmentation['method']='watershed'
parameters_segmentation['threshold']=250
[18]:
# Perform segmentation and save results to files:
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')
<xarray.DataArray 'olr' (time: 54, lat: 131, lon: 184)>
array([[[306.847749, 306.847749, ..., 263.987914, 279.675285],
[306.847749, 306.847749, ..., 284.331736, 285.586608],
...,
[290.318092, 290.318092, ..., 294.654012, 294.654012],
[290.318092, 290.318092, ..., 294.654012, 297.020815]],
[[306.847749, 304.470854, ..., 272.836495, 272.836495],
[304.470854, 304.470854, ..., 281.993749, 285.586608],
...,
[290.318092, 290.318092, ..., 297.020815, 294.654012],
[290.318092, 290.318092, ..., 297.020815, 294.654012]],
...,
[[305.300037, 305.300037, ..., 294.250334, 294.250334],
[306.847749, 306.847749, ..., 294.250334, 294.250334],
...,
[261.385907, 248.390063, ..., 290.903128, 288.574062],
[258.106908, 258.106908, ..., 290.903128, 288.574062]],
[[306.847749, 306.847749, ..., 294.250334, 294.250334],
[306.847749, 306.847749, ..., 294.250334, 294.250334],
...,
[260.32182 , 254.87151 , ..., 290.903128, 288.574062],
[258.106908, 255.910795, ..., 290.903128, 288.574062]]])
Coordinates:
* time (time) datetime64[ns] 2013-06-19T19:02:22 ... 2013-06-20T02:45:18
* lat (lat) float64 28.03 28.07 28.11 28.15 ... 32.88 32.92 32.96 32.99
* lon (lon) float64 -94.99 -94.95 -94.91 -94.87 ... -88.08 -88.04 -88.01
Attributes:
long_name: OLR
units: W m^-2
segmentation OLR performed, start saving results to files
segmentation OLR performed and saved
Trajectory linking:
The 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.
[19]:
# keyword arguments for linking step
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 53: 19 trajectories present.
Visualisation:
[19]:
# Set extent of maps created in the following cells:
axis_extent=[-95,-89,28,32]
[20]:
# 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)
[21]:
# 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,cmap='Blues_r',
plot_outline=True,plot_marker=True,marker_track='x',plot_number=True,plot_features=True)
[23]:
# Display animation:
from IPython.display import HTML, Image, display
HTML(animation_test_tobac.to_html5_video())
[23]:
[ ]:
# # 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}')
[22]:
# 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')
[22]:
Text(0, 0.5, 'counts')
[ ]: