tobac example: Tracking deep convection based on VIS from geostationary satellite retrievals
This example notebook demonstrates the use of tobac to track isolated deep convective clouds based on radiances within the VIS using channel 2 (red light - 600 nm) of the GOES-16 imaging instrument in 5-min resolution. The study area is loacted within the CONUS extent of the GOES-E for investigating the formation of deep convection over the Carribean, following the EUREC4A initiative (http://eurec4a.eu/).
The data used in this example is saved on the cloud of the Amazon Web Services, providing an efficient way of processing satellite data without facing the need of downloading the data.
In this example, the Cloud and Moisture Imagery data (ABI-L2-CMIPC) data set is used for the cloud tracking. The product contains one or more Earth-view images with pixel values identifying brightness values that are scaled to support visual analysis. A resolution of 500 m in the VIS channels of the imager ensures a majority of features to be detectable. Also, the data includes quality information as well as information on the solar zenith angle.
Further information on the AWS and provided data sets can be found here.
Accessing GOES data
Configurations for conducting the example:
[14]:
# Import libraries:
import requests
import netCDF4
import boto3
from botocore import UNSIGNED
from botocore.config import Config
import xarray
import rioxarray
import rasterio
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
[2]:
# Import tobac itself:
import tobac
[3]:
# 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)
Build access to Amazon Web Server where GOES-16 data are stored:
[4]:
# For accessing data from AWS bucket first define specifics:
bucket_name = 'noaa-goes16'
product_name = 'ABI-L2-CMIPC'
year = 2020
day_of_year = 45
hour = 18
band = 2
[5]:
# Initialize an s3 client:
s3_client = boto3.client('s3', config=Config(signature_version=UNSIGNED))
[6]:
# Function for generating file name keys for the S3 bucket:
def get_s3_keys(bucket, s3_client, prefix = ''):
"""
Generate the keys in an S3 bucket.
:param bucket: Name of the S3 bucket.
:param prefix: Only fetch keys that start with this prefix (optional).
"""
kwargs = {'Bucket': bucket}
if isinstance(prefix, str):
kwargs['Prefix'] = prefix
while True:
resp = s3_client.list_objects_v2(**kwargs)
for obj in resp['Contents']:
key = obj['Key']
if key.startswith(prefix):
yield key
try:
kwargs['ContinuationToken'] = resp['NextContinuationToken']
except KeyError:
break
[7]:
# Retrieve the keys for of the file names:
keys = get_s3_keys(bucket_name,
s3_client,
prefix = f'{product_name}/{year}/{day_of_year:03.0f}/{hour:02.0f}/OR_{product_name}-M6C{band:02.0f}'
)
keys = [key for key in keys][0:6]
Request data from AWS S3 server using the keys and store all files to one data set:
[8]:
for k in range(len(keys)):
resp = requests.get(f'https://{bucket_name}.s3.amazonaws.com/{keys[k]}')
file_name = keys[k].split('/')[-1].split('.')[0]
nc4_ds = netCDF4.Dataset(file_name, memory = resp.content)
store = xarray.backends.NetCDF4DataStore(nc4_ds)
if k == 0:
DS = xarray.open_dataset(store)
else:
DS2 = xarray.open_dataset(store)
DS = xarray.combine_nested([DS, DS2], concat_dim=["t"], combine_attrs = "override")
Processing GOES data
Reprojecting image coordinates:
[9]:
# Extract properties to define a projection string:
geos_p = DS['goes_imager_projection']
proj_name = DS.goes_imager_projection.grid_mapping_name[0]
h = DS.goes_imager_projection.perspective_point_height
a = DS.goes_imager_projection.semi_major_axis
b = DS.goes_imager_projection.semi_minor_axis
rf = DS.goes_imager_projection.inverse_flattening
lat_0 = DS.goes_imager_projection.latitude_of_projection_origin
lon_0 = DS.goes_imager_projection.longitude_of_projection_origin
sweep = DS.goes_imager_projection.sweep_angle_axis
pstr = '+proj=geos +h=%f +a=%f +b=%f +rf=%f +lat_0=%f +lon_0=%f' \
% (h, a, b, rf, lat_0, lon_0)
[10]:
DS = DS.rio.write_crs(pstr, inplace=True)
[11]:
# Multiply the original coordinates by the satellite height:
DS["x"] = DS["x"]*h
DS["y"] = DS["y"]*h
[12]:
# Reduce dataset to only variables needed for further analysis:
DS = DS[["CMI"]]
[13]:
# Reproject to WGS84 (EPSG: 4326) using rasterio:
DS_R = DS.rio.reproject("epsg:4326")
Plot CONUS extent:
[16]:
from mpl_toolkits.basemap import Basemap
[17]:
fig = plt.figure(figsize=(15,15))
ax_0 = fig.add_subplot(1, 1, 1)
cmap = Basemap(projection='cyl',llcrnrlon=-153,llcrnrlat=14,urcrnrlon=-53,urcrnrlat=57, resolution='i', ax=ax_0)
con = cmap.pcolormesh(DS_R.x,DS_R.y,DS_R.CMI[0], cmap='cividis', vmin=0,vmax=1)
cmap.drawcountries(color='white', linewidth=.5, ax=ax_0)
cmap.drawcoastlines(linewidth=.5, color='white', ax=ax_0)
cmap.drawparallels(np.arange(10,60,10),labels=[1,0,0,0], linewidth=.5, color="grey")
cmap.drawmeridians(np.arange(-160,-50,10),labels=[0,0,0,1], linewidth=.5, color="grey")
cmap.colorbar(con, shrink=0.6, label='Reflectance', ax=ax_0)
[17]:
<matplotlib.colorbar.Colorbar at 0x7fa55d020670>
Crop the image to a smaller extent, here focussing on the isle of Jamaica:
[18]:
# Define boundaries of mask used for cropping:
min_lon = -80.0
min_lat = 15.0
max_lon = -75.0
max_lat = 20.0
mask_lon = (DS_R.x >= min_lon) & (DS_R.x <= max_lon)
mask_lat = (DS_R.y >= min_lat) & (DS_R.y <= max_lat)
[19]:
cropped_ds = DS_R.where(mask_lon & mask_lat, drop=True)
Plot cropped image showing Jamaica:
[20]:
fig = plt.figure(figsize=(15,15))
ax_1 = fig.add_subplot(1, 1, 1)
bmap = Basemap(projection='cyl',llcrnrlon=-79,llcrnrlat=17.5,urcrnrlon=-76,urcrnrlat=19, resolution='i', ax=ax_1)
dat = bmap.pcolormesh(cropped_ds.x,cropped_ds.y,cropped_ds.CMI[0], cmap='cividis', vmin=0,vmax=1)
bmap.drawcountries(color='white', linewidth=1, ax=ax_1)
bmap.drawcoastlines(linewidth=1, color='white', ax=ax_1)
bmap.drawparallels(np.arange(10,20,1),labels=[1,0,0,0], linewidth=1, color="grey")
bmap.drawmeridians(np.arange(-80,-70,1),labels=[0,0,0,1], linewidth=1, color="grey")
bmap.colorbar(dat, ax=ax_1, shrink=0.6, label='Reflectance')
[20]:
<matplotlib.colorbar.Colorbar at 0x7fa55cd97910>
Tobac Cloud Tracking
[21]:
# 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)
[22]:
# Change standard_name to be a regular argument of cf-conventions:
cropped_ds.CMI.attrs["standard_name"] = "toa_bidirectional_reflectance"
[23]:
# Assign former coordinates as variables to avoid dupplicate naming:
new_data = cropped_ds.reset_coords(names = ["y_image", "x_image"])
Feature identification: Identify features based on VIS field and a set of threshold values.
[24]:
# Determine temporal and spatial sampling of the input data:
dxy,dt=tobac.utils.get_spacings(new_data.CMI,grid_spacing=500)
[26]:
# 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']=8
parameters_features['n_min_threshold']=8
parameters_features['target']='maximum'
parameters_features['threshold']=[0.2,0.4,0.6]
[27]:
# Feature detection and save results to file:
print('starting feature detection')
Features=tobac.themes.tobac_v1.feature_detection_multithreshold(new_data.CMI,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 6 3.967704 199.359137 13 0.2 1
1 0 9 1.384956 308.325801 11 0.2 2
2 0 11 0.973372 340.422488 25 0.2 3
3 0 14 1.055909 479.695983 18 0.2 4
4 0 27 5.538348 3.698972 32 0.2 5
... ... ... ... ... ... ... ...
1487 5 1966 335.889404 285.967179 9 0.6 1488
1488 5 1967 422.424521 34.334552 24 0.6 1489
1489 5 1975 472.002232 31.821198 16 0.6 1490
1490 5 1986 522.294278 142.710449 15 0.6 1491
1491 5 2002 538.749685 148.998903 18 0.6 1492
time timestr latitude longitude \
0 2020-02-14 18:02:17 2020-02-14 18:02:17 19.958121 -78.153419
1 2020-02-14 18:02:17 2020-02-14 18:02:17 19.981960 -77.147658
2 2020-02-14 18:02:17 2020-02-14 18:02:17 19.985759 -76.851406
3 2020-02-14 18:02:17 2020-02-14 18:02:17 19.984997 -75.565914
4 2020-02-14 18:02:17 2020-02-14 18:02:17 19.943624 -79.959360
... ... ... ... ...
1487 2020-02-14 18:27:17 2020-02-14 18:27:17 16.894488 -77.354028
1488 2020-02-14 18:27:17 2020-02-14 18:27:17 16.095770 -79.676594
1489 2020-02-14 18:27:17 2020-02-14 18:27:17 15.638168 -79.699792
1490 2020-02-14 18:27:17 2020-02-14 18:27:17 15.173973 -78.676286
1491 2020-02-14 18:27:17 2020-02-14 18:27:17 15.022090 -78.618244
goes_imager_projection
0 -78.153419
1 -77.147658
2 -76.851406
3 -75.565914
4 -79.959360
... ...
1487 -77.354028
1488 -79.676594
1489 -79.699792
1490 -78.676286
1491 -78.618244
[1492 rows x 12 columns]
['x', 'y', 't', 'goes_imager_projection']
feature detection performed and saved
Segmentation:
Segmentation is performed based on the OLR field and a threshold value to determine the cloud areas.
[28]:
# Keyword arguments for the segmentation step:
parameters_segmentation={}
parameters_segmentation['target']='maximum'
parameters_segmentation['method']='watershed'
parameters_segmentation['threshold']=0.2
[29]:
# Perform segmentation and save results to files:
Mask_VIS,Features_VIS=tobac.themes.tobac_v1.segmentation(Features,new_data.CMI,dxy,**parameters_segmentation)
print('segmentation VIS performed, start saving results to files')
Features_VIS.to_netcdf(os.path.join(savedir,'Features_VIS.nc'))
Mask_VIS.to_netcdf(os.path.join(savedir,'Mask_Segmentation_VIS.nc'))
print('segmentation VIS performed and saved')
<xarray.DataArray 'CMI' (t: 6, y: 542, x: 542)>
array([[[0.17650777, 0.15174589, 0.09904753, ..., 0.07682532,
0.0873015 , 0.08920626],
[0.3238092 , 0.05301582, 0.11460306, ..., 0.0809523 ,
0.08603166, 0.0873015 ],
[0.1206348 , 0.047619 , 0.05142852, ..., 0.09587292,
0.0904761 , 0.08222214],
...,
[0.0396825 , 0.03746028, 0.0380952 , ..., 0.04507932,
0.04253964, 0.04222218],
[0.03714282, 0.03682536, 0.03841266, ..., 0.04412694,
0.04539678, 0.04571424],
[0.03904758, 0.03873012, 0.03936504, ..., 0.04063488,
0.0412698 , 0.04031742]],
[[0.11492053, 0.16380937, 0.06952374, ..., 0.07396818,
0.08825389, 0.09301578],
[0.14698398, 0.10730148, 0.06285708, ..., 0.07587294,
0.0825396 , 0.08444437],
[0.07015866, 0.0714285 , 0.05047614, ..., 0.09079356,
0.08380944, 0.08285706],
...
[0.03746028, 0.0380952 , 0.03746028, ..., 0.0587301 ,
0.03873012, 0.04507932],
[0.03619044, 0.03841266, 0.03873012, ..., 0.04158726,
0.04095234, 0.04063488],
[0.03682536, 0.0365079 , 0.03714282, ..., 0.1158729 ,
0.03777774, 0.0428571 ]],
[[0.04507932, 0.04349202, 0.05650788, ..., 0.07015866,
0.07999992, 0.0825396 ],
[0.04158726, 0.04380948, 0.04507932, ..., 0.07396818,
0.07873008, 0.08190469],
[0.0412698 , 0.0428571 , 0.0412698 , ..., 0.08634912,
0.08317453, 0.07587294],
...,
[0.04222218, 0.03841266, 0.03746028, ..., 0.0984126 ,
0.0730158 , 0.0412698 ],
[0.0380952 , 0.04634916, 0.0682539 , ..., 0.03873012,
0.03999996, 0.0380952 ],
[0.0365079 , 0.03714282, 0.03873012, ..., 0.03873012,
0.04031742, 0.03714282]]], dtype=float32)
Coordinates:
* x (x) float64 -79.99 -79.98 -79.98 ... -75.01 -75.0
* y (y) float64 19.99 19.99 19.98 ... 15.02 15.01 15.0
* t (t) datetime64[ns] 2020-02-14T18:02:17.194589056 ...
goes_imager_projection int64 0
Attributes:
long_name: ABI L2+ Cloud and Moisture Imagery reflectance fa...
standard_name: toa_bidirectional_reflectance
sensor_band_bit_depth: 12
valid_range: [ 0 4095]
units: 1
resolution: y: 0.000014 rad x: 0.000014 rad
grid_mapping: goes_imager_projection
cell_methods: t: point area: point
ancillary_variables: DQF
segmentation VIS performed, start saving results to files
segmentation VIS 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.
[30]:
# 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'
[31]:
# Perform linking and save results to file:
Track=tobac.themes.tobac_v1.linking_trackpy(Features,new_data.CMI,dt=dt,dxy=dxy,**parameters_linking)
Track.to_netcdf(os.path.join(savedir,'Track.nc'))
Frame 5: 257 trajectories present.
Visualisation of the tracks:
[32]:
axis_extent=[-79,-76,17.5,19] #xmin,xmax,ymin,ymax
[33]:
import cartopy.crs as ccrs
fig_map,ax_map=plt.subplots(figsize=(14,8),subplot_kw={'projection': ccrs.PlateCarree()})
ax_map=tobac.plot.map_tracks(Track,axis_extent=axis_extent,axes=ax_map)
[34]:
# Create animation of tracked clouds and outlines with VIS as a background field:
animation_test_tobac=tobac.plot.animation_mask_field(Track,Features,new_data.CMI,Mask_VIS,
axis_extent=axis_extent,figsize=(14,8),#orientation_colorbar='horizontal',pad_colorbar=0.2,
vmin=0,vmax=1,cmap='Blues_r',linewidth_contour=1,
plot_outline=True,plot_marker=True,marker_track='x',plot_number=True,plot_features=True)
[35]:
# Display animation:
from IPython.display import HTML, Image, display
HTML(animation_test_tobac.to_html5_video())
[35]:
[36]:
# 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,35,5),density=False,width_bar=2.5)
ax_lifetime.set_xlabel('lifetime (min)')
ax_lifetime.set_ylabel('counts')
[36]:
Text(0, 0.5, 'counts')
