Winter Arctic sea ice state variability#

Summary: In this notebook, we highlight key additional sea ice variables: snow depth, ice type, snow density, and sea ice drift. We’ll use cartopy and xarray to generate maps and lineplots of the data to demonstrate methods for visualizing the data statically, as opposed to the interactive plotting functions highlighted in the seperate interactive notebooks (which can be slow to render and push the GitHub file size limits!).

Please take a look at the ‘Data variables’ tab in the book_ds cell below to explore the potential variables for analysis in this notebook.

The dataset also includes a few interpolated/smoothed variables (freeboard_int, snow_depth_int, ice_thickness_int) that can be used instead of the raw monthly means to increase coverage in some regions. See the interp_demo notebook for more information on how they are derived.

The analysis presented here was peer-reviewed in this paper in The Cryosphere.

Version history: Version 1 (01/01/2022)

Import notebook dependencies#

# For working with gridded climate data 
import xarray as xr 
# Helper function for reading the data from the bucket
from utils.read_data_utils import read_book_data 
from utils.plotting_utils import static_winter_comparison_lineplot, staticArcticMaps, staticArcticMaps_overlayDrifts, interactiveArcticMaps, compute_gridcell_winter_means # Plotting utils 
import numpy as np
# Plotting dependencies
#%config InlineBackend.figure_format = 'retina'
import matplotlib as mpl
# Sets figure size in the notebook
mpl.rcParams['figure.dpi'] = 150 

# Remove warnings to improve display
import warnings 
warnings.filterwarnings('ignore')

Read in the aggregated monthly gridded winter Arctic sea ice data#

book_ds = read_book_data() # Read/download the data 
book_ds

<xarray.Dataset>
Dimensions:               (time: 30, y: 448, x: 304)
Coordinates:
  * time                  (time) datetime64[ns] 2018-11-01 ... 2021-04-01
    longitude             (y, x) float32 ...
    latitude              (y, x) float32 ...
    xgrid                 (y, x) float32 ...
    ygrid                 (y, x) float32 ...
Dimensions without coordinates: y, x
Data variables: (12/20)
    projection            (time) float64 ...
    ice_thickness         (time, y, x) float32 ...
    ice_thickness_unc     (time, y, x) float32 ...
    num_segments          (time, y, x) float32 ...
    mean_day_of_month     (time, y, x) float32 ...
    snow_depth            (time, y, x) float32 ...
    ...                    ...
    region_mask           (time, y, x) float32 ...
    piomas_ice_thickness  (time, y, x) float32 ...
    t2m                   (time, y, x) float32 ...
    msdwlwrf              (time, y, x) float32 ...
    x_vel                 (time, y, x) float64 ...
    y_vel                 (time, y, x) float64 ...
Attributes:
    contact:      Alek Petty (alek.a.petty@nasa.gov)
    description:  Gridded Oct 2019 Arctic sea ice thickness and ancillary dat...
    reference:    Data doi: 10.5067/CV6JEXEE31HF. Original methodology descri...
    history:      Created 21/01/22

Static winter mean maps#

Compute and map (static) gridcell winter means for given variables

years = [2018,2019,2020] # Years over which to perform analysis

print(compute_gridcell_winter_means.__doc__) # Print docstring
Hide code cell output 
 Compute winter means over the time dimension. Useful for plotting as the grid is maintained. 
    
    Args: 
        da (xr.Dataset or xr.DataArray): data to restrict by time; must contain "time" as a coordinate 
        years (list of str): years over which to compute mean (default to unique years in the dataset)
        year_start (str, optional): year to start time range; if you want Nov 2019 - Apr 2020, set year="2019" (default to the first year in the dataset)
        start_month (str, optional): first month in winter (default to November)
        end_month (str, optional): second month in winter; this is the following calender year after start_month (default to April)
        force_complete_season (bool, optional): require that winter season returns data if and only if all months have data? i.e. if Sep and Oct have no data, return nothing even if Nov-Apr have data? (default to False) 
    
    Returns: 
        merged (xr.DataArray): DataArray with winter means as a time coordinate
    

print(staticArcticMaps.__doc__) # Print docstring
Hide code cell output 
 Show data on a basemap of the Arctic. Can be one month or multiple months of data. 
    Creates an xarray facet grid. For more info, see: http://xarray.pydata.org/en/stable/user-guide/plotting.html
    
    Args: 
        da (xr DataArray): data to plot
        title (str, optional): title string for plot
        dates (str list, option): dates to assign to subtitles, else defaults to whatever cartopy thinks they are
        out_str (str, optional): output string when saving
        cmap (str, optional): colormap to use (default to viridis)
        col (str, optional): coordinate to use for creating facet plot (default to "time")
        col_wrap (int, optional): number of columns of plots to display (default to 3, or None if time dimension has only one value)
        vmin (float, optional): minimum on colorbar (default to 1st percentile)
        vmax (float, optional): maximum on colorbar (default to 99th percentile)
        min_lat (float, optional): minimum latitude to set extent of plot (default to 50 deg lat)
        set_cbarlabel (str, optional): set colorbar label
        savefig (bool): output figure
    
    Returns:
        Figure displayed in notebook 
    
    

freeboard_winter_means = compute_gridcell_winter_means(book_ds.freeboard, years=years)
staticArcticMaps(freeboard_winter_means, dates=freeboard_winter_means.time.values, set_cbarlabel = "Sea ice freeboard (m)", cmap="YlOrRd", vmin=0, vmax=0.8, out_str='freeboard_winter')
../_images/d6939aff97467c6ed605b982f1fedf09fea9d4888e5d61821830cc0bcf6d1846.png 
snow_depth_winter_means = compute_gridcell_winter_means(book_ds.snow_depth, years=years)
staticArcticMaps(snow_depth_winter_means, dates=snow_depth_winter_means.time.values, set_cbarlabel = "Snow depth (m)", cmap="inferno", vmin=0, vmax=0.5, out_str='snowdepth_winter')
../_images/4050bb515dec8b0dcd450fe6203caa76fe8a3cade23785b477bdf4acdb65172a.png 
snow_density_winter_means = compute_gridcell_winter_means(book_ds.snow_density, years=years)
staticArcticMaps(snow_density_winter_means, dates=snow_depth_winter_means.time.values, set_cbarlabel = "Snow density (kg/m3)", cmap="GnBu", vmin=240, vmax=330, out_str='snowdensity_winter')
../_images/8a6ea1b547300566d4fed4adcde06b1fd5307b1e05473f77efad9ecc4f5d5a68.png 
thickness_winter_means = compute_gridcell_winter_means(book_ds.ice_thickness, years=years)
staticArcticMaps(thickness_winter_means, dates=thickness_winter_means.time.values, title="", set_cbarlabel = "Sea ice thickness (m)", cmap="viridis", vmin=0, vmax=5, out_str='thickness_winter')
../_images/fe72e8e4bc596aeed823028200bb327cb67282ec4cf61e7645ed652f8a621652.png 
ice_type_winter_means = compute_gridcell_winter_means(book_ds.ice_type, years=years)
staticArcticMaps(ice_type_winter_means, dates=ice_type_winter_means.time.values, set_cbarlabel = "Sea ice type (0 = FYI, 1 = MYI)", cmap="YlOrRd", vmin=0, vmax=1, out_str='icetype_winter')
../_images/83c93da3091e899e4df55503dc293e61f735f7267f662ecced23b7de036de8ee.png 
ice_conc_winter_means = compute_gridcell_winter_means(book_ds.sea_ice_conc, years=years)
staticArcticMaps(ice_conc_winter_means, dates=ice_conc_winter_means.time.values, set_cbarlabel = "Sea ice concentration", cmap="Blues_r", vmin=0, vmax=1, out_str='iceconc_winter')
../_images/e5c5bcdeb18286504fdc875d8982667dfd13d4fa4b817133fc02fc7a178bd7c6.png

Overlay sea ice drift vectors#

We can use a modified version of the plotting function used above to overlay sea ice drift vectors on any variable of interest. Below, we’ll generate the plot presented in Petty et al. (2022) that shows mean sea ice thickness and mean sea ice drift for the three winters.

print(staticArcticMaps_overlayDrifts.__doc__)
Hide code cell output 
 Show data on a basemap of the Arctic. Can be one month or multiple months of data. Overlay drift vectors on top 
    Creates an xarray facet grid. For more info, see: http://xarray.pydata.org/en/stable/user-guide/plotting.html
    
    Args: 
        da (xr DataArray): data to plot
        drifts_x (xr.DataArray): sea ice drifts along-x component of the ice motion
        drifts_y (xr.DataArray): sea ice drifts along-y component of the ice motion
        alpha (float 0-1, optional): Set this variable if you want da to have a reduced opacity (default to 1)
        res (int, optional): resolution of vectors (default to 6; plot 1 out of every 6 vectors)
        title (str, optional): title string for plot
        out_str (str, optional): output string when saving
        cmap (str, optional): colormap to use (default to viridis)
        col (str, optional): coordinate to use for creating facet plot (default to "time")
        col_wrap (int, optional): number of columns of plots to display (default to 3, or None if time dimension has only one value)
        vmin (float, optional): minimum on colorbar (default to 1st percentile)
        vmax (float, optional): maximum on colorbar (default to 99th percentile)
        min_lat (float, optional): minimum latitude to set extent of plot (default to 50 deg lat)
        set_cbarlabel (str, optional): set colorbar label
        savefig (bool): output figure
    
    Returns:
        Figure displayed in notebook 
    
    

# Compute winter means 

drifts_xvel_winter_means = compute_gridcell_winter_means(book_ds.x_vel, years=years)
drifts_yvel_winter_means = compute_gridcell_winter_means(book_ds.y_vel, years=years)
drifts_mag_winter_means = compute_gridcell_winter_means(np.sqrt(book_ds.x_vel**2+book_ds.y_vel**2).rename('mag'), years=years)

# Generate plot 
staticArcticMaps_overlayDrifts(da=drifts_mag_winter_means,
                                           dates=drifts_mag_winter_means.time.values,
                                           drifts_x=drifts_xvel_winter_means, 
                                           drifts_y=drifts_yvel_winter_means, 
                                           set_cbarlabel="Sea ice drift speed (m)", cmap="cubehelix_r", 
                                           vmin=0, vmax=0.5, alpha=0.8, res = 6,
                                           out_str='driftmag_drifts_overlayed')
../_images/85bcc18e70814787bc8226291cebb52fb321272946882085645869f12abb30e9.png 
# Compute winter means 
thickness_winter_means = compute_gridcell_winter_means(book_ds.ice_thickness, years=years)
drifts_xvel_winter_means = compute_gridcell_winter_means(book_ds.x_vel, years=years)
drifts_yvel_winter_means = compute_gridcell_winter_means(book_ds.y_vel, years=years)

# Generate plot 
pl_drifts = staticArcticMaps_overlayDrifts(da=thickness_winter_means,
                                           dates=thickness_winter_means.time.values,
                                           drifts_x=drifts_xvel_winter_means, 
                                           drifts_y=drifts_yvel_winter_means, 
                                           set_cbarlabel="Sea ice thickness (m)", cmap="viridis", 
                                           vmin=0, vmax=5, alpha=0.8, 
                                           out_str='thickness_drifts_overlayed')
display(pl_drifts)
../_images/41663f34e444bc6dac50a9bca32d45260533c0a28365e4e3e35da306197e55f0.png

You can also plot the sea ice drifts for just a single month. Below, we’ll show November sea ice drift overlayed on November sea ice freeboard for three winter seasons.

# Generate plot 
pl_drifts = staticArcticMaps_overlayDrifts(da=book_ds.freeboard.sel(time=["Nov 2018","Nov 2019","Nov 2020"]),
                                           drifts_x=book_ds.x_vel.sel(time=["Nov 2018","Nov 2019","Nov 2020"]), 
                                           drifts_y=book_ds.y_vel.sel(time=["Nov 2018","Nov 2019","Nov 2020"]), 
                                           set_cbarlabel="Sea ice freeboard (m)", cmap="YlOrRd", 
                                           vmin=0, vmax=0.8, alpha=0.8, 
                                           out_str='nov_freeboard_drifts_overlayed')
display(pl_drifts)
../_images/5e93693c96df5abc2df03187fb2e67e1f9831b64948fdf7b5b08ac2eda45a109.png

We can also change the resolution of the drifts to see more or less detail in the overlaid vectors. You can do this by changing the res argument, set to a default of 6. You also might want to change the opacity of the data to improve the visualization, which can be altered with the alpha argument. Setting alpha=1 indicates you want full opacity, and less than 1 indicates for some degree of transparency.

Below, we’ll plot sea ice type with an alpha of 0.7 and a drift vector resoltion of 11 (display 1 out of every 11 vectors).

# Generate plot 
pl_drifts = staticArcticMaps_overlayDrifts(da=book_ds.ice_type.sel(time=slice("Nov 2020", "April 2021")),
                                           drifts_x=book_ds.x_vel.sel(time=slice("Nov 2020", "April 2021")), 
                                           drifts_y=book_ds.y_vel.sel(time=slice("Nov 2020", "April 2021")), 
                                           set_cbarlabel="Sea ice type (FYI=0, MYI=1)", cmap="YlOrRd", 
                                           vmin=0, vmax=1, alpha=0.7, res = 11,  
                                           out_str='ice_type_drifts_overlayed')
display(pl_drifts)
../_images/cc2d847d3aaf9aaf31f161f328f05eb7f0367679abb27755a9a954c9e1361fce.png

Static monthly maps#

Compute and map (static) data given variables across the three winters for individual months

pl = staticArcticMaps(book_ds.freeboard.sel(time=["Nov 2018","Nov 2019","Nov 2020"]), set_cbarlabel = "Sea ice freeboard", cmap="YlOrRd", vmin=0, vmax=0.8, out_str='freeboard_november')
display(pl)
../_images/453863026779b88030c9750761bd317a76c4e82e78330f6a02158f8a70469942.png 
pl = staticArcticMaps(book_ds.freeboard.sel(time=["Apr 2019","Apr 2020","Apr 2021"]), set_cbarlabel = "Sea ice freeboard", cmap="YlOrRd", vmin=0, vmax=0.8, out_str='freeboard_april')
display(pl)
../_images/2bb405aca022ed586aa9d5e96f5335bc0bdf6cc57bd44b01ef8eae8f5702119b.png 
pl = staticArcticMaps(book_ds.snow_depth.sel(time=["Nov 2018","Nov 2019","Nov 2020"]), set_cbarlabel = "Snow depth", cmap="inferno", vmin=0, vmax=0.5, out_str='snowdepth_november')
display(pl)
../_images/94888d19a12396caa8489d955be82cd90d0452d8469a7e7ff80fe4404728f8ac.png 
pl = staticArcticMaps(book_ds.snow_depth.sel(time=["Apr 2019","Apr 2020","Apr 2021"]), set_cbarlabel = "Snow depth", cmap="inferno", vmin=0, vmax=0.5, out_str='snowdepth_april')
display(pl)
../_images/d6954bded48b00b48a8e2ef18245acdb89092d50f90201a4b6e79b541531cb7c.png 
pl = staticArcticMaps(book_ds.ice_thickness.sel(time=["Nov 2018","Nov 2019","Nov 2020"]), set_cbarlabel = "Sea ice thickness", cmap="viridis", vmin=0, vmax=5, out_str='thickness_november')
display(pl)
../_images/e62facbec29afb1ad95afe3ea768d487436df01bf391e58f1067fd3a1967f89f.png 
pl = staticArcticMaps(book_ds.ice_thickness.sel(time=["Apr 2019","Apr 2020","Apr 2021"]), set_cbarlabel = "Sea ice thickness", cmap="viridis", vmin=0, vmax=5, out_str='thickness_april')
display(pl)
../_images/8acdbd5beb573a66c246e6a7acfbbe374bd481031ca5613a267c6c7e1711215e.png

Thickness anomaly plots#

Compute and map (static) monthly anomaly maps (relative to the mean across the three winters by default) for given variables across the three winters

freeboard_winter_means = compute_gridcell_winter_means(book_ds.freeboard, years=years)
pl = staticArcticMaps(freeboard_winter_means-freeboard_winter_means.mean(axis=0), title="", set_cbarlabel = "Sea ice freeboard anomaly (m)", cmap="BrBG", vmin=-0.2, vmax=0.2, out_str='freeboard_winter_anomaly')
display(pl)
../_images/fd58937bbfa39a0c850792c62cd42ae558fb9d11370974f46bf208a80aa025ee.png 
snow_depth_winter_means = compute_gridcell_winter_means(book_ds.snow_depth, years=years)
pl = staticArcticMaps(snow_depth_winter_means-snow_depth_winter_means.mean(axis=0), set_cbarlabel = "Snow depth anomaly (m)", cmap="PRGn", vmin=-0.2, vmax=0.2, out_str='snowdepth_winter_anomaly')
display(pl)
../_images/4e391b44e65da5f42e3ea172c3f185b7f7a6b3d11b86864a32a4ea5dc1bedbcd.png 
thickness_winter_means = compute_gridcell_winter_means(book_ds.ice_thickness, years=years)
pl = staticArcticMaps(thickness_winter_means-thickness_winter_means.mean(axis=0), set_cbarlabel = "Sea ice thickness anomaly (m)", cmap="RdBu", vmin=-1.5, vmax=1.5, out_str='thickness_winter_anomaly')
display(pl)
../_images/252e88843a33805f0d2d137d0b9d56a86644c2b0dc5060b6a549a58102b44463.png

Do the same for the monthly anomalies (November then April)

pl = staticArcticMaps(book_ds.freeboard[0::12]-book_ds.freeboard[0::12].mean(axis=0), set_cbarlabel = "Sea ice freeboard anomaly (m)", cmap="BrBG", vmin=-0.2, vmax=0.2, out_str='freeboard_november_anomaly')
display(pl)
../_images/922d2a60fedfa80a957017a0f635aaaa66633901fe8c491d37dff6b0a1a9ec0a.png 
pl = staticArcticMaps(book_ds.snow_depth[0::12]-book_ds.snow_depth[0::12].mean(axis=0), set_cbarlabel = "Snow depth anomaly (m)", cmap="PRGn", vmin=-0.2, vmax=0.2, out_str='snow_depth_november_anomaly')
display(pl)
../_images/a95151498169fba3f14e25a2b33a61b32d150f38318311f352a5a1099a864978.png 
pl = staticArcticMaps(book_ds.ice_thickness[5::12]-book_ds.ice_thickness[5::12].mean(axis=0), set_cbarlabel = "Sea ice thickness anomaly (m)", cmap="RdBu", vmin=-1.5, vmax=1.5, out_str='thickness_november_anomaly')
display(pl)
../_images/ef92e782b9e6cff28786d0560d42a3f1cac520c6c64d5bb224c5ca473374b47d.png 
pl = staticArcticMaps(book_ds.freeboard[5::12]-book_ds.freeboard[5::12].mean(axis=0), set_cbarlabel = "Sea ice freeboard anomaly (m)", cmap="BrBG", vmin=-0.2, vmax=0.2, out_str='freeboard_april_anomaly')
display(pl)
../_images/a70ac3117ef512c97913db77855447f1010e2c12aeb53801b2a4881284a49113.png 
pl = staticArcticMaps(book_ds.snow_depth[5::12]-book_ds.snow_depth[5::12].mean(axis=0), set_cbarlabel = "Snow depth anomaly (m)", cmap="PRGn", vmin=-0.2, vmax=0.2, out_str='snow_depth_april_anomaly')
display(pl)
../_images/4848b5bb362b0e69dcf11a2bc70b21c2aec98a2780bdd28c20a71c13c92626d7.png 
pl = staticArcticMaps(book_ds.ice_thickness[5::12]-book_ds.ice_thickness[5::12].mean(axis=0), set_cbarlabel = "Sea ice thickness anomaly (m)", cmap="RdBu", vmin=-1.5, vmax=1.5, out_str='thickness_april_anomaly')
display(pl)
../_images/ef92e782b9e6cff28786d0560d42a3f1cac520c6c64d5bb224c5ca473374b47d.png

Monthly mean timeseries#

Next we’ll compute monthly means by averaging over all gridcells within a given region. We’ll use this to generate a lineplot to compare across the three winter seasons for each variable.

# Here is where we might also want to set a region mask, e.g. to avoid including some of the more uncertain data in the peripheral seas
innerArctic = [1,2,3,4,5,6]
book_ds_innerArctic = book_ds.where(book_ds.region_mask.isin(innerArctic))

# Uncomment out to set an additional ice type mask too and change the save_label accordingly (0 = FYI, 1 = MYI)
book_ds_innerArctic = book_ds_innerArctic.where(book_ds_innerArctic.ice_type==1)

save_label='Inner_Arctic_MYI'

print(static_winter_comparison_lineplot.__doc__) # Print docstring
Hide code cell output 
 Make a lineplot with markers comparing monthly mean data across winter seasons 
    
    Args: 
        da (xr.DataArray): data to plot and compute mean for; must contain "time" as a coordinate 
        years (list of str): list of years for which to plot data. 2020 would correspond to the winter season defined by start month 2020 - end month 2021 (default to all unique years in da)
        title (str, optional): title to give plot (default to no title) 
        set_ylabel (str, optional): prescribed y label string
        set_units (str, optional): prescribed y label unit string
        legend (bool): print legend
        savefig (bool): output figure
        save_label (str, optional): additional string for output
        figsize (tuple, optional): figure size to display in notebook (default to (5,3))
        start_month (str, optional): first month in winter (default to September)
        end_month (str, optional): second month in winter; this is the following calender year after start_month (default to April)
        force_complete_season (bool, optional): require that winter season returns data if and only if all months have data? i.e. if Sep and Oct have no data, return nothing even if Nov-Apr have data? (default to False) 
        loc_pos (int, optional): if greater than one use that, if not default to "best"

       Returns: 
           Figure displayed in notebook
        
    

static_winter_comparison_lineplot(book_ds_innerArctic.freeboard, years=years, start_month="Sep", figsize=(5,3), annotation='(a)', set_ylabel=r'Total freeboard (m)', save_label=save_label, loc_pos=4, legend=True)
../_images/99df8d7d0f4d56731cdb79cdad34bd766130c163f701f92f18cd5a1ac5639278.png 
static_winter_comparison_lineplot(book_ds_innerArctic.snow_depth, years=years, start_month="Sep", figsize=(5,3), annotation='(b)',set_ylabel='Snow depth (m)', save_label=save_label, legend=False)
../_images/4cde7ff31d4c6f64c79cb4ae44e07d37b32e6461c98f867fed5354d18a977393.png 
static_winter_comparison_lineplot(book_ds_innerArctic.snow_density, years=years, start_month="Sep", figsize=(5,3), annotation='(c)',set_ylabel=r'Snow density (kg/m$^3$)', save_label=save_label, legend=False)
../_images/200a83c4bf67e570823ec4c2355f8642a759f06256a97e82b2a9252f89daede6.png 
static_winter_comparison_lineplot(book_ds_innerArctic.ice_thickness, 
                                  da_unc = book_ds_innerArctic.ice_thickness_unc, 
                                  years=years, start_month="Sep", figsize=(5,3), set_ylabel='Sea ice thickness (m)', save_label=save_label, annotation='(d)', legend=False)
../_images/698d9b180cb782b19b4f5a932b21990ee4ddf94d65c47b39c0c03fe1c50b1e31.png 
static_winter_comparison_lineplot(book_ds_innerArctic.sea_ice_conc, years=years, start_month="Sep", figsize=(5,3), annotation='(e)', set_ylabel='Sea ice concentration', save_label=save_label, legend=False)
../_images/692e43375354eb63915c20e4d89d43b98752e0576b453a84c1c0d79b10f23f8a.png 
static_winter_comparison_lineplot(book_ds_innerArctic.ice_type, years=years, start_month="Sep", figsize=(5,3), annotation='(f)',set_ylabel='Multi-year ice fraction', save_label=save_label, legend=False)
../_images/e79308c4304e87818677eeee238dab2a967e9728fd6853c7fb25b54a1a224464.png

Let’s take things a bit further and set some conditions on the plots. Storing our data in an xarray dataset makes this very easy!

# Produce timeseries just using grid-cells of multi-year ice type (ice_type =1)
#static_winter_comparison_lineplot(out, title="Multiyear ice", years=years, figsize=(5,3), start_month="Nov", set_ylabel='Sea ice thickness (m)', save_label=save_label, legend=True)

# Produce timeseries just using grid-cells of first-year ice type (ice_type = 0)
# Use the flag at the start of this section to generate all plots based o the ice type classification
#out = book_ds_innerArctic.ice_thickness.where(book_ds_innerArctic.ice_type==0)
#static_winter_comparison_lineplot(out, years=years, figsize=(5,3), title="First-year ice", start_month="Nov", set_ylabel='Sea ice thickness (m)', save_label=save_label, legend=False)

Play around with some other conditions and see how things look!