Scipy - Feature Tracking Cookbook

Overview¶

Within this notebook, we will cover:

How to filter continuous SLP data to identify features of interest.
How to identify the lowest SLP minima associated with ETCs.
How to use Scipy’s label function to identify coherent objects within the SLP field.
How to track identified ETCs through successive time steps.

Prerequisites¶

Concepts	Importance	Notes
Intro to Cartopy	Helpful
Understanding of NetCDF	Necessary	Familiarity with metadata structure
Intro to Numpy	Necessary	Basic n-dimensional array indexing

Imports¶

import xarray as xr
import numpy as np
from scipy.ndimage import minimum_filter
from scipy.ndimage import label, generate_binary_structure
import matplotlib.pyplot as plt
import matplotlib as mpl
from matplotlib.animation import FuncAnimation
from IPython.display import HTML
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.colors as mcolors
from matplotlib.cm import ScalarMappable
import math

Core functions for this notebook¶

Distance function¶

This function serves to calculate the distance between two given latitudes, and longitudes in degrees.

def distance(lon1, lat1, lon2, lat2):
        radius = 6371
        dlat = math.radians(lat2 - lat1)
        dlon = math.radians(lon2 - lon1)
        a = math.sin(dlat/2) * math.sin(dlat/2) + math.cos(math.radians(lat1)) * math.cos(math.radians(lat2)) * math.sin(dlon/2) * math.sin(dlon/2)
        c = 2 * math.atan2(math.sqrt(a), math.sqrt(1-a))
        d = radius * c
        return d

Crosses dateline function¶

This function ensures that features adjecent to each other, but on the international dateline, are considered as the same feature.

def crosses_dateline(labeled_arr, resolution):
    #Shift 180 degrees
    shifted_arr = np.roll(labeled_arr, 180*resolution, axis = 1)

    #Examine a 3x3 region centered at the dateline
    for lat in range(0, len(shifted_arr) - 1):
    	sample_arr = shifted_arr[lat:lat+2, 179:181]

    	top_left = sample_arr[0, 0]
    	top_right = sample_arr[0, 1]
    	bottom_right = sample_arr[1, 1]
    	bottom_left = sample_arr[1, 0]
		
    	#Compare top left with position next to it
    	if(top_left != 0 and top_right != 0):
    		if(top_left != top_right):
    			num_to_change = top_left
    			shifted_arr = np.where(shifted_arr != num_to_change, shifted_arr, top_right)
    			sample_arr = shifted_arr[lat:lat+2, 179:181]

    	#Compare top left with position diagonal to it
    	if(top_left != 0 and bottom_right != 0):
    		if(top_left != bottom_right):
    			num_to_change = top_left
    			shifted_arr = np.where(shifted_arr != num_to_change, shifted_arr, bottom_right)
    			sample_arr = shifted_arr[lat:lat+2, 179:181]
		
    	#Compate bottom left with position next to it
    	if(bottom_left != 0 and bottom_right != 0):
    		if(bottom_left != bottom_right):
    			num_to_change = bottom_left
    			shifted_arr = np.where(shifted_arr != num_to_change, shifted_arr, bottom_right)
    			sample_arr = shifted_arr[lat:lat+2, 179:181]
    			#print('Changed arr bot left, bot right', sample_arr)
    	
    	#Compare bottom left with position diagonal to it
    	if(bottom_left != 0 and top_right != 0):
    		if(bottom_left != top_right):
    			num_to_change = bottom_left
    			shifted_arr = np.where(shifted_arr != num_to_change, shifted_arr, top_right)
    			sample_arr = shifted_arr[lat:lat+2, 179:181]
    		
    return np.roll(shifted_arr, 180, axis = 1)

Read in and visualize our raw data¶

ds = xr.open_dataset('../notebooks/data/SLP_ex.nc')
SLP = ds['SLP']
SLP_hPa = SLP / 100
lon = SLP_hPa['longitude'].values
lat = SLP_hPa['latitude'].values
lon2d, lat2d = np.meshgrid(lon, lat)
levels = np.arange(950, 1050, 10)

Function to plot an animation over each time step in our data¶

def update(t):
    """
    Animation function, called once for each timestep in the dataset.
    """
    # Reset the map on each timestep
    ax.clear()
    ax.set_extent(extent, crs=ccrs.PlateCarree())
    ax.set_aspect('auto')

    field = SLP_hPa.isel(time=t).values # total precipitation data for current timestep

    # Colored contour lines
    cs = ax.contour(
        lon2d, lat2d, field,
        levels=levels,
        cmap=cmap,         
        linewidths=0.9,
        transform=ccrs.PlateCarree()
    )

    # Create outlines for continents/coastlines and gridlines
    ax.coastlines(color="0.55", linewidth=0.6)
    gl = ax.gridlines(
        draw_labels=True,
        linewidth=0.4,
        color="0.5",
        alpha=0.25,
        linestyle="--"
    )
    
    gl.top_labels = False
    gl.right_labels = False
    ax.set_title(f"Sea Level Pressure | time index = {t}")

# Create map
vals = SLP_hPa.values
vmin = np.nanpercentile(vals, 1)
vmax = np.nanpercentile(vals, 99)

norm = mcolors.Normalize(vmin=vmin, vmax=vmax) # normalize
cmap = 'turbo'

fig, ax = plt.subplots(figsize=(9, 4), subplot_kw={"projection": ccrs.PlateCarree()})
extent = [lon.min(), lon.max(), lat.min(), lat.max()] # geographic boundaries of the map
ax.set_extent(extent, crs=ccrs.PlateCarree())

# Create colorbar
sm = ScalarMappable(norm=norm, cmap=cmap)
sm.set_array([])
cbar = fig.colorbar(
    sm,
    ax=ax,
    pad=0.15,
    orientation = 'horizontal'
)
cbar.set_label("SLP")

# Create animation
anim = FuncAnimation(
    fig,
    update,
    frames=SLP_hPa.sizes["time"],
    interval=300,
    blit=False,
)

plt.close(fig)
HTML(anim.to_jshtml())

Get our desired objects¶

Here is where we convert our xarray dataarray to a numpy array.

SLP_arr = SLP_hPa.to_numpy()

Choose a threshold for SLP¶

The threshold chosen here will serve as a dividing line between what is and is not considered a feature of interest. In the block of code below, we will choose a threshold of 1000 hPa. This means that all of our features of interest must be comprised of grid cells whose SLP values are 1000 hPa or less. Other thresholds could also be used, such as a departure from a long-term mean.

SLP_thresh = 1000 #Units of hPa

SLP_1000 = np.where(SLP_arr < SLP_thresh, SLP, 0)

Create our “binary structure”¶

The binary structure is a boolean array that sets the definition for our smallest possible feature of interest. It essentially looks at all grid cells in our numpy array and determines which neighboring grid cells can be connected to form a coherent spatial object. The first input into generate_binary_structure defines the number of dimensions over which an object can cover. Since our numpy array contains three dimensions (time, latitude, longitude), our desired number of dimensions will be two, corresponding to latitude and longitude. The second input defines the number of grid cells that must be connected to each other to be considered as an object. For our case here, we will define an object minimally as one that has at least two connected grid cells.

s = generate_binary_structure(2,2)

Run the label function¶

Here is where we actually find our desired objects given a 2-dimensional field. When we run the label function, there will be two outputs. The first output is a 2-dimensional array that is the exact same size as our input array. However, this new array contains integer values corresponding to the feature of interest. Typically, a value of zero means that the grid cell is just part of the background (i.e. not a feature of interest). Then all grid cells with a value of 1 correspond to a feature, all grid cells with a value of 2 correspond to another feature, and so on. The second output of label is the integer total of distinct objects that were found in the given 2-dimensional field.

In the block of code below, we first define our grid size resolution in degrees. Next, we set up empty arrays that will store the output from the label for each time step in our original data file. For labeled_arr, we want a 3-dimensional array to store the output (number of time steps, number of latitudes, number of longitudes). Since the second output of the label is only an integer, a 1-dimensional array will suffice.

Lastly, we loop through each time step in our data file. For our case, this is only 13 time steps. After running the label function in the for loop, we send the output to the ‘crosses_dateline’ function, which ensures that objects that cross the international dateline are counted as the same rather than separate objects. We then append our correct outputs to the arrays that we defined earlier.

resolution = 1 #In degrees

labeled_arr = np.zeros_like(SLP_1000)

num_features_arr = np.zeros(len(SLP_1000))

for time_step in np.arange(0, len(SLP_1000)):
    labeled_SLP_field, num_features = label(SLP_1000[time_step], s)
  
    correct_labeled_SLP_field = crosses_dateline(labeled_SLP_field, resolution)
    
    num_features_arr[time_step] = num_features

    labeled_arr[time_step] = correct_labeled_SLP_field

Create a map of our objects¶

Here, we can visualize the objects identified at each time step. The shading on the plot below corresponds to the numerical object identified by the label function.

def update(t):
    """
    Animation function, called once for each timestep in the dataset.
    """
    # Reset the map on each timestep
    ax.clear()
    ax.set_extent(extent, crs=ccrs.PlateCarree())
    ax.set_aspect('auto')

    field = SLP_hPa.isel(time=t).values # total precipitation data for current timestep

    # Colored contour lines
    cs = ax.pcolormesh(
        lon2d, lat2d, labeled_arr[t], shading = 'auto',
        cmap=cmap,         
        linewidths=0.9,
        transform=ccrs.PlateCarree(),
    )

    # Create outlines for continents/coastlines and gridlines
    ax.coastlines(color="0.55", linewidth=0.6)
    gl = ax.gridlines(
        draw_labels=True,
        linewidth=0.4,
        color="0.5",
        alpha=0.25,
        linestyle="--"
    )
    
    gl.top_labels = False
    gl.right_labels = False
    ax.set_title(f"Identified Low-Pressure Objects | time index = {t}")

# Create map
vals = labeled_arr
vmin = np.nanpercentile(vals, 1)
vmax = np.nanpercentile(vals, 99)

norm = mcolors.Normalize(vmin=vmin, vmax=vmax) # normalize
cmap = 'turbo'

fig, ax = plt.subplots(figsize=(9, 4), subplot_kw={"projection": ccrs.PlateCarree()})
extent = [lon.min(), lon.max(), lat.min(), lat.max()] # geographic boundaries of the map
ax.set_extent(extent, crs=ccrs.PlateCarree())

# Create colorbar
sm = ScalarMappable(norm=norm, cmap=cmap)
sm.set_array([])
cbar = fig.colorbar(
    sm,
    ax=ax,
    pad=0.15,
    orientation = 'horizontal'
)
cbar.set_label("Object Number")

# Create animation
anim = FuncAnimation(
    fig,
    update,
    frames=SLP_hPa.sizes["time"],
    interval=300,
    blit=False,
)

plt.close(fig)
HTML(anim.to_jshtml())

Find SLP Minima¶

Here, we utlize scipy’s minimum_filter function to identify the locations of local SLP minima. The minimum_filter function first takes in the n-dimensional array over which we wish to identify local minima. Then, minimum_filter requires a numerical value that defines the size of the grid over which minima will be found. In the following example, we locate minima within 5x5 grid boxes throughout our domain. The modes keyword tells us how to find minima at the edges of the domain. In the y-direction, we simply use the nearest grid cell within our domain, whereas in the x-direction we want to wrap to the other side of our domain since our field is continuous in the x-direction. Lastly, we specify over which axes we wish to calculate minima (latitude, longitude).

minima = minimum_filter(SLP_arr, 5, mode=['nearest', 'wrap'], axes = [1, 2])
minima_3 = (SLP_arr == minima)
min_times, min_lats, min_lons = np.where(1 == minima_3)

Plot SLP and identified minima¶

In the animation produced by the following block of code, you will see a plot of SLP, with all identified local minima shown as blue dots.

def update(t):
    """
    Animation function, called once for each timestep in the dataset.
    """

    where_time = np.where(min_times == t)

    proper_time_lats = min_lats[where_time]
    proper_time_lons = min_lons[where_time]
    # Reset the map on each timestep
    ax.clear()
    ax.set_extent(extent, crs=ccrs.PlateCarree())
    ax.set_aspect('auto')

    field = SLP_hPa.isel(time=t).values # total precipitation data for current timestep

    # Colored contour lines
    cs = ax.contourf(
        lon2d, lat2d, SLP_arr[t],
        cmap=cmap,         
        transform=ccrs.PlateCarree(),
    )

    cs = ax.scatter(
        lon[proper_time_lons], lat[proper_time_lats], 
        zorder = 2,
    )

    # Create outlines for continents/coastlines and gridlines
    ax.coastlines(color="0.55", linewidth=0.6)
    gl = ax.gridlines(
        draw_labels=True,
        linewidth=0.4,
        color="0.5",
        alpha=0.25,
        linestyle="--"
    )
    
    gl.top_labels = False
    gl.right_labels = False
    ax.set_title(f"SLP and its Minima | time index = {t}")

# Create map
vals = SLP_arr
vmin = np.nanpercentile(vals, 1)
vmax = np.nanpercentile(vals, 99)

norm = mcolors.Normalize(vmin=vmin, vmax=vmax) # normalize
cmap = 'turbo'

fig, ax = plt.subplots(figsize=(9, 4), subplot_kw={"projection": ccrs.PlateCarree()})
extent = [lon.min(), lon.max(), lat.min(), lat.max()] # geographic boundaries of the map
ax.set_extent(extent, crs=ccrs.PlateCarree())

# Create colorbar
sm = ScalarMappable(norm=norm, cmap=cmap)
sm.set_array([])
cbar = fig.colorbar(
    sm,
    ax=ax,
    pad=0.15,
    orientation = 'horizontal'
)
cbar.set_label("Object Number")

# Create animation
anim = FuncAnimation(
    fig,
    update,
    frames=SLP_hPa.sizes["time"],
    interval=300,
    blit=False,
)

plt.close(fig)
HTML(anim.to_jshtml())

Find main SLP minimum for each feature¶

Below, we loop through each time step and then loop through each identified feature at that time step. First, we identify the latitudes and longitudes of our features of interest at each time step. We then look at each identified feature and find the lowest SLP minimum that occurs within that feature.

#Find the main SLP center for each object at each time step
deepest_mins_list = []

for time_step in np.arange(0, len(SLP_arr)):
    where_time = np.where(min_times == time_step)
    proper_time_lats = min_lats[where_time]
    proper_time_lons = min_lons[where_time]
    
    deepest_mins = {}
    
    for label_num in range(1, int(labeled_arr[time_step].max() + 1)):

        mask = (labeled_arr[time_step, proper_time_lats, proper_time_lons] == label_num)

        if(not np.any(mask)):
            continue

        obj_lat = proper_time_lats[mask]
        obj_lon = proper_time_lons[mask]

        obj_slp = SLP_arr[time_step][obj_lat, obj_lon]

        lowest = np.argmin(obj_slp)

        deepest_mins[label_num] = (lat[obj_lat[lowest]], obj_lon[lowest])
    deepest_mins_list.append(deepest_mins.copy())
    deepest_mins.clear()

Track features through time¶

We now have all proper minima at each time step. It is time to stitch together the minima through successive time steps to form tracks. The goal is to take all identified minima at the first time step in our dataset and track them through time. We introduce a variable in the block below called still exists. As the name suggests, when this variable is True, the storm track still exists at a certain time step. If it is False, then the track no longer exists at a certain time step.

To track features through time, we take a look at the next time step (t+1) and caculate the distances between our minimum at time step t and all identified features at time step t+1. We then search for a minimum at t+1 that is within 500 km of our minimum at t. This minimum at t+1 is taken as the position of our feature of interest at the next time step.

#We now have all minima at each time step. Now, it's time to stitch the tracks together
#Start with time step 1, check to see what's within 500 km. If none, then the track ends

first_time_step = deepest_mins_list[0]

track_no = 1
for min_point in first_time_step:
    print('----------------------------')
    print('Track Number ' + str(track_no))
    time_ind = 0
    
    #Start with minimum for first object
    min_point_lat = first_time_step[min_point][0]
    min_point_lon = first_time_step[min_point][1]
    print('Time step 1:', min_point_lat, min_point_lon)
    still_exists = True

    while(still_exists):
        time_ind += 1

        if(time_ind > len(deepest_mins_list) - 1):
            break
        
        #Go to next time step
        next_time_step = deepest_mins_list[time_ind]
        
        distances = np.zeros((len(next_time_step.items())))
        
        ind = 0
        for feat in next_time_step.values():
            dist = distance(min_point_lon, min_point_lat, feat[1], feat[0])
            distances[ind] = dist

            ind += 1
        
        min_dist = np.argmin(distances)
        #print(min_dist, next_time_step)
        if(min_dist < 500):
            new_min_point_lat = list(next_time_step.items())[min_dist][1][0]
            new_min_point_lon = list(next_time_step.items())[min_dist][1][1]
            print('Time step ' + str(time_ind + 1) + ':', new_min_point_lat, new_min_point_lon)
        else:
            still_exists = False
    track_no += 1

----------------------------
Track Number 1
Time step 1: 78.0 97
Time step 2: 78.0 97
Time step 3: 78.0 98
Time step 4: 79.0 101
Time step 5: 79.0 102
Time step 6: 79.0 103
Time step 7: 79.0 104
Time step 8: 79.0 104
Time step 9: 79.0 104
Time step 10: 79.0 105
Time step 11: 79.0 105
Time step 12: 80.0 106
Time step 13: 80.0 107
----------------------------
Track Number 2
Time step 1: 60.0 270
Time step 2: 60.0 271
Time step 3: 60.0 271
Time step 4: 60.0 272
Time step 5: 60.0 273
Time step 6: 60.0 274
Time step 7: 60.0 274
Time step 8: 61.0 274
Time step 9: 61.0 275
Time step 10: 61.0 276
Time step 11: 61.0 276
Time step 12: 61.0 276
Time step 13: 61.0 277
----------------------------
Track Number 3
Time step 1: 60.0 179
Time step 2: 60.0 179
Time step 3: 60.0 180
Time step 4: 60.0 180
Time step 5: 60.0 180
Time step 6: 52.0 187
Time step 7: 52.0 187
Time step 8: 60.0 181
Time step 9: 53.0 189
Time step 10: 53.0 189
Time step 11: 53.0 190
Time step 12: 54.0 190
Time step 13: 54.0 191
----------------------------
Track Number 4
Time step 1: 59.0 56
Time step 2: 55.0 58
Time step 3: 59.0 56
Time step 4: 56.0 59
Time step 5: 51.0 54
Time step 6: 51.0 54
Time step 7: 51.0 54
Time step 8: 51.0 54
Time step 9: 51.0 54
Time step 10: 51.0 54
Time step 11: 51.0 54
Time step 12: 52.0 54
Time step 13: 52.0 54
----------------------------
Track Number 5
Time step 1: 59.0 216
Time step 2: 59.0 216
Time step 3: 59.0 216
Time step 4: 60.0 180
Time step 5: 60.0 180
Time step 6: 52.0 187
Time step 7: 52.0 187
Time step 8: 60.0 181
Time step 9: 53.0 189
Time step 10: 53.0 189
Time step 11: 53.0 190
Time step 12: 54.0 190
Time step 13: 54.0 191
----------------------------
Track Number 6
Time step 1: 57.0 339
Time step 2: 57.0 339
Time step 3: 58.0 339
Time step 4: 58.0 339
Time step 5: 58.0 339
Time step 6: 58.0 339
Time step 7: 58.0 339
Time step 8: 59.0 339
Time step 9: 59.0 339
Time step 10: 59.0 339
Time step 11: 59.0 339
Time step 12: 59.0 339
Time step 13: 59.0 340
----------------------------
Track Number 7
Time step 1: 50.0 53
Time step 2: 50.0 53
Time step 3: 50.0 53
Time step 4: 50.0 54
Time step 5: 51.0 54
Time step 6: 51.0 54
Time step 7: 51.0 54
Time step 8: 51.0 54
Time step 9: 51.0 54
Time step 10: 51.0 54
Time step 11: 51.0 54
Time step 12: 52.0 54
Time step 13: 52.0 54
----------------------------
Track Number 8
Time step 1: 42.0 167
Time step 2: 42.0 168
Time step 3: 42.0 168
Time step 4: 42.0 168
Time step 5: 42.0 168
Time step 6: 42.0 168
Time step 7: 42.0 169
Time step 8: 41.0 169
Time step 9: 41.0 169
Time step 10: 41.0 169
Time step 11: 41.0 169
Time step 12: 41.0 169
Time step 13: 41.0 169
----------------------------
Track Number 9
Time step 1: 18.0 116
Time step 2: 19.0 116
Time step 3: 19.0 116
Time step 4: 19.0 116
Time step 5: 19.0 116
Time step 6: 19.0 116
Time step 7: 19.0 116
Time step 8: 19.0 116
Time step 9: 19.0 116
Time step 10: 19.0 116
Time step 11: 19.0 116
Time step 12: 19.0 116
Time step 13: 20.0 116
----------------------------
Track Number 10
Time step 1: -63.0 69
Time step 2: -64.0 71
Time step 3: -64.0 72
Time step 4: -64.0 72
Time step 5: -64.0 72
Time step 6: -64.0 72
Time step 7: -64.0 73
Time step 8: -64.0 73
Time step 9: -64.0 74
Time step 10: -65.0 76
Time step 11: -65.0 76
Time step 12: -65.0 76
Time step 13: -65.0 77
----------------------------
Track Number 11
Time step 1: -42.0 85
Time step 2: -64.0 71
Time step 3: -64.0 72
Time step 4: -64.0 72
Time step 5: -64.0 72
Time step 6: -64.0 72
Time step 7: -64.0 73
Time step 8: -64.0 73
Time step 9: -64.0 74
Time step 10: -65.0 76
Time step 11: -65.0 76
Time step 12: -65.0 76
Time step 13: -65.0 77
----------------------------
Track Number 12
Time step 1: -81.0 77
Time step 2: -80.0 91
Time step 3: -64.0 72
Time step 4: -64.0 72
Time step 5: -64.0 72
Time step 6: -64.0 72
Time step 7: -64.0 73
Time step 8: -64.0 73
Time step 9: -64.0 74
Time step 10: -65.0 76
Time step 11: -65.0 76
Time step 12: -65.0 76
Time step 13: -65.0 77