Comparisons with CryoSat-2 and PIOMAS sea ice thickness estimates

Comparisons with CryoSat-2 and PIOMAS sea ice thickness estimates#

Summary: In this notebook, we compare our ICESat-2-derived sea ice thickness estimates with various estimates derived from ESA’s CryoSat-2 radar altimeter, a combined ICESat-2/CryoSat-2 product and PIOMAS. See the CryoSat-2 wrangling notebook for more insight into the different CryoSat-2 products.

Notes:

Doesn’t include spatial comparison stats, but would be pretty trivial to add that in - working on it!
Also lots more to be done in terms of assessing consistency between CS-2 and IS-2 and where/why that breaks down.

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

Import notebook dependencies#

import xarray as xr 
import pandas as pd
import numpy as np
import itertools
import pyproj 
from netCDF4 import Dataset
import scipy.interpolate 
# Helper function for reading the data from the bucket
from utils.read_data_utils import read_book_data  
from utils.plotting_utils import compute_gridcell_winter_means, interactiveArcticMaps, interactive_winter_mean_maps, interactive_winter_comparison_lineplot # Plotting

# Plotting dependencies
import cartopy.crs as ccrs
from textwrap import wrap
import hvplot.xarray
import holoviews as hv
import matplotlib.pyplot as plt
from matplotlib.axes import Axes
from cartopy.mpl.geoaxes import GeoAxes
GeoAxes._pcolormesh_patched = Axes.pcolormesh # Helps avoid some weird issues with the polar projection 
%config InlineBackend.figure_format = 'retina'
import matplotlib as mpl
mpl.rcParams['figure.dpi'] = 150 # Sets figure size in the notebook

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

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

Below, we’ll read in the aggreated IS-2/CS-2 data and restrict it to the Inner Arctic region.

# Read/download the IS-2 data that includes the various CryoSat-2 products too

book_ds = read_book_data(CS2=True)

# Restrict to the Inner Arctic domain 
innerArctic = [1,2,3,4,5,6]
book_ds = book_ds.where(book_ds.region_mask.isin(innerArctic))

# Set analysis time range

start_date = "Nov 2018"
end_date = "Apr 2021"
date_range = pd.date_range(start=start_date, end=end_date, freq='MS') # MS indicates a time frequency of start of the month
#date_range = date_range[((date_range.month <5) | (date_range.month > 8))]

Choose CS-2 product#

Choose from the following options: AWISMOS, CPOM, GSFC, UBRIS, KK. See the CryoSat-2 wrangling notebook for more insight into the different CryoSat-2 products.

cs2_var='AWISMOS'

Generate a new dataset using a common data coverage mask#

Approach is to mask primarily based on the IS2 data then PIOMAS and CS-2 if there is valid data in those datasets fro the given month (really just to stop cases when CS-2 data is completely missing from that month, PIOMAS has complete seasonal coverage).

Likely a much easier way to do this in xarray but couldn’t figure it out |

book_ds_common_mask_list=[]
for date in date_range:
    masked_month = book_ds.sel(time=date).copy(deep=True)
    # Mask everything where missing ICESat-2 data
    masked_month = masked_month.where(masked_month["ice_thickness_int"]>0.0)
    if np.isfinite(np.nanmean(masked_month["cs2_sea_ice_thickness_"+cs2_var])):
        # If there is valid data in the CS-2 dataset then mask all data using this dataset
        masked_month = masked_month.where(masked_month["cs2_sea_ice_thickness_"+cs2_var]>0.0) 
    if np.isfinite(np.nanmean(masked_month["piomas_ice_thickness"])):
        # If there is valid data in the PIOMAS dataset then mask all data using this dataset (I think this never really happens)
        masked_month = masked_month.where(masked_month["piomas_ice_thickness"]>0.0) 
    
    book_ds_common_mask_list.append(masked_month)

book_ds_common_mask = xr.concat(book_ds_common_mask_list, "time")

# Generate monthly means
mean_monthly_data = book_ds.mean(dim=("x","y"), skipna=True, keep_attrs=True)
mean_monthly_data_common_mask = book_ds_common_mask.mean(dim=("x","y"), skipna=True, keep_attrs=True)
#mean_monthly_gsfc = cs2_gsfc.mean(dim=("x","y"), skipna=True, keep_attrs=True)

# Holoviews requires us to generate individual line plots then combine later
is2_lineplot_p = mean_monthly_data.ice_thickness_int.hvplot.line(grid=True, label="IS2SITMOGR4", line_dash='dashed', color='gray', frame_width=700, frame_height=250) * mean_monthly_data.ice_thickness_int.hvplot.scatter(marker='x', color='gray', size=40) # Overlay scatter plot to add markers
is2_lineplot_p_cm = mean_monthly_data_common_mask.ice_thickness_int.hvplot.line(grid=True, label="IS2SITMOGR4 CM",  line_width=1.5,color='k', frame_width=700, frame_height=250) * mean_monthly_data_common_mask.ice_thickness_int.hvplot.scatter(marker='o', color='k', size=40) # Overlay scatter plot to add markers
cs2_lineplot_p = mean_monthly_data['cs2_sea_ice_thickness_'+cs2_var].hvplot.line(grid=True,  line_dash='dashed', label="CryoSat-2", color='g', frame_width=700, frame_height=250) * mean_monthly_data['cs2_sea_ice_thickness_'+cs2_var].hvplot.scatter(marker='x', color='g', size=40) # Overlay scatter plot to add markers
cs2_lineplot_p_cm = mean_monthly_data_common_mask['cs2_sea_ice_thickness_'+cs2_var].hvplot.line(grid=True, label="CryoSat-2 CM",  line_width=1,color='g', frame_width=700, frame_height=250) * mean_monthly_data_common_mask['cs2_sea_ice_thickness_'+cs2_var].hvplot.scatter(marker='o', color='g', size=40) # Overlay scatter plot to add markers
piomas_lineplot_p = mean_monthly_data.piomas_ice_thickness.hvplot.line(grid=True, label="PIOMAS", line_dash='dashed', color='c',  frame_width=700, frame_height=250) * mean_monthly_data.piomas_ice_thickness.hvplot.scatter(marker='x', color='c', size=30) # Overlay scatter plot to add markers
piomas_lineplot_p_cm = mean_monthly_data_common_mask.piomas_ice_thickness.hvplot.line(grid=True, label="PIOMAS CM", color='c', line_width=1,frame_width=700, frame_height=250) * mean_monthly_data_common_mask.piomas_ice_thickness.hvplot.scatter(marker='o', color='c', size=40) # Overlay scatter plot to add markers

r_squared_str = '%.02f' %(xr.corr(mean_monthly_data_common_mask.ice_thickness_int, mean_monthly_data['cs2_sea_ice_thickness_'+cs2_var]).values**2)
mb_str = '%.02f' %(np.mean(mean_monthly_data_common_mask.ice_thickness_int-mean_monthly_data_common_mask['cs2_sea_ice_thickness_'+cs2_var]).values)
sd_str = '%.02f' %(np.std(mean_monthly_data_common_mask.ice_thickness_int-mean_monthly_data_common_mask['cs2_sea_ice_thickness_'+cs2_var]).values)

r_squared_str_pmas = '%.02f' %(xr.corr(mean_monthly_data_common_mask.ice_thickness_int, mean_monthly_data['piomas_ice_thickness']).values**2)
mb_str_pmas = '%.02f' %(np.mean(mean_monthly_data_common_mask.ice_thickness_int-mean_monthly_data_common_mask['piomas_ice_thickness']).values)
sd_str_pmas = '%.02f' %(np.std(mean_monthly_data_common_mask.ice_thickness_int-mean_monthly_data_common_mask['piomas_ice_thickness']).values)


# Combine plots and display
comb_plot = (is2_lineplot_p*is2_lineplot_p_cm*cs2_lineplot_p*cs2_lineplot_p_cm*piomas_lineplot_p*piomas_lineplot_p_cm).opts(hv.opts.Layout(shared_axes=True, merge_tools=True))
comb_plot.opts(title=" IS2SITMOGR4 vs CS-2 ("+cs2_var+r") r²: "+r_squared_str+" "+mb_str+" +/- "+sd_str+r" (m)    IS2SITMOGR4 vs PMAS r²: "+r_squared_str_pmas+"  "+mb_str_pmas+" +/- "+sd_str_pmas+" (m)", fontsize=10)
#comb_plot.opts(annotate="Mean sea ice thickness (meters)")
comb_plot.opts(ylabel="Mean sea ice thickness (meters)")
#comb_plot.opts(annotate=r_str)
comb_plot.opts(xticks=book_ds.time.values[0::6])
display(comb_plot)
#hv.save(comb_plot, './figs/icesat-2_cryosat-2'+cs2_var+'.png', fmt='png', dpi=300)

# New version of the above, remove some the CM lines for clarity
# Generate monthly means
mean_monthly_data = book_ds.mean(dim=("x","y"), skipna=True, keep_attrs=True)
mean_monthly_data_common_mask = book_ds_common_mask.mean(dim=("x","y"), skipna=True, keep_attrs=True)
#mean_monthly_gsfc = cs2_gsfc.mean(dim=("x","y"), skipna=True, keep_attrs=True)

# Holoviews requires us to generate individual line plots then combine later
is2_lineplot_p = mean_monthly_data.ice_thickness_int.hvplot.line(grid=True, label="IS2SITMOGR4", line_dash='dashed', color='gray', frame_width=700, frame_height=250) * mean_monthly_data.ice_thickness_int.hvplot.scatter(marker='x', color='gray', size=40) # Overlay scatter plot to add markers
is2_lineplot_p_cm = mean_monthly_data_common_mask.ice_thickness_int.hvplot.line(grid=True, label="IS2SITMOGR4",  line_width=1.5,color='k', frame_width=700, frame_height=250) * mean_monthly_data_common_mask.ice_thickness_int.hvplot.scatter(marker='o', color='k', size=40) # Overlay scatter plot to add markers
cs2_lineplot_p = mean_monthly_data['cs2_sea_ice_thickness_'+cs2_var].hvplot.line(grid=True,  line_dash='dashed', label="CryoSat-2", color='g', frame_width=700, frame_height=250) * mean_monthly_data['cs2_sea_ice_thickness_'+cs2_var].hvplot.scatter(marker='x', color='g', size=40) # Overlay scatter plot to add markers
cs2_lineplot_p_cm = mean_monthly_data_common_mask['cs2_sea_ice_thickness_'+cs2_var].hvplot.line(grid=True, label="CryoSat-2/SMOS",  line_width=1,color='g', frame_width=700, frame_height=250) * mean_monthly_data_common_mask['cs2_sea_ice_thickness_'+cs2_var].hvplot.scatter(marker='o', color='g', size=40) # Overlay scatter plot to add markers
piomas_lineplot_p = mean_monthly_data.piomas_ice_thickness.hvplot.line(grid=True, label="PIOMAS (unmasked)", line_dash='dashed', color='c',  frame_width=700, frame_height=250) * mean_monthly_data.piomas_ice_thickness.hvplot.scatter(marker='x', color='c', size=30) # Overlay scatter plot to add markers
piomas_lineplot_p_cm = mean_monthly_data_common_mask.piomas_ice_thickness.hvplot.line(grid=True, label="PIOMAS", color='c', line_width=1,frame_width=700, frame_height=250) * mean_monthly_data_common_mask.piomas_ice_thickness.hvplot.scatter(marker='o', color='c', size=40) # Overlay scatter plot to add markers

r_squared_str = '%.02f' %(xr.corr(mean_monthly_data_common_mask.ice_thickness_int, mean_monthly_data['cs2_sea_ice_thickness_'+cs2_var]).values**2)
mb_str = '%.02f' %(np.mean(mean_monthly_data_common_mask.ice_thickness_int-mean_monthly_data_common_mask['cs2_sea_ice_thickness_'+cs2_var]).values)
sd_str = '%.02f' %(np.std(mean_monthly_data_common_mask.ice_thickness_int-mean_monthly_data_common_mask['cs2_sea_ice_thickness_'+cs2_var]).values)

r_squared_str_pmas = '%.02f' %(xr.corr(mean_monthly_data_common_mask.ice_thickness_int, mean_monthly_data['piomas_ice_thickness']).values**2)
mb_str_pmas = '%.02f' %(np.mean(mean_monthly_data_common_mask.ice_thickness_int-mean_monthly_data_common_mask['piomas_ice_thickness']).values)
sd_str_pmas = '%.02f' %(np.std(mean_monthly_data_common_mask.ice_thickness_int-mean_monthly_data_common_mask['piomas_ice_thickness']).values)


# Combine plots and display
#comb_plot = (is2_lineplot_p*is2_lineplot_p_cm*cs2_lineplot_p*cs2_lineplot_p_cm*piomas_lineplot_p*piomas_lineplot_p_cm).opts(hv.opts.Layout(shared_axes=True, merge_tools=True))
comb_plot = (is2_lineplot_p_cm*cs2_lineplot_p_cm*piomas_lineplot_p_cm*piomas_lineplot_p).opts(hv.opts.Layout(shared_axes=True, merge_tools=True))

comb_plot.opts(title=" IS2 vs CS-2/SMOS  r²: "+r_squared_str+", bias: "+mb_str+" +/- "+sd_str+r" (m)    IS2 vs PMAS  r²: "+r_squared_str_pmas+", bias: "+mb_str_pmas+" +/- "+sd_str_pmas+" (m)", fontsize=10)
#comb_plot.opts(annotate="Mean sea ice thickness (meters)")
comb_plot.opts(ylabel="Mean Arctic Ocean sea ice thickness (m)")
#comb_plot.opts(annotate=r_str)
comb_plot.opts(xticks=book_ds.time.values[0::6])
display(comb_plot)
#hv.save(comb_plot, './figs/icesat-2_cryosat-2'+cs2_var+'.png', fmt='png', dpi=300)

# Interactive map plot for a given month

date='2019-03'

# Use either all data or data on a common mask
data1 = book_ds['cs2_sea_ice_thickness_'+cs2_var].sel(time=date[0:7])
data2 = book_ds.ice_thickness_int.sel(time=date[0:7])
#data1 = book_ds_common_mask['cs2_sea_ice_thickness_'+cs2_var].sel(time=date[0:7])
#data2 = book_ds_common_mask.ice_thickness_int.sel(time=date[0:7])
pl1 = interactiveArcticMaps(data1, vmin=0, vmax=4, colorbar=False, title=date[0:7]+'  CS-2 ('+cs2_var+')')
pl2 = interactiveArcticMaps(data2, vmin=0, vmax=4, colorbar=True, clabel="sea ice thickness (m)", title='ICESat-2') 
diff_pl = interactiveArcticMaps((data2-data1), vmin=-1.5, vmax=1.5, cmap="coolwarm", clabel="difference", title="Difference (b - a)") 
data_pl = (pl1+pl2+diff_pl).opts(hv.opts.Layout(shared_axes=False, merge_tools=True))
#data_pl.opts(title="Winter "+str(years[i])+"-"+str(years[i]+1))
display(data_pl)

    