All-season Arctic sea ice thickness (preliminary, using new summer ICESat-2 ice thickness estimates)#
Summary: In this notebook, we provide a preliminary look at all-season thickness estimates from ICESat-2, including comparisons with a CryoSat-2 thickness product and a preliminary ICESat-2/CryoSat-2 snow depth product.
Author: Alek Petty
Version history: Version 1 (05/01/2025)
# Regular Python library imports
import xarray as xr
import numpy as np
import pandas as pd
import pyproj
import scipy.interpolate
import matplotlib.pyplot as plt
import glob
from datetime import datetime
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import os
# Helper functions for reading the data from the bucket and plotting
from utils.read_data_utils import read_IS2SITMOGR4, read_book_data
from utils.plotting_utils import static_winter_comparison_lineplot, staticArcticMaps, interactiveArcticMaps, compute_gridcell_winter_means, interactive_winter_comparison_lineplot # Plotting utils
from extra_funcs import read_IS2SITMOGR4S, regrid_ubris_to_is2, get_cs2is2_snow, apply_interpolation_time
# Plotting dependencies
#%config InlineBackend.figure_format = 'retina'
import matplotlib as mpl
# Remove warnings to improve display
import warnings
warnings.filterwarnings('ignore')
# Get the current working directory
current_directory = os.getcwd()
# Set some plotting parameters
mpl.rcParams.update({
"text.usetex": False, # Use LaTeX for rendering
"font.family": "sans-serif",
"lines.linewidth": 0.8,
"font.size": 9,
#"lines.alpha": 0.8,
"axes.labelsize": 8,
"xtick.labelsize": 8,
"ytick.labelsize": 8,
"legend.fontsize": 7
})
mpl.rcParams['font.sans-serif'] = ['Arial']
mpl.rcParams['figure.dpi'] = 300
# add '_int' if we want to mask using the interpolated/smoothed IS-2 data.
# Annoyingly we dont have the interpolated data for the IS-2 SM-LG data in winter, so I default to not using that so all variables are consistent in coverage.
int_str='_int'
# Test loading the saved file
IS2_CS2_allseason = xr.open_dataset('./data/book_data_allseason.nc')
print("Successfully loaded all-season wrangled dataset")
Successfully loaded all-season wrangled dataset
# Get some map proj info needed for later functions
out_proj = 'EPSG:3411'
mapProj = pyproj.Proj("+init=" + out_proj)
xIS2 = IS2_CS2_allseason.x.values
yIS2 = IS2_CS2_allseason.y.values
xptsIS2, yptsIS2 = np.meshgrid(xIS2, yIS2)
out_lons = IS2_CS2_allseason.longitude.values
out_lats = IS2_CS2_allseason.latitude.values
IS2_CS2_allseason
<xarray.Dataset>
Dimensions: (time: 33, y: 448, x: 304)
Coordinates:
* time (time) datetime64[ns] 2018-11-15 ... 2021...
* x (x) float32 -3.838e+06 ... 3.738e+06
* y (y) float32 5.838e+06 ... -5.338e+06
longitude (y, x) float32 168.3 168.1 ... -10.18 -9.999
latitude (y, x) float32 31.1 31.2 ... 34.58 34.47
Data variables: (12/47)
crs (time) float64 ...
ice_thickness_sm (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_sm (time, y, x) float32 ...
... ...
ice_density_j22_int (time, y, x) float32 ...
ice_thickness_j22_int (time, y, x) float32 ...
ice_thickness_cs2_ubris (time, y, x) float64 ...
cs2_sea_ice_type_UBRIS (time, y, x) float64 ...
cs2_sea_ice_density_UBRIS (time, y, x) float64 ...
cs2is2_snow_depth (time, y, x) float64 ...
Attributes:
contact: Alek Petty (akpetty@umd.edu)
description: Aggregated IS2SITMOGR4 summer V0 dataset.
history: Created 20/12/23# Plot the region mask just as a quick sanity check
im = IS2_CS2_allseason.region_mask.plot(x="longitude", y="latitude", transform=ccrs.PlateCarree(), cmap="viridis", zorder=3,
cbar_kwargs={'pad':0.01,'shrink': 0.8,'extend':'both', 'label':'Region Mask', 'location':'left'},
vmin=0, vmax=25,
subplot_kws={'projection':ccrs.NorthPolarStereo(central_longitude=-45)})
ax = im.axes
ax.coastlines()
ax.set_extent([-179, 179, 68, 90], crs=ccrs.PlateCarree())
ax.gridlines(draw_labels=True)
plt.show()
# Inner Arctic domain (Central_Arctic Beaufort_Sea Chukchi_Sea East_Siberian_Sea Laptev_Sea Kara_Sea)
innerArctic = [1,2,3,4,5,6]
# Inner Arcic plus Barents and Greenland
#innerArctic = [1,2,3,4,5,6,7,8]
# Plot the mean summer thickness data just for broad context/check all is good.
# Also add BGEP/ULS locations for reference
# Select only the summer months
summer_months = [5, 6, 7, 8] # May=5, June=6, July=7, August=8
IS2_CS2_summer = IS2_CS2_allseason.sel(time=IS2_CS2_allseason.time.dt.month.isin(summer_months))
# Calculate the mean ice thickness over the summer months
mean_thickness = IS2_CS2_summer.ice_thickness_sm_int.mean(dim='time')
plt.figure(figsize=(4, 4))
# Plot the region mask
im =mean_thickness.plot(x="longitude", y="latitude", transform=ccrs.PlateCarree(), cmap="viridis", zorder=3,
cbar_kwargs={'pad':0.01,'shrink': 0.7,'extend':'both', 'label':'Sea ice thickness (m)', 'location':'left'},
vmin=0, vmax=4,
subplot_kws={'projection':ccrs.NorthPolarStereo(central_longitude=-45)})
ax = im.axes
ax.coastlines()
ax.set_extent([-179, 179, 65, 90], crs=ccrs.PlateCarree())
ax.gridlines(draw_labels=True, linewidth=0.3, zorder=5)
# Coordinates for BGEP/ULS in the Beaufort Sea
beaufort_sea_lona = -150
beaufort_sea_lata = 75
beaufort_sea_lonb = -150.5
beaufort_sea_latb = 77.6
beaufort_sea_lond = -140.2
beaufort_sea_latd = 73.6
# Add markers for BGEP/ULS
ax.scatter(beaufort_sea_lona, beaufort_sea_lata, color='red', s=50, transform=ccrs.PlateCarree(), zorder=5)
ax.scatter(beaufort_sea_lonb, beaufort_sea_latb, color='blue', s=50, transform=ccrs.PlateCarree(), zorder=5)
ax.scatter(beaufort_sea_lond, beaufort_sea_latd, color='black', s=50, transform=ccrs.PlateCarree(), zorder=5)
ax.text(beaufort_sea_lona - 2, beaufort_sea_lata-3, 'BGEP/ULS', transform=ccrs.PlateCarree(), fontsize=9, color='black')
# Create a binary region mask
region_mask = IS2_CS2_allseason.region_mask.isin(innerArctic).astype(int)
# Cartopy contour issue with NPS plots so split the data across the meridian
lon = mean_thickness.longitude
lat = mean_thickness.latitude
mask_west = lon < 0
mask_east = lon >= 0
ax.contour(lon.where(mask_west), lat, region_mask.where(mask_west), levels=[0.5], colors='m', linewidths=0.9, transform=ccrs.PlateCarree(), zorder=5)
ax.contour(lon.where(mask_east), lat, region_mask.where(mask_east), levels=[0.5], colors='m', linewidths=0.9, transform=ccrs.PlateCarree(), zorder=5)
plt.subplots_adjust(left=0.02, right=0.95, top=0.95, bottom=0.05)
plt.savefig('./figs/summer_mean_is2.png', dpi=300, facecolor="white")
plt.show()
# Plot the individual monthly IS-2 summer thickness data
# Select only the summer months
summer_months = [5, 6, 7, 8] # May=5, June=6, July=7, August=8
IS2_CS2_summer = IS2_CS2_allseason.sel(time=IS2_CS2_allseason.time.dt.month.isin(summer_months))
# Group by year for plotting purposes
grouped_by_year = IS2_CS2_summer.groupby('time.year')
fig, axes = plt.subplots(nrows=len(grouped_by_year.groups), ncols=len(summer_months),
subplot_kw={'projection': ccrs.NorthPolarStereo(central_longitude=-45)},
figsize=(7, 6.1))
# Iterate over each year
for year_idx, (year, year_data) in enumerate(grouped_by_year):
#print(year_idx)# Group by month within each year
grouped_by_month = year_data.groupby('time.month')
for month_idx, (month, month_data) in enumerate(grouped_by_month):
ax = axes[year_idx, month_idx]
im = month_data.isel(time=0)['ice_thickness_sm_int'].plot(ax=ax, x="longitude", y="latitude",
transform=ccrs.PlateCarree(), cmap="viridis",
vmin=0, vmax=5, add_colorbar=False)
ice_conc = month_data.isel(time=0)['sea_ice_conc']
lons=ice_conc.longitude
lats=ice_conc.latitude
#print(lons.shape)
mask_pos = lons >= 0
mask_neg = lons < 0
# Plot contours for positive longitudes
cs_pos = ax.contour(lons.where(mask_pos), lats, ice_conc.where(mask_pos),
levels=[0.5], colors='c', linewidths=0.7,
transform=ccrs.PlateCarree())
# Plot contours for negative longitudes
cs_neg = ax.contour(lons.where(mask_neg), lats, ice_conc.where(mask_neg),
levels=[0.5], colors='c', linewidths=0.7 ,
transform=ccrs.PlateCarree())
ax.set_extent([-179, 179, 66, 90], crs=ccrs.PlateCarree())
ax.text(0.98, 0.9, f'{year}-{month:02d}', transform=ax.transAxes, va='bottom', ha='right')
ax.set_title('')
ax.coastlines(linewidth=0.15, color = 'black', zorder = 2) # Coastlines
ax.add_feature(cfeature.LAND, color ='0.95', zorder = 1) # Land
#ax.add_feature(cfeature.LAKES, color = 'grey', zorder = 5) # Lakes
ax.gridlines(draw_labels=False, linewidth=0.25, color='gray', alpha=0.7, linestyle='--', zorder=3) # Gridlines
axes[2, 3].set_visible(False)
# Add a colorbar at the bottom
cbar = fig.colorbar(im, ax=axes, orientation='horizontal', extend='both', fraction=0.03, pad=0.1)
cbar.set_label('Sea ice thickness (m)')
plt.subplots_adjust(left=0.02, right=0.98, top=0.98, bottom=0.15, wspace=0.01, hspace=0.03)
plt.savefig('./figs/maps_is2_summer_months_int.png', dpi=300, facecolor="white")
plt.show()
# Plot the individual monthly CS-2 summer thickness data
fig, axes = plt.subplots(nrows=len(grouped_by_year.groups), ncols=len(summer_months),
subplot_kw={'projection': ccrs.NorthPolarStereo(central_longitude=-45)},
figsize=(7, 6.1))
# Iterate over each year
for year_idx, (year, year_data) in enumerate(grouped_by_year):
#print(year_idx)# Group by month within each year
grouped_by_month = year_data.groupby('time.month')
for month_idx, (month, month_data) in enumerate(grouped_by_month):
ax = axes[year_idx, month_idx]
im = month_data.isel(time=0).ice_thickness_cs2_ubris.plot(ax=ax, x="longitude", y="latitude",
transform=ccrs.PlateCarree(), cmap="viridis",
vmin=0, vmax=5, add_colorbar=False)
ax.set_extent([-179, 179, 66, 90], crs=ccrs.PlateCarree())
ax.text(0.98, 0.9, f'{year}-{month:02d}', transform=ax.transAxes, va='bottom', ha='right')
ax.set_title('')
ax.coastlines(linewidth=0.15, color = 'black', zorder = 2) # Coastlines
ax.add_feature(cfeature.LAND, color ='0.95', zorder = 1) # Land
#ax.add_feature(cfeature.LAKES, color = 'grey', zorder = 5) # Lakes
ax.gridlines(draw_labels=False, linewidth=0.25, color='gray', alpha=0.7, linestyle='--', zorder=3) # Gridlines
axes[2, 3].set_visible(False)
# Add a colorbar at the bottom
cbar = fig.colorbar(im, ax=axes, orientation='horizontal', extend='both', fraction=0.03, pad=0.1)
cbar.set_label('Sea ice thickness (m)')
plt.subplots_adjust(left=0.02, right=0.98, top=0.98, bottom=0.15, wspace=0.01, hspace=0.03)
plt.savefig('./figs/maps_cs2_summer_months.png', dpi=300, facecolor="white")
plt.show()
# Plot the individual monthly IS-2 vs CS-2 summer thickness data differences
fig, axes = plt.subplots(nrows=len(grouped_by_year.groups), ncols=len(summer_months),
subplot_kw={'projection': ccrs.NorthPolarStereo(central_longitude=-45)},
figsize=(7, 6.1))
# Iterate over each year
for year_idx, (year, year_data) in enumerate(grouped_by_year):
#print(year_idx)# Group by month within each year
grouped_by_month = year_data.groupby('time.month')
for month_idx, (month, month_data) in enumerate(grouped_by_month):
ax = axes[year_idx, month_idx]
im = (month_data.isel(time=0).ice_thickness_sm_int - month_data.isel(time=0).ice_thickness_cs2_ubris).plot(ax=ax, x="longitude", y="latitude",
transform=ccrs.PlateCarree(), cmap="RdBu",
vmin=-3, vmax=3, add_colorbar=False)
ax.set_extent([-179, 179, 66, 90], crs=ccrs.PlateCarree())
ax.text(0.98, 0.9, f'{year}-{month:02d}', transform=ax.transAxes, va='bottom', ha='right')
ax.set_title('')
ax.coastlines(linewidth=0.15, color = 'black', zorder = 2) # Coastlines
ax.add_feature(cfeature.LAND, color ='0.95', zorder = 1) # Land
#ax.add_feature(cfeature.LAKES, color = 'grey', zorder = 5) # Lakes
ax.gridlines(draw_labels=False, linewidth=0.25, color='gray', alpha=0.7, linestyle='--', zorder=3) # Gridlines
axes[2, 3].set_visible(False)
# Add a colorbar at the bottom
cbar = fig.colorbar(im, ax=axes, orientation='horizontal', extend='both', fraction=0.03, pad=0.1)
cbar.set_label('Sea ice thickness difference (m)')
plt.subplots_adjust(left=0.02, right=0.98, top=0.98, bottom=0.15, wspace=0.01, hspace=0.03)
plt.savefig('./figs/maps_is2_cs2_summer_anoms.png', dpi=300, facecolor="white")
plt.show()
Time series analysis#
Aim here is to generate time-series plots of the data. Note the extra region masking and also common data masking, where we only analyze grid-cells with consistent data across the given variables in the given month.
Play around with different options here as desired (a few opiotns provided in the commented out lines). Note that later we apply a different masking method, only looking at data where both variables have data across ALL months (not just for a given month).
# Again we can apply a region or ice type mask as desired
IS2_CS2_allseason = IS2_CS2_allseason.where(IS2_CS2_allseason.region_mask.isin(innerArctic))
# Could also apply ice type masking here too if desired
#IS2_CS2_allseason = IS2_CS2_allseason.where(IS2_CS2_allseason.region_mask.isin([1]))
#IS2_CS2_allseason = IS2_CS2_allseason.where(IS2_CS2_allseason.ice_type==1)
def mask_and_concat_ice_thickness(IS2_CS2_allseason, start_date="Oct 2018", end_date="July 2021", int_str=''):
"""
Masks and concatenates ice thickness data over a specified date range.
Args:
IS2_CS2_allseason (xarray.Dataset): The dataset containing ice thickness data.
start_date (str): The start date for the date range (default is "Oct 2018").
end_date (str): The end date for the date range (default is "Apr 2023").
Returns:
xarray.Dataset: Concatenated dataset with applied masks.
"""
date_range = pd.date_range(start=start_date, end=end_date, freq='MS') + pd.Timedelta(days=14)
IS2_CS2_allseason_region_common_list = []
vars_to_check = ["ice_thickness_cs2_ubris", "ice_thickness_sm"+int_str]
for date in date_range:
#print(date)
try:
#print('try')
monthly_data_to_mask = IS2_CS2_allseason.sel(time=date).copy(deep=True)
for var_check in vars_to_check:
#print(var_check)
if np.isfinite(np.nanmean(monthly_data_to_mask[var_check])):
monthly_data_to_mask = monthly_data_to_mask.where(monthly_data_to_mask[var_check] > 0.0)
IS2_CS2_allseason_region_common_list.append(monthly_data_to_mask)
except:
print('except')
continue
return xr.concat(IS2_CS2_allseason_region_common_list, "time")
# Apply function to mask and concatenate ice thickness data where all variables have valid data for the given month
IS2_CS2_allseason_region_common = mask_and_concat_ice_thickness(IS2_CS2_allseason, int_str=int_str)
except
# plot monthly maps again to see how things look with the common mask applied (mostly that CS-2 data is removed where we have no IS-2 data)
summer_months = [5, 6, 7, 8] # May=5, June=6, July=7, August=8
IS2_CS2_summer = IS2_CS2_allseason_region_common.sel(time=IS2_CS2_allseason_region_common.time.dt.month.isin(summer_months))
# Group by year
grouped_by_year = IS2_CS2_summer.groupby('time.year')
# Create a figure with subplots
fig, axes = plt.subplots(nrows=len(grouped_by_year.groups), ncols=len(summer_months),
subplot_kw={'projection': ccrs.NorthPolarStereo(central_longitude=-45)},
figsize=(7, 6.1))
for year_idx, (year, year_data) in enumerate(grouped_by_year):
#print(year_idx)# Group by month within each year
grouped_by_month = year_data.groupby('time.month')
for month_idx, (month, month_data) in enumerate(grouped_by_month):
ax = axes[year_idx, month_idx]
im = month_data.isel(time=0)['ice_thickness_sm'+int_str].plot(ax=ax, x="longitude", y="latitude",
transform=ccrs.PlateCarree(), cmap="viridis",
vmin=0, vmax=5, add_colorbar=False)
ax.set_extent([-179, 179, 66, 90], crs=ccrs.PlateCarree())
ax.text(0.98, 0.9, f'{year}-{month:02d}', transform=ax.transAxes, va='bottom', ha='right')
ax.set_title('')
ax.coastlines(linewidth=0.15, color = 'black', zorder = 2) # Coastlines
ax.add_feature(cfeature.LAND, color ='0.95', zorder = 1) # Land
#ax.add_feature(cfeature.LAKES, color = 'grey', zorder = 5) # Lakes
ax.gridlines(draw_labels=False, linewidth=0.25, color='gray', alpha=0.7, linestyle='--', zorder=3) # Gridlines
axes[2, 3].set_visible(False)
cbar = fig.colorbar(im, ax=axes, orientation='horizontal', extend='both', fraction=0.03, pad=0.1)
cbar.set_label('IS-2 sea ice thickness (m)')
plt.subplots_adjust(left=0.02, right=0.98, top=0.98, bottom=0.15, wspace=0.01, hspace=0.03)
plt.savefig('./figs/maps_is2_summer_months_cm'+int_str+'.png', dpi=300, facecolor="white")
plt.show()
# Do same for CS-2 data
IS2_CS2_summer = IS2_CS2_allseason_region_common.sel(time=IS2_CS2_allseason_region_common.time.dt.month.isin(summer_months))
# Group by year
grouped_by_year = IS2_CS2_summer.groupby('time.year')
fig, axes = plt.subplots(nrows=len(grouped_by_year.groups), ncols=len(summer_months),
subplot_kw={'projection': ccrs.NorthPolarStereo(central_longitude=-45)},
figsize=(7, 6.1))
for year_idx, (year, year_data) in enumerate(grouped_by_year):
#print(year_idx)# Group by month within each year
grouped_by_month = year_data.groupby('time.month')
for month_idx, (month, month_data) in enumerate(grouped_by_month):
ax = axes[year_idx, month_idx]
im = month_data.isel(time=0).ice_thickness_cs2_ubris.plot(ax=ax, x="longitude", y="latitude",
transform=ccrs.PlateCarree(), cmap="viridis",
vmin=0, vmax=5, add_colorbar=False)
ax.set_extent([-179, 179, 66, 90], crs=ccrs.PlateCarree())
ax.text(0.98, 0.9, f'{year}-{month:02d}', transform=ax.transAxes, va='bottom', ha='right')
ax.set_title('')
ax.coastlines(linewidth=0.15, color = 'black', zorder = 2) # Coastlines
ax.add_feature(cfeature.LAND, color ='0.95', zorder = 1) # Land
#ax.add_feature(cfeature.LAKES, color = 'grey', zorder = 5) # Lakes
ax.gridlines(draw_labels=False, linewidth=0.25, color='gray', alpha=0.7, linestyle='--', zorder=3) # Gridlines
axes[2, 3].set_visible(False)
cbar = fig.colorbar(im, ax=axes, orientation='horizontal', extend='both', fraction=0.03, pad=0.1)
cbar.set_label('CS-2 sea ice thickness (m)')
plt.subplots_adjust(left=0.02, right=0.98, top=0.98, bottom=0.15, wspace=0.01, hspace=0.03)
plt.savefig('./figs/maps_cs2_summer_months_cm'+int_str+'.png', dpi=300, facecolor="white")
plt.show()
# Generate Arctic basin-mean time series
IS2_CS2_allseason_region_means = IS2_CS2_allseason_region_common.mean(dim=("x","y"), skipna=True, keep_attrs=True).resample(skipna=True, time='1M', label='left',loffset='15D').asfreq()
# Commented out as I now do this earlier, but nif not...
# Drop later data as no public SM-LG data. Drop if you want to show the extended IS-2 data on the plot.
#end_date = "2021-07-21"
#IS2_CS2_allseason_region_means = IS2_CS2_allseason_region_means.sel(time=slice(None, end_date))
# Plot Arctic basin-mean time series
import matplotlib.dates as mdates
fig, axes = plt.subplots(nrows=6, ncols=1, figsize=(6.8, 7.6))
# Calculate correlation coefficient and mean bias for both IS-2 products
corr_sm = np.corrcoef(IS2_CS2_allseason_region_means["ice_thickness_cs2_ubris"],
IS2_CS2_allseason_region_means["ice_thickness_sm"+int_str])[0,1]
bias_sm = np.mean(IS2_CS2_allseason_region_means["ice_thickness_sm"+int_str] -
IS2_CS2_allseason_region_means["ice_thickness_cs2_ubris"])
# Plot sea ice thickness
axes[0].plot(IS2_CS2_allseason_region_means.time, IS2_CS2_allseason_region_means["ice_thickness"+int_str], 'ko-', markersize=3, label='IS-2/NSIM')
axes[0].plot(IS2_CS2_allseason_region_means.time, IS2_CS2_allseason_region_means["ice_thickness_sm"+int_str], 'bo-', markersize=3, label='IS-2/SM-LG')
axes[0].plot(IS2_CS2_allseason_region_means.time, IS2_CS2_allseason_region_means["ice_thickness_cs2_ubris"], 'ro-', markersize=3, label='CS-2/SM-LG')
axes[0].set_ylabel('Sea ice thickness (m)')
axes[0].legend(loc='upper right', ncol=4, frameon=False)
axes[0].set_yticks([1, 2, 3])
axes[0].set_ylim([0.4,3.9])
axes[0].grid(axis='y') # Only show grid lines for the x-axis
axes[0].tick_params(axis='x', which='both', bottom=False, top=False, labelbottom=False) # Remove x-axis ticks and labels
# Add statistics as text annotations
axes[0].text(0.01, 0.6, f'r: {corr_sm:.2f}, MB: {bias_sm:.2f} m', transform=axes[0].transAxes, color='blue')
# Plot thickness anomaly relative to CS-2_UIT
axes[1].plot(IS2_CS2_allseason_region_means.time, IS2_CS2_allseason_region_means["ice_thickness_sm"+int_str] - IS2_CS2_allseason_region_means["ice_thickness_cs2_ubris"], 'bo-', markersize=3, label='')
axes[1].plot(IS2_CS2_allseason_region_means.time, IS2_CS2_allseason_region_means["ice_thickness"+int_str] - IS2_CS2_allseason_region_means["ice_thickness_cs2_ubris"], 'ko-', markersize=3, label='')
axes[1].set_ylabel('Thickness anomaly (m)')
axes[1].set_yticks([-0.5, 0, 0.5])
axes[1].set_ylim([-0.5, 0.9])
axes[1].legend(loc='upper right', ncol=2, frameon=False)
axes[1].grid(axis='y') # Only show grid lines for the x-axis
axes[1].tick_params(axis='x', which='both', bottom=False, top=False, labelbottom=False) # Remove x-axis ticks and labels
# Plot snow depth
axes[2].plot(IS2_CS2_allseason_region_means.time, IS2_CS2_allseason_region_means["snow_depth"+int_str], 'ko-', markersize=3, label='NSIM')
axes[2].plot(IS2_CS2_allseason_region_means.time, IS2_CS2_allseason_region_means["snow_depth_sm"+int_str], 'bo-', markersize=3, label='SM-LG')
axes[2].plot(IS2_CS2_allseason_region_means.time, IS2_CS2_allseason_region_means["cs2is2_snow_depth"], 'mo-', markersize=3, label='IS-2/CS-2')
axes[2].plot(IS2_CS2_allseason_region_means.time, IS2_CS2_allseason_region_means["snow_depth_mw99"+int_str], 'go-', markersize=3, label='MW99')
axes[2].set_ylabel('Snow depth (m)')
axes[2].set_yticks([0, 0.1, 0.2])
axes[2].set_ylim([0,0.3])
axes[2].legend(loc='upper right', ncol=4, frameon=False)
axes[2].grid(axis='y') # Only show grid lines for the y-axis
axes[2].tick_params(axis='x', which='both', bottom=False, top=False, labelbottom=False) # Remove x-axis ticks and labels
# Plot snow density
axes[3].plot(IS2_CS2_allseason_region_means.time, IS2_CS2_allseason_region_means["snow_density"+int_str], 'ko-', markersize=3, label='NSIM')
axes[3].plot(IS2_CS2_allseason_region_means.time, IS2_CS2_allseason_region_means["snow_density_sm"+int_str], 'bo-', markersize=3, label='SM-LG')
axes[3].plot(IS2_CS2_allseason_region_means.time, IS2_CS2_allseason_region_means["snow_density_w99"+int_str], 'go-', markersize=3, label='MW99')
axes[3].set_ylabel('Snow density (kg/m³)')
axes[3].set_yticks([180, 240, 300, 360])
axes[3].set_ylim([160,400])
axes[3].legend(loc='upper right', ncol=4,frameon=False)
axes[3].grid(axis='y') # Only show grid lines for the x-axis
axes[3].tick_params(axis='x', which='both', bottom=False, top=False, labelbottom=False) # Remove x-axis ticks and labels
# Plot ice density
#axes[4].hlines(y=916, xmin=np.datetime64('2018-11-15'), xmax=np.datetime64('2021-07-15'), color='b', linestyle='-', label='IS-2', linewidth=1)
axes[4].plot(IS2_CS2_allseason_region_means.time, IS2_CS2_allseason_region_means["ice_density"], 'bo-', label='IS-2', markersize=3)
axes[4].plot(IS2_CS2_allseason_region_means.time, IS2_CS2_allseason_region_means["cs2_sea_ice_density_UBRIS"], 'ro-', label='CS-2', markersize=3)
axes[4].plot(IS2_CS2_allseason_region_means.time, IS2_CS2_allseason_region_means["ice_density_j22"+int_str], 'yo-', label='J22', markersize=3)
axes[4].set_ylabel('Ice density (kg/m³)')
axes[4].grid(axis='y') # Only show grid lines for the x-axis
axes[4].set_ylim([870,970])
#axes[4].set_xlim([np.datetime64('2018-09'), np.datetime64('2021-09')])
axes[4].tick_params(axis='x', which='both', bottom=False, top=False, labelbottom=False) # Remove x-axis ticks and labels
axes[4].legend(loc='upper right', ncol=4, frameon=False)
# Plot ice concentration
axes[5].plot(IS2_CS2_allseason_region_means.time, IS2_CS2_allseason_region_means.sea_ice_conc, 'mo-', markersize=3)
axes[5].set_ylabel('Ice concentration')
axes[5].grid(axis='y') # Only show grid lines for the x-axis
axes[5].set_ylim([0.8,1.])
# Set x-axis limits and ticks
start_date = np.datetime64('2018-11-01')
end_date = np.datetime64('2021-07-30')
for ax in axes:
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.set_xlim([start_date, end_date]) # Apply same x-limits to all subplots
axes[5].set_xlabel('Date')
axes[5].xaxis.set_major_locator(mdates.MonthLocator(bymonthday=15,interval=3)) # Show ticks every 6 months
axes[5].xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m')) # Format as YYYY-MM
panel_letters = ['(a)', '(b)', '(c)', '(d)', '(e)', '(f)']
for ax, letter in zip(axes, panel_letters):
ax.text(0.01, 0.93, letter, transform=ax.transAxes, va='bottom', ha='left')
plt.subplots_adjust(left=0.08, right=0.99, top=0.99, bottom=0.06, hspace=0.11) # Adjust these values to reduce whitespace
plt.savefig('./figs/allseason_timeseries_6rows'+int_str+'r'+str(innerArctic[0])+'-'+str(innerArctic[1])+'.pdf', dpi=300)
plt.show()
No artists with labels found to put in legend. Note that artists whose label start with an underscore are ignored when legend() is called with no argument.
Monthly masking#
Now I want to mask the data by month as things are a bit messy with averaging summer months across the IS-2 thickness dataset as lost of missing data which impacts the averages a lot. So the approach here is stricter masking, using just grid-cells where both datasets (IS-2 and CS-2) have data every month. A perennial masking of sorts.
def mask_and_concat_ice_thickness_monthly(IS2_CS2_allseason, int_str='_int'):
"""
Creates a new dataset with just the two key variables and masks entire time series
for grid cells that have any NaN values in either variable at any time.
"""
vars_to_keep = ["ice_thickness_cs2_ubris", "ice_thickness_sm"+int_str]
# Create new dataset with just the two variables
dataset = IS2_CS2_allseason[vars_to_keep]
# Find indices where any time is null for either variable
var1_null = dataset[vars_to_keep[0]].isnull().any(dim='time').values
var2_null = dataset[vars_to_keep[1]].isnull().any(dim='time').values
# Combine masks - True where either variable has nulls
has_nulls = var1_null | var2_null
# Apply mask (keep only points where has_nulls is False)
masked_dataset = dataset.where(~has_nulls)
#print(f"Number of valid grid cells: {(~has_nulls).sum().values}")
return masked_dataset, var1_null, var2_null, has_nulls
# Assuming 'IS2_CS2_allseason' is your dataset and 'ice_thickness_sm_int' is the variable of interest
summer_months = [5, 6, 7, 8] # May=5, June=6, July=7, August=8
# Exclude July 2019
exclude_time = '2019-07'
IS2_CS2_summer = IS2_CS2_allseason.sel(
time=IS2_CS2_allseason.time.dt.month.isin(summer_months)
)
IS2_CS2_summer = IS2_CS2_summer.sel(
time=~((IS2_CS2_summer.time.dt.year == 2019) & (IS2_CS2_summer.time.dt.month == 7))
)
# Apply the masking
IS2_CS2_summer_month_mask, var1_null, var2_null, has_nulls = mask_and_concat_ice_thickness_monthly(IS2_CS2_summer)
# Generate the same monthly maps as before but with the stricter common/perennial mask applied.
grouped_by_year = IS2_CS2_summer_month_mask.groupby('time.year')
# Create a figure with subplots
fig, axes = plt.subplots(nrows=len(grouped_by_year.groups), ncols=len(summer_months),
subplot_kw={'projection': ccrs.NorthPolarStereo(central_longitude=-45)},
figsize=(7, 6.1))
# Iterate over each year
for year_idx, (year, year_data) in enumerate(grouped_by_year):
grouped_by_month = year_data.groupby('time.month')
for month_idx, (month, month_data) in enumerate(grouped_by_month):
print(year_idx, month_idx)
if (int(year_idx)==0) & (int(month_idx)==2):
print('shift axes')
ax = axes[year_idx, month_idx+1]
else:
ax = axes[year_idx, month_idx]
im = month_data.isel(time=0).ice_thickness_sm_int.plot(ax=ax, x="longitude", y="latitude",
transform=ccrs.PlateCarree(), cmap="viridis",
vmin=0, vmax=5, add_colorbar=False)
ax.set_extent([-179, 179, 66, 90], crs=ccrs.PlateCarree())
ax.text(0.98, 0.9, f'{year}-{month:02d}', transform=ax.transAxes, va='bottom', ha='right')
ax.set_title('')
ax.coastlines(linewidth=0.15, color = 'black', zorder = 2) # Coastlines
ax.add_feature(cfeature.LAND, color ='0.95', zorder = 1) # Land
#ax.add_feature(cfeature.LAKES, color = 'grey', zorder = 5) # Lakes
ax.gridlines(draw_labels=False, linewidth=0.25, color='gray', alpha=0.7, linestyle='--', zorder=3) # Gridlines
# After your main plotting loop, add an empty map at [0, 2]
ax = axes[0, 2]
ax.set_extent([-179, 179, 66, 90], crs=ccrs.PlateCarree())
ax.set_title('')
ax.coastlines(linewidth=0.15, color='black', zorder=2)
ax.add_feature(cfeature.LAND, color='0.95', zorder=1)
ax.gridlines(draw_labels=False, linewidth=0.25, color='gray', alpha=0.7, linestyle='--', zorder=3)
ax.text(0.98, 0.82, '2019-07\n(no data)', transform=ax.transAxes, va='bottom', ha='right')
#axes[0, 2].set_visible(False)
axes[2, 3].set_visible(False)
cbar = fig.colorbar(im, ax=axes, orientation='horizontal', extend='both', fraction=0.03, pad=0.1)
cbar.set_label('IS-2 sea ice thickness (m)')
plt.subplots_adjust(left=0.02, right=0.98, top=0.98, bottom=0.15, wspace=0.01, hspace=0.03)
plt.savefig('./figs/maps_is2_summer_months_monthlymask.png', dpi=300, facecolor="white")
plt.show()
0 0
0 1
0 2
shift axes
1 0
1 1
1 2
1 3
2 0
2 1
2 2
# Again group the data by year and plot the same maps as before but with the stricter mask applied.
yearly_means = IS2_CS2_summer_month_mask.ice_thickness_sm_int.groupby('time.year').mean(skipna=False)
yearly_means_CS2 = IS2_CS2_summer_month_mask.ice_thickness_cs2_ubris.groupby('time.year').mean(skipna=False)
yearly_means = yearly_means.assign_coords(source='IS2')
yearly_means_CS2 = yearly_means_CS2.assign_coords(source='CS2')
combined_data = xr.concat([yearly_means, yearly_means_CS2], dim='source')
#combined_data
# Plot IS-2 and CS-2 summer mean maps with the stricter mask applied.
panel_letters = ['(a) 2019 IS-2/SM-LG', '(b) 2020 IS-2/SM-LG', '(c) 2021 IS-2/SM-LG',
'(d) 2019 CS-2/SM-LG', '(e) 2020 CS-2/SM-LG', '(f) 2021 CS-2/SM-LG']
fig_width=6.8
fig_height=5.4
im1 = combined_data.plot(x="longitude", y="latitude", col='year', row='source', transform=ccrs.PlateCarree(),
cmap="viridis", zorder=3, add_labels=False,
subplot_kws={'projection':ccrs.NorthPolarStereo(central_longitude=-45)},
cbar_kwargs={'pad':0.01,'shrink': 0.3,'extend':'both',
'label':'Sea ice thickness (m)', 'location':'bottom'},
vmin=0, vmax=5,
figsize=(fig_width, fig_height))
i=0
for i, ax in enumerate(im1.axes.flatten()):
ax.coastlines(linewidth=0.15, color='black', zorder=2)
ax.add_feature(cfeature.LAND, color='0.95', zorder=1)
ax.gridlines(draw_labels=False, linewidth=0.25, color='gray', alpha=0.7, linestyle='--', zorder=3)
ax.set_extent([-179, 179, 65, 90], crs=ccrs.PlateCarree())
ax.text(0.02, 0.93, panel_letters[i], transform=ax.transAxes, fontsize=7, va='bottom', ha='left')
ax.set_title('', fontsize=8)
i+=1
plt.savefig('./figs/maps_summer_thickness_comps_monthmask.png', dpi=300, facecolor="white")
plt.show()
plt.close()
# plot the annual anomalies for each variable using the strict monthly masking
combined_data_anoms = xr.concat([yearly_means-yearly_means.mean(dim='year'), yearly_means_CS2-yearly_means_CS2.mean(dim='year')], dim='source')
panel_letters = ['(a) 2019 IS-2/SM-LG', '(b) 2020 IS-2/SM-LG', '(c) 2021 IS-2/SM-LG',
'(d) 2019 CS-2/SM-LG', '(e) 2020 CS-2/SM-LG', '(f) 2021 CS-2/SM-LG']
fig_width=6.8
fig_height=5.4
im1 = combined_data_anoms.plot(x="longitude", y="latitude", col='year', row='source', transform=ccrs.PlateCarree(),
cmap="RdBu", zorder=3, add_labels=False,
subplot_kws={'projection':ccrs.NorthPolarStereo(central_longitude=-45)},
cbar_kwargs={'pad':0.03,'shrink': 0.3,'extend':'both', 'ticks': np.arange(-1.5, 1.6, 0.5),
'label':'thickness anomaly (m)', 'location':'bottom'},
vmin=-1.25, vmax=1.25,
figsize=(fig_width, fig_height))
# Add map features to all subplots
i=0
for i, ax in enumerate(im1.axes.flatten()):
ax.coastlines(linewidth=0.15, color='black', zorder=2)
ax.add_feature(cfeature.LAND, color='0.95', zorder=1)
ax.gridlines(draw_labels=False, linewidth=0.25, color='gray', alpha=0.7, linestyle='--', zorder=3)
ax.set_extent([-179, 179, 65, 90], crs=ccrs.PlateCarree())
ax.text(0.02, 0.93, panel_letters[i], transform=ax.transAxes, fontsize=7, va='bottom', ha='left')
ax.set_title('', fontsize=8)
i+=1
plt.gcf().subplots_adjust(bottom=0.17, wspace=0.03)
plt.savefig('./figs/maps_summer_thickness_anoms_monthmask.png', dpi=300, facecolor="white")
plt.show()
plt.close()
# Plot IS-2 and CS-2 summer anomaly mean maps using the stricter mask applied.
anoms = yearly_means - yearly_means_CS2
fig_width=6.8
fig_height=3.4
panel_letters = ['(g) 2019 IS-2 - CS-2', '(h) 2020 IS-2 - CS-2', '(i) 2021 IS-2 - CS-2']
im = anoms.plot(x="longitude", y="latitude", col_wrap=3, col='year', transform=ccrs.PlateCarree(), cmap="RdBu", zorder=3,
cbar_kwargs={'pad':0.01,'shrink': 0.3,'extend':'both', 'label':'thickness difference (m)', 'location':'bottom'},
vmin=-2.5, vmax=2.5,
subplot_kws={'projection':ccrs.NorthPolarStereo(central_longitude=-45)},
figsize=(fig_width, fig_height))
ax_iter = im.axes
i=0
for ax in ax_iter.flatten():
ax.coastlines(linewidth=0.15, color = 'black', zorder = 2) # Coastlines
ax.add_feature(cfeature.LAND, color ='0.95', zorder = 1) # Land
#ax.add_feature(cfeature.LAKES, color = 'grey', zorder = 5) # Lakes
ax.gridlines(draw_labels=False, linewidth=0.25, color='gray', alpha=0.7, linestyle='--', zorder=3) # Gridlines
ax.set_extent([-179, 179, 65, 90], crs=ccrs.PlateCarree()) # Set extent to zoom in on Arctic
ax.set_title('', fontsize=8, horizontalalignment="right",verticalalignment="top", x=0.97, y=0.95, zorder=4)
ax.text(0.02, 0.93, panel_letters[i], fontsize=7, transform=ax.transAxes, va='bottom', ha='left')
i+=1
plt.subplots_adjust(left=0.07, right=0.96, wspace=0.07, bottom=0.2)
plt.savefig('./figs/maps_summer_thickness_diff_monthmask.png', dpi=300, facecolor="white", bbox_inches='tight')
Let’s do the same for winter (Jan through April) to see if we get similar results.#
# Assuming 'IS2_CS2_allseason' is your dataset and 'ice_thickness_sm_int' is the variable of interest
winter_months = [1, 2, 3, 4] # May=5, June=6, July=7, August=8
# Select only the summer months
IS2_CS2_winter = IS2_CS2_allseason.sel(time=IS2_CS2_allseason.time.dt.month.isin(winter_months))
# Apply the masking
IS2_CS2_winter_month_mask, var1_null, var2_null, has_nulls = mask_and_concat_ice_thickness_monthly(IS2_CS2_winter)
# Generate the same monthly maps as before but with the stricter mask applied.
grouped_by_year = IS2_CS2_winter_month_mask.groupby('time.year')
# Create a figure with subplots
fig, axes = plt.subplots(nrows=len(grouped_by_year.groups), ncols=len(winter_months),
subplot_kw={'projection': ccrs.NorthPolarStereo(central_longitude=-45)},
figsize=(7, 6.1))
# Iterate over each year
for year_idx, (year, year_data) in enumerate(grouped_by_year):
#print(year_idx)# Group by month within each year
grouped_by_month = year_data.groupby('time.month')
for month_idx, (month, month_data) in enumerate(grouped_by_month):
ax = axes[year_idx, month_idx]
im = month_data.isel(time=0).ice_thickness_sm_int.plot(ax=ax, x="longitude", y="latitude",
transform=ccrs.PlateCarree(), cmap="viridis",
vmin=0, vmax=5, add_colorbar=False)
ax.set_extent([-179, 179, 66, 90], crs=ccrs.PlateCarree())
ax.text(0.98, 0.9, f'{year}-{month:02d}', transform=ax.transAxes, va='bottom', ha='right')
ax.set_title('')
ax.coastlines(linewidth=0.15, color = 'black', zorder = 2) # Coastlines
ax.add_feature(cfeature.LAND, color ='0.95', zorder = 1) # Land
#ax.add_feature(cfeature.LAKES, color = 'grey', zorder = 5) # Lakes
ax.gridlines(draw_labels=False, linewidth=0.25, color='gray', alpha=0.7, linestyle='--', zorder=3) # Gridlines
axes[2, 3].set_visible(False)
cbar = fig.colorbar(im, ax=axes, orientation='horizontal', extend='both', fraction=0.03, pad=0.1)
cbar.set_label('IS-2 sea ice thickness (m)')
plt.subplots_adjust(left=0.02, right=0.98, top=0.98, bottom=0.15, wspace=0.01, hspace=0.03)
plt.savefig('./figs/maps_is2_winter_months_monthlymask.png', dpi=300, facecolor="white")
plt.show()
# Again group the data by year and plot the same maps as before but with the stricter mask applied.
yearly_means = IS2_CS2_winter_month_mask.ice_thickness_sm_int.groupby('time.year').mean(skipna=False)
yearly_means_CS2 = IS2_CS2_winter_month_mask.ice_thickness_cs2_ubris.groupby('time.year').mean(skipna=False)
yearly_means = yearly_means.assign_coords(source='IS2')
yearly_means_CS2 = yearly_means_CS2.assign_coords(source='CS2')
combined_data = xr.concat([yearly_means, yearly_means_CS2], dim='source')
#combined_data
# Plot IS-2 and CS-2 summer mean maps with the stricter mask applied.
panel_letters = ['(a) 2019 IS-2/SM-LG', '(b) 2020 IS-2/SM-LG', '(c) 2021 IS-2/SM-LG',
'(d) 2019 CS-2/SM-LG', '(e) 2020 CS-2/SM-LG', '(f) 2021 CS-2/SM-LG']
fig_width=6.8
fig_height=5.4
im1 = combined_data.plot(x="longitude", y="latitude", col='year', row='source', transform=ccrs.PlateCarree(),
cmap="viridis", zorder=3, add_labels=False,
subplot_kws={'projection':ccrs.NorthPolarStereo(central_longitude=-45)},
cbar_kwargs={'pad':0.01,'shrink': 0.3,'extend':'both',
'label':'Sea ice thickness (m)', 'location':'bottom'},
vmin=0, vmax=5,
figsize=(fig_width, fig_height))
i=0
for i, ax in enumerate(im1.axes.flatten()):
ax.coastlines(linewidth=0.15, color='black', zorder=2)
ax.add_feature(cfeature.LAND, color='0.95', zorder=1)
ax.gridlines(draw_labels=False, linewidth=0.25, color='gray', alpha=0.7, linestyle='--', zorder=3)
ax.set_extent([-179, 179, 65, 90], crs=ccrs.PlateCarree())
ax.text(0.02, 0.93, panel_letters[i], transform=ax.transAxes, fontsize=7, va='bottom', ha='left')
ax.set_title('', fontsize=8)
i+=1
plt.savefig('./figs/maps_winter_thickness_comps_monthmask.png', dpi=300, facecolor="white")
plt.show()
plt.close()
# plot the annual anomalies for each variable using the strict monthly masking
combined_data_anoms = xr.concat([yearly_means-yearly_means.mean(dim='year'), yearly_means_CS2-yearly_means_CS2.mean(dim='year')], dim='source')
panel_letters = ['(a) 2019 IS-2/SM-LG', '(b) 2020 IS-2/SM-LG', '(c) 2021 IS-2/SM-LG',
'(d) 2019 CS-2/SM-LG', '(e) 2020 CS-2/SM-LG', '(f) 2021 CS-2/SM-LG']
fig_width=6.8
fig_height=5.4
im1 = combined_data_anoms.plot(x="longitude", y="latitude", col='year', row='source', transform=ccrs.PlateCarree(),
cmap="RdBu", zorder=3, add_labels=False,
subplot_kws={'projection':ccrs.NorthPolarStereo(central_longitude=-45)},
cbar_kwargs={'pad':0.03,'shrink': 0.3,'extend':'both', 'ticks': np.arange(-1.5, 1.6, 0.5),
'label':'thickness anomaly (m)', 'location':'bottom'},
vmin=-1.25, vmax=1.25,
figsize=(fig_width, fig_height))
# Add map features to all subplots
i=0
for i, ax in enumerate(im1.axes.flatten()):
ax.coastlines(linewidth=0.15, color='black', zorder=2)
ax.add_feature(cfeature.LAND, color='0.95', zorder=1)
ax.gridlines(draw_labels=False, linewidth=0.25, color='gray', alpha=0.7, linestyle='--', zorder=3)
ax.set_extent([-179, 179, 65, 90], crs=ccrs.PlateCarree())
ax.text(0.02, 0.93, panel_letters[i], transform=ax.transAxes, fontsize=7, va='bottom', ha='left')
ax.set_title('', fontsize=8)
i+=1
plt.gcf().subplots_adjust(bottom=0.17, wspace=0.03)
plt.savefig('./figs/maps_winter_thickness_anoms_monthmask.png', dpi=300, facecolor="white")
plt.show()
plt.close()
# Plot IS-2 and CS-2 summer anomaly mean maps using the stricter mask applied.
anoms = yearly_means - yearly_means_CS2
fig_width=6.8
fig_height=3.4
panel_letters = ['(g) 2019 IS-2 - CS-2', '(h) 2020 IS-2 - CS-2', '(i) 2021 IS-2 - CS-2']
im = anoms.plot(x="longitude", y="latitude", col_wrap=3, col='year', transform=ccrs.PlateCarree(), cmap="RdBu", zorder=3,
cbar_kwargs={'pad':0.01,'shrink': 0.3,'extend':'both', 'label':'thickness difference (m)', 'location':'bottom'},
vmin=-2.5, vmax=2.5,
subplot_kws={'projection':ccrs.NorthPolarStereo(central_longitude=-45)},
figsize=(fig_width, fig_height))
ax_iter = im.axes
i=0
for ax in ax_iter.flatten():
ax.coastlines(linewidth=0.15, color = 'black', zorder = 2) # Coastlines
ax.add_feature(cfeature.LAND, color ='0.95', zorder = 1) # Land
#ax.add_feature(cfeature.LAKES, color = 'grey', zorder = 5) # Lakes
ax.gridlines(draw_labels=False, linewidth=0.25, color='gray', alpha=0.7, linestyle='--', zorder=3) # Gridlines
ax.set_extent([-179, 179, 65, 90], crs=ccrs.PlateCarree()) # Set extent to zoom in on Arctic
ax.set_title('', fontsize=8, horizontalalignment="right",verticalalignment="top", x=0.97, y=0.95, zorder=4)
ax.text(0.02, 0.93, panel_letters[i], fontsize=7, transform=ax.transAxes, va='bottom', ha='left')
i+=1
plt.subplots_adjust(left=0.07, right=0.96, wspace=0.07, bottom=0.2)
plt.savefig('./figs/maps_winter_thickness_diff_monthmask.png', dpi=300, facecolor="white", bbox_inches='tight')
