Tutorial

Welcome to the STEP tutorial! Here we give some examples of how to use the package to help users get started. Please see the Implementation Details for more on function specifics and some general usage tips.

First, let's create a new project directory.

$ mkdir step_tutorial

And a file tutorial.py.

$ touch tutorial.py

Next, we'll download the tutorial_files directory, place it in step_tutorial, and open up tutorial.py in our favorite code editor.

To begin, let's import STEP and some additional dependencies we'll need.

from step import identification as idf
from step import quantification as qu
from step import tracking as tr
from step import visualization as viz

import copy
import matplotlib.pyplot as plt
import netCDF4 as nc
import numpy as np
from tutorial_files import plot_with_map as map_plot
from skimage import draw

Once we've done this, we can load in the precipitation data as shown below.

precip_data = np.load('tutorial_files/precip_1996.npy', allow_pickle=True)

We can also read-in data from an .nc file like so. In this case, we extract the corresponding latitude and longitude data.

file = 'tutorial_files/lat_long_1996.nc'
data = nc.Dataset(file)
lat_data = data['XLAT'][:]
long_data = data['XLONG'][:]

Before we move on, it's important to get a good idea of the geographical orientation of our data. If we take a look at the latitude and longitude data, we'll find that our precipitation is situated in and around the United States and Canada. Importantly though, both are given with the common orientation:

		South
	West		East
		North

This doesn't affect how we use STEP, but it's important to know when computing the central location of each storm and visualizing our results. But enough about that for now, let's start using the package!

If we need to find an appropriate precipitation threshold, we can create a histogram from the precip_data to aid in finding a precipitation threshold.

viz.histogram(precip_data, 3, (0, 2), 'show')

Let's use 0.6.

THRESHOLD = 0.6 

With our previously defined precipitation threshold, we can use Numpy to help us remove the data we don't wish to use in our calculations.

trimmed_data = np.where(precip_data < THRESHOLD, 0, precip_data)

If we want, we can see and save a plot of the precipitation intensities using the visualization function intensities. But first, we need to give decide on a Matplotlib colormap for our intensity plot.

Let's go with this one!

colormap = plt.get_cmap('gnuplot').reversed()

Now that we're all set to plot, let's specify a title, unit, and start time. Also, let's just show the image for now.

viz.intensities(precip_data, colormap, 'Precipitation Intensities 1996',
                           'mm (3h\N{SUPERSCRIPT MINUS}\N{SUPERSCRIPT ONE})', start_time=1, show_save='show')

Great! Now, to identify storms we first need to construct a custom structural set (most likely a disk). We can use scikit-image to do this. For more information on creating such a disk, check out the scikit-image docs.

struct = np.zeros((16, 16))
rr, cc = draw.disk((7.5, 7.5), radius=8.5)
struct[rr, cc] = 1
# https://stackoverflow.com/a/41495033

The algorithms take 3D (Time x Rows x Cols) arrays, so if we have a 2D array, we need to reshape it to include a phony time dimension; we'll be fine without doing so here though.

if trimmed_data.ndim == 2:
   trimmed_data = trimmed_data.reshape(1, trimmed_data.shape[0], trimmed_data.shape[1])

To compute the identification algorithm, we simply supply our data and our morphological structure.

labeled_maps = idf.identify(trimmed_data, struct)

To see the results of our work, we can use the visualization function storms, but first let's try a different colormap.

colormap = plt.get_cmap('hsv')

Let's provide a title, a start time, and we'll just show the images again.

viz.storms(labeled_maps, colormap, 'Identified Storms 1996', 1, show_save='show')

To save the result as a handy Numpy array, we can use Numpy save.

np.save('labeled_maps.npy', labeled_maps)

And load it, just like we did above to read-in our precipitation data.

labeled_maps = np.load('labeled_maps.npy', allow_pickle=True)

Now, let's track our identified storms, but first, we might want to keep a copy of labeled_maps if we haven't saved it before tracking storms, since it will be passed by reference.

id_result = copy.deepcopy(labeled_maps)

We can track storms by supplying the necessary data and user-specified values –– again, see Implementation Details for more on this.

If we're still figuring out the optimal parameter values for our data, we can turn test on, which prints information about the decisions the algorithm is making, which in turn can help us tune our parameters. Note: the result of this calculation is provided for the user below, since it is recommended that storm tracking should only be undertaken on machines designed specifically for tasks of such computational weight.

tracked_storms = tr.track(labeled_maps, trimmed_data, 0.7, 0.003, 18.6, test=True)

If we've decided to track the storms, we can save this if we're happy with the result –– this is probably a good idea.

np.save('tracked_storms_1996.npy', tracked_storms)

Otherwise, let's load the result of the calculation above to save time.

tracked_storms = np.load('tutorial_files/tracked_storms_1996.npy', allow_pickle=True)

Now, we can use storms again for a sample plot of our results.

viz.storms(tracked_storms, colormap, 'Tracked Storms 1996', 1, show_save='show')

If we wish to plot over a map, we can do so using a function like the example one used below.

But first, let's give our data a flip just for practice. It's important to remember that when transforming data, we need to make sure we perform the same transformation on any data we might use in the future. We'll do so on the following.

tracked_storms = np.flip(tracked_storms, 1)
lat_data = np.flip(lat_data, 1)
trimmed_data = np.flip(trimmed_data, 1)

Now we can plot over a map!

map_plot.storms_with_map(tracked_storms, 'Tracked Storms with Map 1996', lat_data, long_data, 1)

Lastly, we can also compute characteristics of these storms by supplying some additional information about the tutorial data to quantify. Good thing we flipped our trimmed precipitation data too.

storm_chars = qu.quantify(tracked_storms, trimmed_data, lat_data, long_data, 3, 16)

Now that we have the quantification results, to view the associated metrics for the storm labeled 2, we can index the tuple returned by quantify.

To view the duration of storm 2, we just have to specify the duration metric and the storm label.

print(f'Duration: {storm_chars[0][2]}')

To view its size, average intensity, or central location in time slice 3, we just have to specify the metric, time slice, and storm of interest.

print(f'Size: {storm_chars[1][3][2]}')
print(f'Intensity: {storm_chars[2][3][2]}')
print(f'Central location: {storm_chars[3][3][2]}')

Hope this helps! For function specifics, see the Implementation Details.