It's easy to create rich visualizations of ice and snow thickness, growth, melt, and temperature profile using SIMB3 data. In this tutorial, we’ll show you how to do it using the Python programming language.

‍

Before we get started, you'll need: 

• An SIMB3 datasheet (or sample below)
• A code editor of choice (we recommend Visual Studio)

‍

If you'd like to skip the tutorial and get straight to plotting, you can download the script that we build in this tutorial here: 

• massBalancePlot.py
• SIMB3sampleDataSheet.csv

‍

Note: this tutorial assumes you have an up-to-date version of Python 3 installed on your computer and that you have basic familiarity with using a shell terminal. If you are unfamiliar with how to install or upgrade Python, here is an excellent tutorial that will walk you through how to get Python setup on MacOS, Windows, or Linux operating systems.‍‍‍

‍

A quick introduction

One of the best parts about working with SIMB3 data is that the results are highly visual. With the surface and bottom rangefinder data alone we can see how the ice (and snow!) change thickness through the season. By incorporating data from the vertical temperature string, we can examine the vertical conductive heat flux. If we combine all three measurements, we create a comprehensive picture of the ice state through time that we refer to as a plot a mass balance plot.

‍

The mass balance plot

A mass balance plot shows ice and snow growth, melt, and internal temperature on a single, easily understandable figure (Figure 1). On the x-axis is time and on the y-axis is thickness or depth. As the ice grows and melts, the distances recorded by the SIMB3 surface and bottom rangefinders lengthen and shorten. With a little knowledge of the SIMB3 dimensions and the deployment conditions, we can turn these distances into ice and snow thickness values. Over time, changes in these values represent ice and snow growth and melt.

‍

‍

In this tutorial, we'll be covering how to make the plot in Figure 1 in Python.

‍

Time-series line plots

SIMB3 is also equipped with sensors that measure upper-ocean temperature, battery voltage, and meteorological data such as air temperature and barometric pressure. Creating these plots is straightforward compared to the mass balance plot, so we're not going to cover it in this article.

‍

Getting started with Python

Okay, let's get right into. From the Real Time Data portal, navigate to the page for the SIMB3 which you’d like to analyze. Click the blue “Download CSV” button to download the most recent datasheet for the instrument. Alternatively, you can use your own file or this sample for this tutorial.

‍‍

Create a new Python file

Inside of your code editor or terminal, navigate to a directory where you'd like to store your files for this project. Create a new file and give it a name and .py extension. We created a directory named <bash inlineCode>SIMB3_Python<bash inlineCode> and a file named <bash inlineCode>massBalancePlot.py<bash inlineCode>. While not required, I'd recommend moving your freshly downloaded datasheet into this directory.

‍

In a MacOS/Linux terminal:

cd ~
mkdir SIMB3-Python
cd SIMB3-Python
touch massBalancePlot.py
 

‍

Note that the prefix <bash inlineCode>sudo<bash inlineCode> might be required if you are not logged in with administrative access.

Installing required libraries

We'll be using the Numpy, Pandas, and Matplotlib libraries to structure, manipulate, and eventually plot our data. To import them, open your Python file and add the following lines at the very top.

import numpy as np
import pandas as pd
from matplotlib import pyplot as plt   
 

‍

If you haven't used any or all of these libraries before, you will need to install them. This is done simply using <bash inlineCode>pip3 install<bash inlineCode>.

pip3 install numpy pandas matplotlib
 

‍

With the libraries installed, it's now possible to create the functions which will layer the filled-area plots required to generate the mass balance plot. We'll accomplish this by making a plotting function that takes in a datasheet as input.

‍

Creating the mass balance plot function

For the mass balance plot, copy the following code directly below your <python inlineCode>import<python inlineCode> statements. Don't worry, we'll explain the whole thing block-by-block below.

def plot_simb(dtc: np.ndarray, surface_distance, bottom_distance, timestamp, initial_snow_depth, distance_between_rangefinders):

    #calculate height of surface rangefinder above the ice surface
    height_above_ice = surface_distance[0] + initial_snow_depth

    #create arrays that store values for each timestamp
    snow_height = height_above_ice - surface_distance
    ice_thickness = distance_between_rangefinders - height_above_ice - bottom_distance
    ice_thickness = - ice_thickness #make negative for plotting because it's below waterline

    #define temperature string position (height above ice)
    temp_string_length = 3.84 #standard SIMB3 temperature string has 2 cm spacing and 192 thermistors
    temp_string_offset = 0 #the axial distance between the surface rangefinder and the first thermister on the DTC
    temp_string_top = height_above_ice - temp_string_offset #distance from the top of the ice to the first thermistor in the temperature string
    temp_string_bottom = temp_string_length - temp_string_top
    temp_string_bottom = temp_string_bottom * (-1) #make negative for plotting
    

    #shift excel serial date (1900 epoch) to the Python 1970 epoch by adding the number of days between 01/01/1900 and 01/01/1970 
    #then subtract 1 to account for the famous 1900 excel leap year timestamp bug
    timestamp = timestamp - 25568 - 1 
    
    #set plot bounds
    upper_bound = 1
    lower_bound = -3 
    
    #create mass balance plot
    fig, ax = plt.subplots()
    im = ax.imshow(dtc, cmap='jet', vmin=-30, vmax=0, extent=[timestamp[0], timestamp[-1], temp_string_bottom, temp_string_top],
                   aspect='auto', zorder=0)     #format ice temperature heatmap from dtc (temperature string) data
    plt.colorbar(im, ax=ax, label='Temperature (C)')

    #color in between the lines
    ax.fill_between(timestamp, snow_height, 0, color="grey", zorder=1)                  #plot snow height on top of ice
    ax.fill_between(timestamp, temp_string_top, snow_height, color='white', zorder=2)       #create whitespace above snow
    ax.fill_between(timestamp, ice_thickness, lower_bound, color='b', zorder=3)         #plot lower boundary of ice

    #bound and label y axis
    ax.set_ylim(lower_bound, upper_bound)
    plt.ylabel('Depth (m)')

    #format x-axis as date
    formatter = mdates.DateFormatter("%Y-%m-%d")
    ax.xaxis.set_major_formatter(formatter)

    #only label each month
    locator = mdates.MonthLocator()
    ax.xaxis.set_major_locator(locator)

    plt.show()
    

‍

Taking a look at the function declaration, we see that it takes values from the temperature string (<python inlineCode>dtc: np.ndarray<python inlineCode>), the surface rangefinder (<python inlineCode>surface_distance<python inlineCode>), the bottom rangefinder (<python inlineCode>bottom_distance<python inlineCode>), the time stamp (<python inlineCode>timestamp<python inlineCode>), the initial snow depth (<python inlineCode>initial_snow_depth)<python inlineCode>), and the rangefinders distance (<python inlineCode>distance_between_rangefiners<python inlineCode>) as inputs.

‍

Of these inputs, the <python inlineCode>surface_distance<python inlineCode>, <python inlineCode>bottom_distance<python inlineCode>, <python inlineCode>timestamp<python inlineCode> are 1D arrays of length 1 x N, where N is the length of datasheet (often the number of transmissions from the buoy). The temperature string array (<python inlineCode>dtc: np.ndarray<python inlineCode>) is 2D with a size of 192 x N. The <python inlineCode>initial_snow_depth<python inlineCode> and <python inlineCode>distance_between_sounders<python inlineCode> values are both scalars.

‍

Define reference height from ice surface

We tried to give the variables in this function intuitive names, but let's still go through it block-by-block to make sure there's no confusion. In the first row, we calculate the fixed distance between the surface rangefinder and the ice surface by summing the first value in the <python inlineCode>surface_distance<python inlineCode> array (i.e., the value immediately after the buoy was deployed) and the initial snow height. You can also declare this value manually using measurements from deployment.

#calculate height of surface rangefinder above the ice surface
height_above_ice = surface_distance[0] + initial_snow_depth
   

‍

Snow height and ice thickness

We then create two arrays that define the snow height (i.e., the snow thickness) and the ice thickness. The snow height is simply the height of the surface rangefinder above the ice (does not change once the SIMB3 is frozen in) minus the reading from the surface rangefinder.

#create arrays that store values for each timestamp
snow_height = height_above_ice - surface_distance
ice_thickness = distance_between_rangefinders - height_above_ice - bottom_distance
ice_thickness = - ice_thickness #make negative for plotting because it's below waterline

‍

We then compute the ice thickness by subtracting the height of the surface rangefinder above the ice (<python inlineCode>height_above_ice<python inlineCode>) and the bottom rangefinder value (<python inlineCode>bottom_distance<python inlineCode> from the total distance between rangefinders, <python inlineCode>distance_between_rangefinders<python inlineCode>. For a standard SIMB3, <python inlineCode>distance_between_rangefinders = 4.051<python inlineCode> meters.

‍

It should be noted that the "ice thickness" variable name is slightly inaccurate; <python inlineCode>ice_thickness<python inlineCode> as computed here is strictly the thickness of ice below freeboard. In the summer, surface melt will contribute to ice thickness reduction and is not captured by this variable.

‍

Temperature string height reference

Moving along to the second block, we create several variables which define the position of the temperature string relative to the SIMB3 surface rangefinder and to the ice surface.

#define temperature string position (height above ice)
temp_string_length = 3.84 #standard SIMB3 temperature string has 2 cm spacing and 192 thermistors
temp_string_offset = 0 #the axial distance between the surface rangefinder and the first thermister on the DTC
temp_string_top = height_above_ice - temp_string_offset #distance from the top of the ice to the first thermistor in the temperature string
temp_string_bottom = temp_string_length - temp_string_top
temp_string_bottom = temp_string_bottom * (-1) #make negative for plotting

‍

The standard SIMB3 Bruncin digital temperature string is 3.84 meters long and measures temperature at 2 cm intervals along the length (192 total values). During buoy manufacturing, the position of the first thermistor is placed 0.2 cm below the surface rangefinder. After installation and freeze-in, this becomes a fixed position above the ice (<python inlineCode>temp_string_top<python inlineCode>).

‍

Quick note: the SIMB3sampleDataSheet.csv used in this tutorial was generated by a past-generation SIMB3 which had a <inlineCode class=matlab>temp_string_offset = 0<inlineCode class=matlab>. All modern SIMB3s have a <inlineCode class=matlab>temp_string_offset = 0.2<inlineCode class=matlab>

‍

The temperature string bottom (i.e., the distance between the ice surface and the final thermistor on the chain) is simply the length of the temperature string minus the height of the first thermistor above the ice surface.

‍

Realign Excel timestamp

While this is not strictly necessary, to have dates to plot along the x-axis, we need to realign the Excel serial date to the 1970 epoch used by Python. To do this, just add the number of days (including leap years!) between 01/01/1900 and 01/01/1970. Note that you must also subtract one extra day to account for the infamous Excel leap year bug. For more on how dates are handled in excel, check this out.

#shift excel serial date (1900 epoch) to the Python 1970 epoch by adding the number of days between 01/01/1900 and 01/01/1970 
#then subtract 1 to account for the famous 1900 excel leap year timestamp bug
timestamp = timestamp - 25568 - 1 

‍

Plot window parameters

Finally, we set the bounds for our plot window. For most buoys, 1 m above the ice surface and 3 meters below the surface is good.

#set plot bounds
upper_bound = 1
lower_bound = -3

‍

Plot temperature string values

We're now ready to assemble the mass balance figure. We'll start by defining a <python inlineCode>plt.subplots()<python inlineCode> object, and then we'll pass it our temperature string data.

#create mass balance plot
fig, ax = plt.subplots()
im = ax.imshow(dtc, cmap='jet', vmin=-30, vmax=0, extent=[timestamp[0], timestamp[-1], temp_string_bottom, temp_string_top],
                   aspect='auto', zorder=0)     #format ice temperature heatmap from dtc (temperature string) data
plt.colorbar(im, ax=ax)
    

‍

Here, the <python inlineCode>imshow()<python inlineCode> function takes in our digital temperature string array (<python inlineCode>dtc<python inlineCode>) and the bounds of the mesh over which we want to plot (in this case, on the x-axis from the first to last entry in <python inlineCode>timestamp<python inlineCode> and on the y-axis from <python inlineCode>temp_string_bottom<python inlineCode> to <python inlineCode>temp_string_top<python inlineCode>).

‍

<python inlineCode>imshow()<python inlineCode> plots the temperature string data and <python inlineCode>plt.colorbar(im, ax=ax)<python inlineCode> adds a colorbar with bounds defined by <python inlineCode>vmin<python inlineCode> and <python inlineCode>vmax<python inlineCode>. Both <python inlineCode>vmin<python inlineCode> and <python inlineCode>vmax<python inlineCode> should be set so that they encompass any temperature value present in the temperature string array.

‍

Plot rangefinder values

As the final step in the creation of our mass balance plot function, we'll fill in the regions of our plot corresponding to air, snow, and water using the rangefinder values. This is straightforward using the <python inlineCode>fill_between()<python inlineCode> function.

#color in between the lines
    ax.fill_between(timestamp, snow_height, 0, color="grey", zorder=1)                  #plot snow height on top of ice
    ax.fill_between(timestamp, temp_string_top, snow_height, color='white', zorder=2)       #preate whitespace above snow
    ax.fill_between(timestamp, ice_thickness, lower_bound, color='b', zorder=3)         #plot lower boundary of ice
    

‍

We can then set the y-limits of the plot, format the x-axis as a date, and command the plot to show when the <python inlineCode>plot_simb<python inlineCode> function is called.

#bound y axis
ax.set_ylim(lower_bound, upper_bound)       

#format x-axis as date
formatter = mdates.DateFormatter("%Y-%m-%d")
ax.xaxis.set_major_formatter(formatter)

#only label each month
locator = mdates.MonthLocator()
ax.xaxis.set_major_locator(locator)
    

‍

Calling the mass balance plot function

To call our mass balance plot, we need to parse our CSV datasheet and pass the required parameters into the <python inlineCode>plot_simb()<python inlineCode> function. To do this, let's make a another function that takes our data path, initial snow depth, and plot range in as arguments and then calls <python inlineCode>plot_simb()<python inlineCode>.

def plot_from_csv(path_to_data: str, initial_snow_depth: float, plot_range):

    #try to read csv at input location
    try:
        simb_df = pd.read_csv(path_to_data,skiprows=plot_range[0], nrows=plot_range[1])      #create pandas dataframe from csv
        timestamp = simb_df.time_stamp.to_numpy()       #extracts time stamp from simb_df dataframe, converts to numpy array
    except FileNotFoundError:
        print('Invalid Path, please check file location.')
        return

    # Distance between sounders in meters, used for bounding the Ice Mass Balance Plot
    sounder_dist = 4.051

    # Extract relevant data from dataframe, convert to numpy array
    snow_dist = simb_df.surface_distance.to_numpy()                     # distance from upper sounder to snow
    water_depth = simb_df.bottom_distance.to_numpy()                    # distance from lower sounder to bottom of ice
    dtc_values = simb_df.filter(like='dtc_values_').to_numpy().T        # temperature string data at 2 cm intervals

    # Call plot function for ice mass balance plot
    plot_simb(dtc_values, snow_dist, water_depth, timestamp, initial_snow_depth, sounder_dist)
    

This function just extracts data from our datasheet as Pandas dataframes and passes it to our <python inlineCode>plot_simb()<python inlineCode> function. When using it, the only thing to watch out for is the <python inlineCode>plot_range<python inlineCode> parameter which cannot be longer than our datasheet (i.e., it cannot have more rows than our datasheet).

‍

As the final line in our script, we call our our <python inlineCode>plot_from_csv()<python inlineCode> function. For the SIMB3sampleDataSheet.csv used in this tutorial, we'll use 0.2 m as the <python inlineCode>initial_snow_depth<python inlineCode> and [0 3378] as the <python inlineCode>plot_range<python inlineCode>. Note: you'll need to use the correct path to your datasheet if it's not in your working directory.

plot_from_csv(r'SIMB3sampleDataSheet.csv', .2, [0, 3378])
    

‍

Running the script

To run the script and create the plot, go back to your terminal and run:

python3 massBalancePlot.py
    

You should see the plot in Figure 2 appear.

‍

‍

Just like that, we've created a mass balance plot from SIMB3 data using Python! Because we build this script using functions, it's easy to call for multiple SIMB3 datasets or for multiple iterations of the same dataset.

‍

I hope you found this tutorial helpful, and if you have any questions please reach out and we'll get back to you as soon as we can!

‍

Lastly, if Python isn't your jam, check out our other tutorial on how to Visualize SIMB3 data using Matlab.