Comparisons with BGEP ULS sea ice draft estimates#

Summary: In this notebook, we produce comparisons of the gridded ICESat-2 sea ice thickness data with draft measurements obtained from Upward Looking Sonar moorings deployed in the Beaufort Sea.

Raw data link: https://www2.whoi.edu/site/beaufortgyre/data/mooring-data/2018-2021-mooring-data-from-the-bgep-project/

Notes:

  • With more years of data we hope to explore seasonal de-trending of the data. Current tools (e..g seasonal_decompose from statsmodel) require more years of data then currently provided by the ICESat2-/BGEP overlap period (2018 to 2021)

  • Some previous studies have removed zero draft values from ULS before undertaking these comparisons. Hard to know what is best here. ATL10 are available for passive microwave concentrations > 50%. Our gridded thickness data includes all this (baring some additional anomaly filters) so could also include reasonably large stretches of thin ice/open water too. I have toyed with included a ‘thickness where we have ice’ variable too, which I could easily include later, then we would just compare this to positive values of ULS draft?

  • Could require a minimum number of IS2 grid-cells, but would perhaps need to depend on the chosen comp_res. Note how in some September/October months there appears to be no data overlap, I included some maps at the end to highlight this more, the moorings are right on the edge of the ice pack in September 2019 for instance.

  • We do a fair bit of averaging to try and reduce noise from various sources, including sampling differences/representation errors. One more direct way of dealing with this would be to use the day_of_the_month information in the IS-2 thickness data to see what actual day the nearest grid-cells to the moorings best represent and comapre that to the daily ULS data. This leads to using the along-track data too, hard to know when to stop..

  • We could have made better use of Pandas in this analysis but I was a bit time pressured.

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 
from utils.read_data_utils import read_book_data # Helper function for reading the data from the bucket
from utils.plotting_utils import compute_gridcell_winter_means, interactiveArcticMaps, interactive_winter_mean_maps, interactive_winter_comparison_lineplot # Plotting
from scipy import stats
import datetime
# Plotting dependencies
import cartopy.crs as ccrs
from textwrap import wrap
import hvplot.pandas
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'] = 200 # Sets figure size in the notebook

# Remove warnings to improve display
import warnings 
warnings.filterwarnings('ignore') 
# Read/download the IS-2 thickness data 
# NG: don't need the CryoSat-2 data here but could if you want to explore those comparisons too...
book_ds = read_book_data()

# Issue in how region mask was assigned
book_ds['region_mask'] = book_ds['region_mask'].isel(time=0, drop=True)

# Estimate ice draft from the ICESat-2 data for more direct comparison with ULS draft measurements
# NB: Another option (not taken) is to simply multiply ULS drafts by 1.1 to convert to thickness (AWI do this)
# Include extra var = '_int' to use interpolated values, not recommended for validation purposes. 
int_str=''
book_ds['ice_draft'] = book_ds['ice_thickness'+int_str] - book_ds['freeboard'+int_str]+ book_ds['snow_depth'+int_str]
# Initialize map projection and project data to it

out_proj = 'EPSG:3411'
out_lons = book_ds.longitude.values
out_lats = book_ds.latitude.values

mapProj = pyproj.Proj("+init=" + out_proj)
xptsIS2, yptsIS2 = mapProj(out_lons, out_lats)

start_date = "Nov 2018"
end_date = "Apr 2021"
IS2_date_range = pd.date_range(start=start_date, end=end_date, freq='MS') # MS indicates a time frequency of start of the month
IS2_date_range = IS2_date_range[((IS2_date_range.month <5) | (IS2_date_range.month > 8))]
IS2_date_range_strs=[str(date.year)+'-%02d'%(date.month) for date in IS2_date_range]
# This is the value in meters of the aggreation length-scale for the ICESat-2 data
comp_res=25000
# Wrange ULS data
# Could really convert this to another data wrangling notebook and store the derived output as netcdfs

dataPathULS='./data/'

def get_ULS_dates(uls_mean_monthly_draft, uls_dates, date):
    #print(date)
    a = uls_mean_monthly_draft[uls_dates==date].values[0]
    return a

def get_uls(letter):
    if letter=='a':
        print('Mooring A (75.0 N, 150 W)')
        uls_x, uls_y = mapProj(-150., 75.)
    if letter=='b':
        print('Mooring B (78.4 N, 150.0 W)')
        uls_x, uls_y = mapProj(-150., 78.4)
    if letter=='d':
        print('Mooring D (74.0 N, 140.0 W)')
        uls_x, uls_y = mapProj(-140., 74.)
 
    uls = pd.read_csv(dataPathULS+'uls18'+letter+'_draft.dat', sep='\s+',names = ['date', 'time', 'draft'], header=2)
    utc_datetime_uls = pd.to_datetime(uls['date'], format='%Y%m%d')
    uls_mean_daily_draft = uls['draft'].groupby([utc_datetime_uls.dt.date]).mean() 

    uls_mean_monthly_draft = uls['draft'].groupby([utc_datetime_uls.dt.to_period('m')]).mean() 
    #print(uls_mean_monthly_draft)
    
    return uls_mean_daily_draft, uls_mean_monthly_draft, uls_x, uls_y

uls_mean_daily_draft_a, uls_mean_monthly_draft_a, uls_x_a, uls_y_a = get_uls('a')
uls_mean_daily_draft_b, uls_mean_monthly_draft_b, uls_x_b, uls_y_b = get_uls('b')
uls_mean_daily_draft_d, uls_mean_monthly_draft_d, uls_x_d, uls_y_d = get_uls('d')
Mooring A (75.0 N, 150 W)
Mooring B (78.4 N, 150.0 W)
Mooring D (74.0 N, 140.0 W)
def grid_IS2_nearby(date, uls_x, uls_y, res=75000):
    #print(date)
    IS2 = book_ds.ice_draft.sel(time=date)
    xptsIS2=IS2.xgrid.values
    yptsIS2=IS2.ygrid.values
    #dist = np.sqrt( (xptsIS2 - uls_x)**2 + (yptsIS2 - uls_y)**2 )
    
    #IS2_uls = IS2.where(dist<res).mean()
    
    #print('Number of valid IS-2 grid cells in month '+str(date)[0:7]+':', np.count_nonzero(~np.isnan(IS2.where(dist<res))))
    
    #Another option I first explored, coarsen the data then do nearest neighbor...provided similar results but above is more flexible.
    #if coarse_res>1:
    #Coarsen array by coarse_res in x/y directions (note that each grid-cell represents 25 km so 4 = 100 km)
    #    IS2 = IS2.coarsen(x=res, y=res, boundary='pad').mean()
    IS2_uls = scipy.interpolate.griddata((xptsIS2.flatten(),yptsIS2.flatten()), IS2.values.flatten(), (uls_x, uls_y), method = 'nearest')
    
    return IS2_uls


monthly_IS2_at_ULS_a  = [grid_IS2_nearby(date, uls_x_a, uls_y_a, res=comp_res) for date in IS2_date_range]
monthly_IS2_at_ULS_b  = [grid_IS2_nearby(date, uls_x_b, uls_y_b, res=comp_res) for date in IS2_date_range]
monthly_IS2_at_ULS_d  = [grid_IS2_nearby(date, uls_x_d, uls_y_d, res=comp_res) for date in IS2_date_range]

monthly_IS2_at_ULS_all = monthly_IS2_at_ULS_a+monthly_IS2_at_ULS_b+monthly_IS2_at_ULS_d
uls_dates=uls_mean_monthly_draft_a.index.astype(str)
uls_mean_monthly_draft_a_IS2_period  = [get_ULS_dates(uls_mean_monthly_draft_a, uls_dates, date) for date in IS2_date_range_strs]

uls_dates=uls_mean_monthly_draft_b.index.astype(str)
uls_mean_monthly_draft_b_IS2_period  = [get_ULS_dates(uls_mean_monthly_draft_b, uls_dates, date) for date in IS2_date_range_strs]

uls_dates=uls_mean_monthly_draft_d.index.astype(str)
uls_mean_monthly_draft_d_IS2_period  = [get_ULS_dates(uls_mean_monthly_draft_d, uls_dates, date) for date in IS2_date_range_strs]
                                        
uls_mean_monthly_draft_IS2_period = uls_mean_monthly_draft_a_IS2_period+uls_mean_monthly_draft_b_IS2_period+uls_mean_monthly_draft_d_IS2_period
# Seasonal plot

# Holoviews requires us to generate individual line plots then combine later
# Needed to reindex the pandas series to proper monthly datetime and offsetting to middle of the month

# Daily
uls_a_daily_lineplot_p = uls_mean_daily_draft_a.hvplot.line(grid=True, label="uls_daily", line_dash='solid', color='gray', frame_width=700, frame_height=200).opts( title='ULS A', xaxis=None, ylabel="Ice draft (meters)")
uls_b_daily_lineplot_p = uls_mean_daily_draft_b.hvplot.line(grid=True,  label="uls_daily", line_dash='solid', color='gray', frame_width=700, frame_height=200).opts(title='ULS B',xaxis=None, ylabel="Ice draft (meters)")
uls_d_daily_lineplot_p = uls_mean_daily_draft_d.hvplot.line(grid=True,  label="uls_daily", line_dash='solid', color='gray', frame_width=700, frame_height=200).opts(title='ULS D',ylabel="Ice draft (meters)")

# Monthly
ULS_date_range_a = pd.date_range(start=uls_mean_monthly_draft_a.index[0].to_timestamp().to_datetime64(), end=uls_mean_monthly_draft_a.index[-1].to_timestamp().to_datetime64(), freq='MS')+ datetime.timedelta(days=14)
uls_mean_monthly_draft_a_reindex = pd.Series(data=uls_mean_monthly_draft_a.values, index=ULS_date_range_a) 
uls_a_monthly_lineplot_p = uls_mean_monthly_draft_a_reindex.hvplot.line(grid=True, label="uls_monthly", line_dash='dashed', color='k', frame_width=700, frame_height=200).opts(ylabel="Ice draft (meters)")

ULS_date_range_b = pd.date_range(start=uls_mean_monthly_draft_b.index[0].to_timestamp().to_datetime64(), end=uls_mean_monthly_draft_b.index[-1].to_timestamp().to_datetime64(), freq='MS')+ datetime.timedelta(days=14)
uls_mean_monthly_draft_b_reindex = pd.Series(data=uls_mean_monthly_draft_b.values, index=ULS_date_range_b) 
uls_b_monthly_lineplot_p = uls_mean_monthly_draft_b_reindex.hvplot.line(grid=True, label="uls_monthly", line_dash='dashed', color='k', frame_width=700, frame_height=200).opts(ylabel="Ice draft (meters)")

ULS_date_range_d = pd.date_range(start=uls_mean_monthly_draft_d.index[0].to_timestamp().to_datetime64(), end=uls_mean_monthly_draft_d.index[-1].to_timestamp().to_datetime64(), freq='MS')+ datetime.timedelta(days=14)
uls_mean_monthly_draft_d_reindex = pd.Series(data=uls_mean_monthly_draft_d.values, index=ULS_date_range_d) 
uls_d_monthly_lineplot_p = uls_mean_monthly_draft_d_reindex.hvplot.line(grid=True, label="uls_monthly", line_dash='dashed', color='k', frame_width=700, frame_height=200).opts(ylabel="Ice draft (meters)")

# ICESat-2
is2_uls_mean_monthly_draft_a = pd.Series(data=monthly_IS2_at_ULS_a, index=IS2_date_range) 
is2_uls_a_monthly_lineplot_p = is2_uls_mean_monthly_draft_a.hvplot.scatter(grid=True, label="IS-2_monthly", marker='x', s=100, color='m', frame_width=700, frame_height=200)

is2_uls_mean_monthly_draft_b = pd.Series(data=monthly_IS2_at_ULS_b, index=IS2_date_range) 
is2_uls_b_monthly_lineplot_p = is2_uls_mean_monthly_draft_b.hvplot.scatter(grid=True, label="IS-2_monthly", marker='x', s=100, color='m', frame_width=700, frame_height=200)

is2_uls_mean_monthly_draft_d = pd.Series(data=monthly_IS2_at_ULS_d, index=IS2_date_range) 
is2_uls_d_monthly_lineplot_p = is2_uls_mean_monthly_draft_d.hvplot.scatter(grid=True, label="IS-2_monthly", marker='x', s=100, color='m', frame_width=700, frame_height=200)

# Combine plots and display
comb_plot = ((uls_a_daily_lineplot_p*uls_a_monthly_lineplot_p*is2_uls_a_monthly_lineplot_p).opts(legend_position='top_left') 
             + (uls_b_daily_lineplot_p*uls_b_monthly_lineplot_p*is2_uls_b_monthly_lineplot_p).opts(legend_position='top_left')
             + (uls_d_daily_lineplot_p*uls_d_monthly_lineplot_p*is2_uls_d_monthly_lineplot_p).opts(legend_position='top_left')).cols(1)

display(comb_plot)
#hv.save(comb_plot, './figs/BGEP_seasonal_'+str(comp_res)+int_str+'.png', fmt='png', dpi=300)
# Validation analysis

# ULS A
mask = ~np.isnan(monthly_IS2_at_ULS_a)
res = stats.linregress(np.array(monthly_IS2_at_ULS_a)[mask], np.array(uls_mean_monthly_draft_a_IS2_period)[mask])

r_str_a = '%.02f' %(res[2]**2)
mb_str_a = '%.02f' %(np.nanmean(np.array(monthly_IS2_at_ULS_a)-np.array(uls_mean_monthly_draft_a_IS2_period)))
sd_str_a = '%.02f' %(np.nanstd(np.array(monthly_IS2_at_ULS_a)-np.array(uls_mean_monthly_draft_a_IS2_period)))

# ULS B
mask = ~np.isnan(monthly_IS2_at_ULS_b)
res = stats.linregress(np.array(monthly_IS2_at_ULS_b)[mask], np.array(uls_mean_monthly_draft_b_IS2_period)[mask])

r_str_b = '%.02f' %(res[2]**2)
mb_str_b = '%.02f' %(np.nanmean(np.array(monthly_IS2_at_ULS_b)-np.array(uls_mean_monthly_draft_b_IS2_period)))
sd_str_b = '%.02f' %(np.nanstd(np.array(monthly_IS2_at_ULS_b)-np.array(uls_mean_monthly_draft_b_IS2_period)))

# ULS D
mask = ~np.isnan(monthly_IS2_at_ULS_d)
res = stats.linregress(np.array(monthly_IS2_at_ULS_d)[mask], np.array(uls_mean_monthly_draft_d_IS2_period)[mask])

r_str_d = '%.02f' %(res[2]**2)
mb_str_d = '%.02f' %(np.nanmean(np.array(monthly_IS2_at_ULS_d)-np.array(uls_mean_monthly_draft_d_IS2_period)))
sd_str_d = '%.02f' %(np.nanstd(np.array(monthly_IS2_at_ULS_d)-np.array(uls_mean_monthly_draft_d_IS2_period)))

# ULS ALL
mask = ~np.isnan(monthly_IS2_at_ULS_all)
res = stats.linregress(np.array(monthly_IS2_at_ULS_all)[mask], np.array(uls_mean_monthly_draft_IS2_period)[mask])

r_str_all = '%.02f' %(res[2]**2)
mb_str_all = '%.02f' %(np.nanmean(np.array(monthly_IS2_at_ULS_all)-np.array(uls_mean_monthly_draft_IS2_period)))
sd_str_all = '%.02f' %(np.nanstd(np.array(monthly_IS2_at_ULS_all)-np.array(uls_mean_monthly_draft_IS2_period)))
# Validation plot

plt.rc('font',**{'family':'sans-serif','sans-serif':['Arial']})
fig=plt.figure(figsize=(4, 4))
plt.scatter(monthly_IS2_at_ULS_a, uls_mean_monthly_draft_a_IS2_period, color='b')
plt.scatter(monthly_IS2_at_ULS_b, uls_mean_monthly_draft_b_IS2_period, color='tab:red')
plt.scatter(monthly_IS2_at_ULS_d, uls_mean_monthly_draft_d_IS2_period, color='tab:orange')
#plt.scatter(IS2_ULS_all, ULS_MONTHLY_IS2_all, color='k')

plt.annotate(r"BGEP A  r$^2$: "+r_str_a+"  "+mb_str_a+" +/- "+sd_str_a+" (m)", color='b', xy=(0.02, 0.98),xycoords='axes fraction', horizontalalignment='left', verticalalignment='top', fontsize=9)
plt.annotate(r"BGEP B  r$^2$: "+r_str_b+"  "+mb_str_b+" +/- "+sd_str_b+" (m)", color='tab:red', xy=(0.02, 0.92),xycoords='axes fraction', horizontalalignment='left', verticalalignment='top', fontsize=9)
plt.annotate(r"BGEP D  r$^2$: "+r_str_d+"  "+mb_str_d+" +/- "+sd_str_d+" (m)", color='tab:orange', xy=(0.02, 0.86),xycoords='axes fraction', horizontalalignment='left', verticalalignment='top', fontsize=9)
plt.annotate(r"ALL  r$^2$: "+r_str_all+"  "+mb_str_all+" +/- "+sd_str_all+" (m)", color='k', xy=(0.02, 0.80),xycoords='axes fraction', horizontalalignment='left', verticalalignment='top', fontsize=9)
        
plt.ylabel('BGEP/ULS ice draft (m)')
plt.xlabel('IS-2 derived ice draft (m)')
plt.xlim([0, 2.5])
plt.ylim([0, 2.5])
ax=plt.gca()
lims = [
    np.min([ax.get_xlim(), ax.get_ylim()]),  # min of both axes
    np.max([ax.get_xlim(), ax.get_ylim()]),  # max of both axes
]

# now plot both limits against each other
ax.plot(lims, lims, 'k--', alpha=0.5, zorder=0)
ax.set_aspect('equal')
ax.set_xlim(lims)
ax.set_ylim(lims)
plt.savefig('./figs/BGEP_scatter_res'+str(comp_res)+int_str+'_nn.pdf', dpi=200)
plt.show()
../_images/39391ec09a05d1aa9047407260fbd7a04715490d3c4a93245e304b0dad221f90.png

As a futher check, take a look at some maps to assess IS-2 coverage comapred to ULS locations, especially in the early season#

uls_x_a, uls_y_a = mapProj(-150., 75.)
uls_x_b, uls_y_b = mapProj(-150., 78.4)
uls_x_d, uls_y_d = mapProj(-140., 74.)

date = "Sep 2019"
var = "ice_thickness"+int_str
data_one_month = book_ds[var].sel(time=date)
p = data_one_month.plot(x="longitude", y="latitude", # horizontal coordinates 
                        vmin=0, vmax=4, # min and max on the colorbar 
                        subplot_kws={'projection':ccrs.NorthPolarStereo(central_longitude=-45)}, # Set the projection 
                        transform=ccrs.PlateCarree())

plt.plot(uls_x_a, uls_y_a, '.', color='b', transform=ccrs.NorthPolarStereo(central_longitude=-45))
plt.plot(uls_x_b, uls_y_b, '.', color='m', transform=ccrs.NorthPolarStereo(central_longitude=-45))
plt.plot(uls_x_d, uls_y_d, '.', color='k', transform=ccrs.NorthPolarStereo(central_longitude=-45))

plt.text(uls_x_a, uls_y_a, 'A', color='b', transform=ccrs.NorthPolarStereo(central_longitude=-45))
plt.text(uls_x_b, uls_y_b, 'B', color='m', transform=ccrs.NorthPolarStereo(central_longitude=-45))
plt.text(uls_x_d, uls_y_d, 'D', color='k', transform=ccrs.NorthPolarStereo(central_longitude=-45))


p.axes.coastlines() # Add coastlines 
plt.show()
../_images/ba932893f4564142b191c95f43cfa15676204f966d2e9aa6e2d4e3890eb1fe3a.png