ICESat-2–SMOS–SMAP fusion: intro and V4 comparison#
Summary: Introduces the daily IS2SMGPSIT-V1 (ICESat-2–SMOS–SMAP fused) product and compares pan-Arctic thickness, freeboard, and snow depth with IS2SITMOGR4-V4 (monthly, ICESat-2 only).
Author: Alek Petty
Version history: Version 1 (08/2026)
import xarray as xr
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
import os
import warnings
warnings.filterwarnings('ignore')
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from utils.read_data_utils import read_IS2SITMOGR4, read_is2smgpsitv1_zarr
import matplotlib as mpl
%config InlineBackend.figure_format = 'retina'
mpl.rcParams['figure.dpi'] = 300
mpl.rcParams.update({
'text.usetex': False,
'font.family': 'sans-serif',
'mathtext.fontset': 'stixsans',
'font.size': 8,
'axes.labelsize': 8,
'xtick.labelsize': 8,
'ytick.labelsize': 8,
'legend.fontsize': 8,
})
mpl.rcParams['font.sans-serif'] = ['Arial']
Load datasets#
Load IS2SITMOGR4-V4 (monthly) and IS2SMGPSIT-V1 (daily) from S3 Zarr.
# Monthly V4 product (ICESat-2 only)
IS2_v4 = read_IS2SITMOGR4(data_type='zarr-s3-v4', persist=False)
# Shift V4 monthly timestamps to the 15th of each month for consistency
_v4_time = pd.to_datetime(IS2_v4['time'].values)
_v4_start = _v4_time.to_period('M').to_timestamp(how='start')
_v4_mid = _v4_start + pd.Timedelta(days=14)
IS2_v4 = IS2_v4.assign_coords(time=_v4_mid.values)
print('IS2SITMOGR4-V4 (monthly, time shifted to 15th)')
print(IS2_v4)
load zarr from S3 bucket
zarr_path: s3://icesat-2-sea-ice-us-west-2/IS2SITMOGR4_V4/zarr/IS2SITMOGR4_V4_201811-202504.zarr
IS2SITMOGR4-V4 (monthly, time shifted to 15th)
<xarray.Dataset> Size: 1GB
Dimensions: (time: 54, y: 448, x: 304)
Coordinates:
latitude (y, x) float32 545kB 31.1 31.2 ... 34.47
longitude (y, x) float32 545kB 168.3 168.1 ... -9.999
* x (x) float32 1kB -3.838e+06 ... 3.738e+06
* y (y) float32 2kB 5.838e+06 ... -5.338e+06
* time (time) datetime64[ns] 432B 2018-11-15 ......
Data variables: (12/38)
crs (time) int32 216B dask.array<chunksize=(54,), meta=np.ndarray>
freeboard (time, y, x) float32 29MB dask.array<chunksize=(1, 448, 304), meta=np.ndarray>
freeboard_int (time, y, x) float32 29MB dask.array<chunksize=(1, 448, 304), meta=np.ndarray>
grid_cell_area (y, x) float32 545kB dask.array<chunksize=(448, 304), meta=np.ndarray>
ice_density (time, y, x) float32 29MB dask.array<chunksize=(1, 448, 304), meta=np.ndarray>
ice_density_j22 (time, y, x) float32 29MB dask.array<chunksize=(1, 448, 304), meta=np.ndarray>
... ...
snow_depth_mw99 (time, y, x) float32 29MB dask.array<chunksize=(1, 448, 304), meta=np.ndarray>
snow_depth_mw99_int (time, y, x) float32 29MB dask.array<chunksize=(1, 448, 304), meta=np.ndarray>
snow_depth_sm_e5 (time, y, x) float32 29MB dask.array<chunksize=(1, 448, 304), meta=np.ndarray>
snow_depth_sm_e5_int (time, y, x) float32 29MB dask.array<chunksize=(1, 448, 304), meta=np.ndarray>
snow_depth_sm_m2 (time, y, x) float32 29MB dask.array<chunksize=(1, 448, 304), meta=np.ndarray>
snow_depth_sm_m2_int (time, y, x) float32 29MB dask.array<chunksize=(1, 448, 304), meta=np.ndarray>
Attributes:
contact: Alek Petty (akpetty@umd.edu)
description: IS2SITMOGR4, version 004. Gridded Nov 2018 winter Arctic se...
history: Created 20/10/25
reference: Official NSIDC data doi: 10.5067/CV6JEXEE31HF. Derived data...
# Read IS2SMGPSIT-V1 from S3 and rebuild the local cache (producer notebook)
import time
t0 = time.perf_counter()
IS2_SMOS_SMAP = read_is2smgpsitv1_zarr(
persist=True,
cache=True,
cache_refresh=False,
)
IS2_SMOS_SMAP.info()
print(f"Elapsed: {time.perf_counter() - t0:.1f} s")
Loading IS2-SMOS-SMAP (is2smsitgp) Zarr from local cache
cache_path: ./data/cache/GPSat_multivar_20181101-20250430.zarr
xarray.Dataset {
dimensions:
time = 1635 ;
y = 448 ;
x = 304 ;
variables:
int32 crs(time) ;
crs:GeoTransform = -3850000.0 25000.0 0 5850000.0 0 -25000.0 ;
crs:false_easting = 0.0 ;
crs:false_northing = 0.0 ;
crs:grid_mapping_name = polar_stereographic ;
crs:latitude_of_projection_origin = 90.0 ;
crs:long_name = NSIDC Sea Ice Polar Stereographic North ;
crs:proj4text = +proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +k=1 +x_0=0 +y_0=0 +a=6378273 +b=6356889.449 +units=m +no_defs ;
crs:srid = urn:ogc:def:crs:EPSG::3411 ;
crs:standard_parallel = 70.0 ;
crs:straight_vertical_longitude_from_pole = -45.0 ;
crs:units = meters ;
float32 freeboard(time, y, x) ;
freeboard:description = Mean freeboard across the grid cell (m). Uses ICESat-2 total freeboard and CDR sea ice concentration (open water) in the solution; gridded via GPSat interpolation. ;
freeboard:grid_mapping = crs ;
freeboard:long_name = sea ice freeboard (ICESat-2 total freeboard) ;
freeboard:units = meters ;
float32 grid_cell_area(y, x) ;
grid_cell_area:description = Grid cell area from NSIDC Polar Stereographic North (EPSG:3411). ;
grid_cell_area:grid_mapping = crs ;
grid_cell_area:long_name = grid cell area ;
grid_cell_area:source = https://doi.org/10.5067/N6INPBT8Y104 ;
grid_cell_area:units = m2 ;
float32 ice_thickness(time, y, x) ;
ice_thickness:description = Mean sea ice thickness across the grid cell (m). From fusion of ICESat-2 and SMAP (thin ice) via GP interpolation (petty). ;
ice_thickness:grid_mapping = crs ;
ice_thickness:long_name = sea ice thickness ;
ice_thickness:units = meters ;
float32 ice_thickness_unc(time, y, x) ;
ice_thickness_unc:description = Standard deviation of thickness (sqrt of predictive variance from GP). Thickness from ICESat-2/SMAP fusion. ;
ice_thickness_unc:grid_mapping = crs ;
ice_thickness_unc:long_name = sea ice thickness uncertainty (standard deviation) ;
ice_thickness_unc:units = meters ;
float32 latitude(y, x) ;
latitude:grid_mapping = crs ;
latitude:long_name = latitude ;
latitude:units = degrees_north ;
float32 longitude(y, x) ;
longitude:grid_mapping = crs ;
longitude:long_name = longitude ;
longitude:units = degrees_east ;
float32 mean_freeboard_regions_1_5(time) ;
mean_freeboard_regions_1_5:description = Area-weighted mean freeboard over grid cells in regions 1–5. Freeboard uses ICESat-2 and CDR SIC with GPSat interpolation. ;
mean_freeboard_regions_1_5:grid_mapping = crs ;
mean_freeboard_regions_1_5:long_name = mean total freeboard (area-weighted, NSIDC regions 1–5 only) ;
mean_freeboard_regions_1_5:units = m ;
float32 mean_sea_ice_conc_regions_1_5(time) ;
mean_sea_ice_conc_regions_1_5:description = Mean sea ice concentration over grid cells in regions 1–5. ;
mean_sea_ice_conc_regions_1_5:grid_mapping = crs ;
mean_sea_ice_conc_regions_1_5:long_name = mean sea ice concentration (NSIDC regions 1–5 only) ;
mean_sea_ice_conc_regions_1_5:units = 1 ;
float32 mean_sea_ice_thickness_regions_1_5(time) ;
mean_sea_ice_thickness_regions_1_5:description = Area-weighted mean ice thickness over grid cells in regions 1–5: sum(ice_thickness × grid_cell_area) / sum(grid_cell_area). Thickness from fusion of ICESat-2 and SMAP (thin ice). ;
mean_sea_ice_thickness_regions_1_5:grid_mapping = crs ;
mean_sea_ice_thickness_regions_1_5:long_name = mean sea ice thickness (area-weighted, NSIDC regions 1–5 only) ;
mean_sea_ice_thickness_regions_1_5:units = m ;
float32 mean_snow_depth_regions_1_5(time) ;
mean_snow_depth_regions_1_5:description = Area-weighted mean snow depth over grid cells in regions 1–5. Snow depth uses CDR SIC in the solution; from NESOSIM v1.1. ;
mean_snow_depth_regions_1_5:grid_mapping = crs ;
mean_snow_depth_regions_1_5:long_name = mean snow depth (area-weighted, NSIDC regions 1–5 only) ;
mean_snow_depth_regions_1_5:units = m ;
float32 region_mask(y, x) ;
region_mask:description = NSIDC-0780 sea ice region mask (sea_ice_region_surface_mask). CAA (Canadian Archipelago, value 12) is retained in gridded fields; see global attribute caa_uncertainty_caution. ;
region_mask:flag_values = 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 30, 31, 32 ;
region_mask:grid_mapping = crs ;
region_mask:long_name = Northern Hemisphere region mask ;
region_mask:source = https://doi.org/10.5067/CYW3O8ZUNIWC ;
float32 snow_depth(time, y, x) ;
snow_depth:description = Mean snow depth across the grid cell (m). Uses CDR sea ice concentration in the solution; from NESOSIM v1.1 (IS-2 sub-sampled + GPSat interpolation). ;
snow_depth:grid_mapping = crs ;
snow_depth:long_name = snow depth on sea ice ;
snow_depth:units = meters ;
datetime64[ns] time(time) ;
time:long_name = time ;
time:standard_name = time ;
float32 total_sea_ice_volume(time) ;
total_sea_ice_volume:description = Sum of sea_ice_volume over all grid cells (km³). ;
total_sea_ice_volume:grid_mapping = crs ;
total_sea_ice_volume:long_name = total sea ice volume ;
total_sea_ice_volume:units = km3 ;
float32 total_sea_ice_volume_regions_1_5(time) ;
total_sea_ice_volume_regions_1_5:description = Sum of sea_ice_volume over grid cells in NSIDC-0780 regions 1–5 (Central Arctic through Laptev Sea), in km³. ;
total_sea_ice_volume_regions_1_5:grid_mapping = crs ;
total_sea_ice_volume_regions_1_5:long_name = total sea ice volume (NSIDC regions 1–5 only) ;
total_sea_ice_volume_regions_1_5:units = km3 ;
float32 x(x) ;
x:long_name = projection x coordinate ;
x:units = meters ;
float32 y(y) ;
y:long_name = projection y coordinate ;
y:units = meters ;
// global attributes:
:caa_uncertainty_caution = Canadian Arctic Archipelago (CAA; NSIDC-0780 sea_ice_region_surface_mask value 12) is included in gridded fields but has elevated uncertainty relative to the open Arctic Ocean: complex coastlines, land-fast ice, and limited altimetry/auxiliary coverage. Prefer open-ocean regions for quantitative applications, and treat CAA values with extra caution. ;
:contact = Alek Petty (akpetty@umd.edu) ;
:description = ICESat-2–SMOS–SMAP fused daily gridded product (GPSat multivar). ;
:history = Created 2026-07-29 ;
:inner_arctic_regions_note = Inner-Arctic scalars (mean_*_regions_1_5, total_sea_ice_volume_regions_1_5) use NSIDC-0780 regions 1–5 only (Central Arctic through Laptev Sea). Thickness/freeboard/snow means are area-weighted. total_sea_ice_volume is the pan-Arctic sum over all valid grid cells. ;
:n_timesteps = 1635 ;
:source = Combined from monthly IS2_interp_test_petty files (thickness/freeboard/snow_depth); grid area NSIDC0771; SIC CDR. ;
:time_coverage_end = 2025-04-30 ;
:time_coverage_start = 2018-11-01 ;
:training_season_caution = GP training is limited to a seasonal window (e.g. 1 Sep–30 Apr); there is no training data before 1 September or after 30 April in that setup. September and April sit at the edges of that window, so they have more limited in-window context than mid-season months—treat September and April fields with extra caution. ;
}Elapsed: 1.2 s
from utils.read_data_utils import is2smgpsit_domain_masks, area_weighted_spatial_mean, static_grid_cell_area, is2smgpsit_monthly_midmonth
# Pan-Arctic totals exclude CAA (region 12); IAO = NSIDC regions 1–5.
_masks_fused = is2smgpsit_domain_masks(IS2_SMOS_SMAP)
mask_pan_arctic_fused = _masks_fused['pan_arctic']
mask_iao_fused = _masks_fused['iao']
_masks_v4 = is2smgpsit_domain_masks(IS2_v4)
mask_pan_arctic_v4 = _masks_v4['pan_arctic']
mask_iao_v4 = _masks_v4['iao']
innerArctic = [1, 2, 3, 4, 5]
print('Pan-Arctic mask (excl. CAA) fused cells:', int(mask_pan_arctic_fused.sum().values))
print('IAO mask (regions 1–5) fused cells:', int(mask_iao_fused.sum().values))
_area_fused = static_grid_cell_area(IS2_SMOS_SMAP['grid_cell_area'])
_area_v4 = static_grid_cell_area(IS2_v4['grid_cell_area'])
Pan-Arctic mask (excl. CAA) fused cells: 37244
IAO mask (regions 1–5) fused cells: 11144
IS2_SMOS_SMAP.keys()
KeysView(<xarray.Dataset> Size: 4GB
Dimensions: (time: 1635, y: 448, x: 304)
Coordinates:
latitude (y, x) float32 545kB dask.array<chunksize=(448, 304), meta=np.ndarray>
longitude (y, x) float32 545kB dask.array<chunksize=(448, 304), meta=np.ndarray>
* time (time) datetime64[ns] 13kB 2018-11-01...
* x (x) float32 1kB -3.838e+06 ... 3.738e+06
* y (y) float32 2kB 5.838e+06 ... -5.338e+06
Data variables: (12/13)
crs (time) int32 7kB dask.array<chunksize=(100,), meta=np.ndarray>
freeboard (time, y, x) float32 891MB dask.array<chunksize=(100, 448, 304), meta=np.ndarray>
grid_cell_area (y, x) float32 545kB dask.array<chunksize=(448, 304), meta=np.ndarray>
ice_thickness (time, y, x) float32 891MB dask.array<chunksize=(100, 448, 304), meta=np.ndarray>
ice_thickness_unc (time, y, x) float32 891MB dask.array<chunksize=(100, 448, 304), meta=np.ndarray>
mean_freeboard_regions_1_5 (time) float32 7kB dask.array<chunksize=(100,), meta=np.ndarray>
... ...
mean_sea_ice_thickness_regions_1_5 (time) float32 7kB dask.array<chunksize=(100,), meta=np.ndarray>
mean_snow_depth_regions_1_5 (time) float32 7kB dask.array<chunksize=(100,), meta=np.ndarray>
region_mask (y, x) float32 545kB dask.array<chunksize=(448, 304), meta=np.ndarray>
snow_depth (time, y, x) float32 891MB dask.array<chunksize=(100, 448, 304), meta=np.ndarray>
total_sea_ice_volume (time) float32 7kB dask.array<chunksize=(100,), meta=np.ndarray>
total_sea_ice_volume_regions_1_5 (time) float32 7kB dask.array<chunksize=(100,), meta=np.ndarray>
Attributes:
caa_uncertainty_caution: Canadian Arctic Archipelago (CAA; NSIDC-0780 ...
contact: Alek Petty (akpetty@umd.edu)
description: ICESat-2–SMOS–SMAP fused daily gridded produc...
history: Created 2026-07-29
inner_arctic_regions_note: Inner-Arctic scalars (mean_*_regions_1_5, tot...
n_timesteps: 1635
source: Combined from monthly IS2_interp_test_petty f...
time_coverage_end: 2025-04-30
time_coverage_start: 2018-11-01
training_season_caution: GP training is limited to a seasonal window (...)
IS2_SMOS_SMAP.time
<xarray.DataArray 'time' (time: 1635)> Size: 13kB
array(['2018-11-01T00:00:00.000000000', '2018-11-02T00:00:00.000000000',
'2018-11-03T00:00:00.000000000', ..., '2025-04-28T00:00:00.000000000',
'2025-04-29T00:00:00.000000000', '2025-04-30T00:00:00.000000000'],
dtype='datetime64[ns]')
Coordinates:
* time (time) datetime64[ns] 13kB 2018-11-01 2018-11-02 ... 2025-04-30
Attributes:
long_name: time
standard_name: timeCompare structure and time coverage#
V4 is monthly; the fused product is daily. Both use a similar polar grid (y, x).
print('V4 time range:', IS2_v4.time.values.min(), 'to', IS2_v4.time.values.max())
print('V4 dimensions:', dict(IS2_v4.dims))
print('V4 data variables:', list(IS2_v4.data_vars))
print()
print('Fused product time range:', IS2_SMOS_SMAP.time.values.min(), 'to', IS2_SMOS_SMAP.time.values.max())
print('Fused dimensions:', dict(IS2_SMOS_SMAP.dims))
print('Fused data variables:', list(IS2_SMOS_SMAP.data_vars))
V4 time range: 2018-11-15T00:00:00.000000000 to 2025-04-15T00:00:00.000000000
V4 dimensions: {'time': 54, 'y': 448, 'x': 304}
V4 data variables: ['crs', 'freeboard', 'freeboard_int', 'grid_cell_area', 'ice_density', 'ice_density_j22', 'ice_thickness', 'ice_thickness_int', 'ice_thickness_j22', 'ice_thickness_j22_int', 'ice_thickness_mw99', 'ice_thickness_mw99_int', 'ice_thickness_sm_e5', 'ice_thickness_sm_e5_int', 'ice_thickness_sm_m2', 'ice_thickness_sm_m2_int', 'ice_thickness_unc', 'ice_thickness_unc_freeboard', 'ice_thickness_unc_ice_density', 'ice_thickness_unc_snow_density', 'ice_thickness_unc_snow_depth', 'ice_type', 'mean_day_of_month', 'num_segments', 'region_mask', 'sea_ice_conc', 'snow_density', 'snow_density_sm_e5', 'snow_density_sm_m2', 'snow_density_w99', 'snow_depth', 'snow_depth_int', 'snow_depth_mw99', 'snow_depth_mw99_int', 'snow_depth_sm_e5', 'snow_depth_sm_e5_int', 'snow_depth_sm_m2', 'snow_depth_sm_m2_int']
Fused product time range: 2018-11-01T00:00:00.000000000 to 2025-04-30T00:00:00.000000000
Fused dimensions: {'time': 1635, 'y': 448, 'x': 304}
Fused data variables: ['crs', 'freeboard', 'grid_cell_area', 'ice_thickness', 'ice_thickness_unc', 'mean_freeboard_regions_1_5', 'mean_sea_ice_conc_regions_1_5', 'mean_sea_ice_thickness_regions_1_5', 'mean_snow_depth_regions_1_5', 'region_mask', 'snow_depth', 'total_sea_ice_volume', 'total_sea_ice_volume_regions_1_5']
IS2_SMOS_SMAP
<xarray.Dataset> Size: 4GB
Dimensions: (time: 1635, y: 448, x: 304)
Coordinates:
latitude (y, x) float32 545kB dask.array<chunksize=(448, 304), meta=np.ndarray>
longitude (y, x) float32 545kB dask.array<chunksize=(448, 304), meta=np.ndarray>
* time (time) datetime64[ns] 13kB 2018-11-01...
* x (x) float32 1kB -3.838e+06 ... 3.738e+06
* y (y) float32 2kB 5.838e+06 ... -5.338e+06
Data variables: (12/13)
crs (time) int32 7kB dask.array<chunksize=(100,), meta=np.ndarray>
freeboard (time, y, x) float32 891MB dask.array<chunksize=(100, 448, 304), meta=np.ndarray>
grid_cell_area (y, x) float32 545kB dask.array<chunksize=(448, 304), meta=np.ndarray>
ice_thickness (time, y, x) float32 891MB dask.array<chunksize=(100, 448, 304), meta=np.ndarray>
ice_thickness_unc (time, y, x) float32 891MB dask.array<chunksize=(100, 448, 304), meta=np.ndarray>
mean_freeboard_regions_1_5 (time) float32 7kB dask.array<chunksize=(100,), meta=np.ndarray>
... ...
mean_sea_ice_thickness_regions_1_5 (time) float32 7kB dask.array<chunksize=(100,), meta=np.ndarray>
mean_snow_depth_regions_1_5 (time) float32 7kB dask.array<chunksize=(100,), meta=np.ndarray>
region_mask (y, x) float32 545kB dask.array<chunksize=(448, 304), meta=np.ndarray>
snow_depth (time, y, x) float32 891MB dask.array<chunksize=(100, 448, 304), meta=np.ndarray>
total_sea_ice_volume (time) float32 7kB dask.array<chunksize=(100,), meta=np.ndarray>
total_sea_ice_volume_regions_1_5 (time) float32 7kB dask.array<chunksize=(100,), meta=np.ndarray>
Attributes:
caa_uncertainty_caution: Canadian Arctic Archipelago (CAA; NSIDC-0780 ...
contact: Alek Petty (akpetty@umd.edu)
description: ICESat-2–SMOS–SMAP fused daily gridded produc...
history: Created 2026-07-29
inner_arctic_regions_note: Inner-Arctic scalars (mean_*_regions_1_5, tot...
n_timesteps: 1635
source: Combined from monthly IS2_interp_test_petty f...
time_coverage_end: 2025-04-30
time_coverage_start: 2018-11-01
training_season_caution: GP training is limited to a seasonal window (...Basic thickness comparison#
Inner-Arctic-Ocean (NSIDC regions 1–5) area-weighted mean ice_thickness_int (V4) vs the fused mean_sea_ice_thickness_regions_1_5 scalar, with May–August masked.
# Calendar-month means from all daily fields in each represented month.
# Confirm that every represented month contains its full complement of daily outputs.
_daily_index = pd.DatetimeIndex(IS2_SMOS_SMAP.time.values)
_month_counts = pd.Series(1, index=_daily_index).groupby(_daily_index.to_period('M')).sum()
_expected_counts = pd.Series(_month_counts.index.days_in_month, index=_month_counts.index)
assert _month_counts.equals(_expected_counts), 'At least one represented month has missing daily outputs'
print(f'Confirmed complete daily coverage for all {_month_counts.size} represented months')
# Resample to calendar-month means, then label each mean at mid-month.
IS2_SMOS_SMAP_monthly = IS2_SMOS_SMAP.sortby('time').resample(time='1ME').mean()
_fused_time = pd.to_datetime(IS2_SMOS_SMAP_monthly['time'].values)
_fused_start = _fused_time.to_period('M').to_timestamp(how='start')
_fused_mid = _fused_start + pd.Timedelta(days=14)
IS2_SMOS_SMAP_monthly = IS2_SMOS_SMAP_monthly.assign_coords(time=_fused_mid.values)
# Carry static lat/lon grid into monthly product (no time dimension)
if 'latitude' in IS2_SMOS_SMAP.coords:
IS2_SMOS_SMAP_monthly = IS2_SMOS_SMAP_monthly.assign_coords(latitude=IS2_SMOS_SMAP['latitude'])
if 'longitude' in IS2_SMOS_SMAP.coords:
IS2_SMOS_SMAP_monthly = IS2_SMOS_SMAP_monthly.assign_coords(longitude=IS2_SMOS_SMAP['longitude'])
IS2_SMOS_SMAP_monthly
Confirmed complete daily coverage for all 54 represented months
<xarray.Dataset> Size: 256MB
Dimensions: (x: 304, y: 448, time: 78)
Coordinates:
* x (x) float32 1kB -3.838e+06 ... 3.738e+06
* y (y) float32 2kB 5.838e+06 ... -5.338e+06
latitude (y, x) float32 545kB dask.array<chunksize=(448, 304), meta=np.ndarray>
longitude (y, x) float32 545kB dask.array<chunksize=(448, 304), meta=np.ndarray>
* time (time) datetime64[ns] 624B 2018-11-15...
Data variables: (12/13)
crs (time) float64 624B dask.array<chunksize=(1,), meta=np.ndarray>
freeboard (time, y, x) float32 42MB dask.array<chunksize=(1, 448, 304), meta=np.ndarray>
grid_cell_area (time, y, x) float32 42MB dask.array<chunksize=(1, 448, 304), meta=np.ndarray>
ice_thickness (time, y, x) float32 42MB dask.array<chunksize=(1, 448, 304), meta=np.ndarray>
ice_thickness_unc (time, y, x) float32 42MB dask.array<chunksize=(1, 448, 304), meta=np.ndarray>
mean_freeboard_regions_1_5 (time) float32 312B dask.array<chunksize=(1,), meta=np.ndarray>
... ...
mean_sea_ice_thickness_regions_1_5 (time) float32 312B dask.array<chunksize=(1,), meta=np.ndarray>
mean_snow_depth_regions_1_5 (time) float32 312B dask.array<chunksize=(1,), meta=np.ndarray>
region_mask (time, y, x) float32 42MB dask.array<chunksize=(1, 448, 304), meta=np.ndarray>
snow_depth (time, y, x) float32 42MB dask.array<chunksize=(1, 448, 304), meta=np.ndarray>
total_sea_ice_volume (time) float32 312B dask.array<chunksize=(1,), meta=np.ndarray>
total_sea_ice_volume_regions_1_5 (time) float32 312B dask.array<chunksize=(1,), meta=np.ndarray>
Attributes:
caa_uncertainty_caution: Canadian Arctic Archipelago (CAA; NSIDC-0780 ...
contact: Alek Petty (akpetty@umd.edu)
description: ICESat-2–SMOS–SMAP fused daily gridded produc...
history: Created 2026-07-29
inner_arctic_regions_note: Inner-Arctic scalars (mean_*_regions_1_5, tot...
n_timesteps: 1635
source: Combined from monthly IS2_interp_test_petty f...
time_coverage_end: 2025-04-30
time_coverage_start: 2018-11-01
training_season_caution: GP training is limited to a seasonal window (...def _growth_season_only(da):
"""Mask May--August while retaining gaps between growth seasons."""
return da.where(~da.time.dt.month.isin([5, 6, 7, 8]))
# Thickness maps are included in the compact nine-panel comparison below.
# March thickness maps are included in the compact nine-panel comparison below.
Freeboard comparison#
Inner-Arctic-Ocean (NSIDC regions 1–5) area-weighted mean time series for freeboard (V4: freeboard_int, fused: freeboard). The freeboard maps are consolidated with thickness and snow depth in the nine-panel comparison below.
# Inner Arctic Ocean area-weighted mean freeboard (NSIDC regions 1–5)
v4_fb = IS2_v4['freeboard_int']
v4_fb_valid = v4_fb.where(v4_fb > -1e10)
v4_mean_fb = area_weighted_spatial_mean(v4_fb_valid, mask_iao_v4, _area_v4)
fused_fb = IS2_SMOS_SMAP_monthly['freeboard']
fused_fb_valid = fused_fb.where(fused_fb > -1e10)
fused_mean_fb_daily = IS2_SMOS_SMAP['mean_freeboard_regions_1_5']
fused_mean_fb_monthly = is2smgpsit_monthly_midmonth(fused_mean_fb_daily)
# Freeboard maps are included in the compact nine-panel comparison below.
# March freeboard maps are included in the compact nine-panel comparison below.
Snow depth comparison#
Inner-Arctic-Ocean (NSIDC regions 1–5) area-weighted mean time series for snow depth (V4: snow_depth_int, fused: snow_depth). The cells below create one compact 3×3 comparison figure per example month, with panels (a)–(i).
# Inner Arctic Ocean area-weighted mean snow depth (NSIDC regions 1–5)
v4_sd = IS2_v4['snow_depth_int']
v4_sd_valid = v4_sd.where(v4_sd > -1e10)
v4_mean_sd = area_weighted_spatial_mean(v4_sd_valid, mask_iao_v4, _area_v4)
fused_sd = IS2_SMOS_SMAP_monthly['snow_depth']
fused_sd_valid = fused_sd.where(fused_sd > -1e10)
fused_mean_sd_daily = IS2_SMOS_SMAP['mean_snow_depth_regions_1_5']
fused_mean_sd_monthly = is2smgpsit_monthly_midmonth(fused_mean_sd_daily)
# Compact nine-panel comparison: thickness, freeboard, and snow depth
from utils.plotting_utils import plot_is2_v4_vs_fused_nine_panel
compare_month = np.datetime64('2019-10-15')
v4_comparison_fields = [
IS2_v4[name].sel(time=compare_month, method='nearest')
for name in ('ice_thickness_int', 'freeboard_int', 'snow_depth_int')
]
fused_comparison_fields = [
IS2_SMOS_SMAP_monthly[name].sel(time=compare_month, method='nearest')
for name in ('ice_thickness', 'freeboard', 'snow_depth')
]
fig, axes = plot_is2_v4_vs_fused_nine_panel(
dataarrays1=v4_comparison_fields,
dataarrays2=fused_comparison_fields,
)
fig.savefig(
'../paper/figures/inner_arctic_seasonal_comparison_oct2019.png',
dpi=300,
bbox_inches='tight',
)
plt.show()
# Repeat the compact comparison for March 2020
compare_month = np.datetime64('2020-03-15')
v4_comparison_fields = [
IS2_v4[name].sel(time=compare_month, method='nearest')
for name in ('ice_thickness_int', 'freeboard_int', 'snow_depth_int')
]
fused_comparison_fields = [
IS2_SMOS_SMAP_monthly[name].sel(time=compare_month, method='nearest')
for name in ('ice_thickness', 'freeboard', 'snow_depth')
]
fig, axes = plot_is2_v4_vs_fused_nine_panel(
dataarrays1=v4_comparison_fields,
dataarrays2=fused_comparison_fields,
)
fig.savefig(
'../paper/figures/inner_arctic_seasonal_comparison_march2020.png',
dpi=300,
bbox_inches='tight',
)
plt.show()
IS2SMGPSIT-V1 GP uncertainty (example months)#
Three-panel Arctic maps of the monthly mean daily GP thickness predictive standard deviation for October~2019, March~2020, and March~2021. These are averages of the daily posterior standard deviations, not posterior uncertainties calculated for a monthly-mean field.
from utils.plotting_utils import plot_is2smgpsit_thickness_unc_three_months
_da_unc_oct = IS2_SMOS_SMAP_monthly['ice_thickness_unc'].sel(time=np.datetime64('2019-10-15'), method='nearest')
_da_unc_mar_1 = IS2_SMOS_SMAP_monthly['ice_thickness_unc'].sel(time=np.datetime64('2020-03-15'), method='nearest')
_da_unc_mar_2 = IS2_SMOS_SMAP_monthly['ice_thickness_unc'].sel(time=np.datetime64('2021-03-15'), method='nearest')
_unc_fpath = '../paper/figures/is2smgpsit_v1_monthly_uncertainty_thickness_oct_mar.png'
plot_is2smgpsit_thickness_unc_three_months(
_da_unc_oct,
_da_unc_mar_1,
_da_unc_mar_2,
vmax=0.6,
)
plt.savefig(_unc_fpath, dpi=300, bbox_inches='tight')
plt.show()
Common-footprint V4–V1 Inner Arctic Ocean comparison#
We compare IS2SITMOGR4-V4 and IS2SMGPSIT-V1 over the same monthly footprint within the inner Arctic Ocean. The common footprint contains grid cells with valid thickness in both products. For context, we also show the IS2SMGPSIT-V1 mean over all valid inner-AO cells; an equivalent full-domain V4 mean is not available because V4 coverage is incomplete and varies by month. Standalone V1 seasonal cycles and anomalies are presented in notebook 11d.
# Inner Arctic Ocean (NSIDC regions 1–5)
# Fused: precomputed area-weighted scalars (mean_*_regions_1_5) from the Zarr.
# V4: grid-based area-weighted means (no equivalent scalar in IS2SITMOGR4).
from utils.read_data_utils import is2smgpsit_monthly_midmonth
v4_thick_iao = IS2_v4['ice_thickness_int'].where(mask_iao_v4)
fused_thick_iao = IS2_SMOS_SMAP_monthly['ice_thickness'].where(mask_iao_fused)
v4_mean_iao_available = area_weighted_spatial_mean(
v4_thick_iao.where(v4_thick_iao > -1e10), mask_iao_v4, _area_v4
)
fused_mean_daily_iao = IS2_SMOS_SMAP['mean_sea_ice_thickness_regions_1_5']
fused_mean_iao_full = is2smgpsit_monthly_midmonth(fused_mean_daily_iao)
v4_fb_iao = IS2_v4['freeboard_int'].where(mask_iao_v4)
fused_fb_iao = IS2_SMOS_SMAP_monthly['freeboard'].where(mask_iao_fused)
fused_mean_fb_iao_full = is2smgpsit_monthly_midmonth(
IS2_SMOS_SMAP['mean_freeboard_regions_1_5']
)
v4_mean_fb_iao_available = area_weighted_spatial_mean(
v4_fb_iao.where(v4_fb_iao > -1e10), mask_iao_v4, _area_v4
)
v4_sd_iao = IS2_v4['snow_depth_int'].where(mask_iao_v4)
fused_sd_iao = IS2_SMOS_SMAP_monthly['snow_depth'].where(mask_iao_fused)
fused_mean_sd_iao_full = is2smgpsit_monthly_midmonth(
IS2_SMOS_SMAP['mean_snow_depth_regions_1_5']
)
v4_mean_sd_iao_available = area_weighted_spatial_mean(
v4_sd_iao.where(v4_sd_iao > -1e10), mask_iao_v4, _area_v4
)
# Build common-footprint inner-AO monthly means (area-weighted; same mask for all three fields)
# Both products are already monthly with 15th-of-month timestamps
common_time = np.intersect1d(v4_mean_iao_available['time'].values, fused_mean_iao_full['time'].values)
v4_iao_common_vals = []
fused_iao_common_vals = []
v4_fb_iao_common_vals = []
fused_fb_iao_common_vals = []
v4_sd_iao_common_vals = []
fused_sd_iao_common_vals = []
common_area_fraction_vals = []
common_times = []
_iao_area_v4 = float(_area_v4.where(mask_iao_v4).sum())
for t in common_time:
v4_field = v4_thick_iao.sel(time=t)
fused_field = fused_thick_iao.sel(time=t)
v4_fb_field = v4_fb_iao.sel(time=t)
fused_fb_field = fused_fb_iao.sel(time=t)
v4_sd_field = v4_sd_iao.sel(time=t)
fused_sd_field = fused_sd_iao.sel(time=t)
# Common spatial mask: valid in both products (use thickness to define mask for all three)
mask = (
(v4_field > -1e-10) & np.isfinite(v4_field) &
(fused_field > -1e-10) & np.isfinite(fused_field)
)
if not mask.any().compute():
continue
def _aw_scalar(field, m, area):
f = field.where(m)
a = area.where(m)
return float((f * a).sum() / a.sum())
v4_iao_common_vals.append(_aw_scalar(v4_field, mask, _area_v4))
fused_iao_common_vals.append(_aw_scalar(fused_field, mask, _area_fused))
v4_fb_iao_common_vals.append(_aw_scalar(v4_fb_field, mask, _area_v4))
fused_fb_iao_common_vals.append(_aw_scalar(fused_fb_field, mask, _area_fused))
v4_sd_iao_common_vals.append(_aw_scalar(v4_sd_field, mask, _area_v4))
fused_sd_iao_common_vals.append(_aw_scalar(fused_sd_field, mask, _area_fused))
common_area_fraction_vals.append(float(_area_v4.where(mask).sum() / _iao_area_v4))
common_times.append(pd.to_datetime(t))
v4_mean_iao_common = xr.DataArray(
v4_iao_common_vals, coords={'time': common_times}, dims='time', name='v4_mean_iao_common'
)
fused_mean_iao_common = xr.DataArray(
fused_iao_common_vals, coords={'time': common_times}, dims='time', name='fused_mean_iao_common'
)
v4_mean_fb_iao_common = xr.DataArray(
v4_fb_iao_common_vals, coords={'time': common_times}, dims='time', name='v4_mean_fb_iao_common'
)
fused_mean_fb_iao_common = xr.DataArray(
fused_fb_iao_common_vals, coords={'time': common_times}, dims='time', name='fused_mean_fb_iao_common'
)
v4_mean_sd_iao_common = xr.DataArray(
v4_sd_iao_common_vals, coords={'time': common_times}, dims='time', name='v4_mean_sd_iao_common'
)
fused_mean_sd_iao_common = xr.DataArray(
fused_sd_iao_common_vals, coords={'time': common_times}, dims='time', name='fused_mean_sd_iao_common'
)
common_area_fraction = xr.DataArray(
common_area_fraction_vals, coords={'time': common_times}, dims='time', name='common_area_fraction'
)
_autumn_coverage = common_area_fraction.where(
common_area_fraction.time.dt.month.isin([9, 10]), drop=True
).to_series().rename('IAO area fraction')
print('Common-footprint coverage in September and October:')
print(_autumn_coverage.map(lambda value: f'{value:.1%}').to_string())
Common-footprint coverage in September and October:
time
2019-09-15 48.1%
2019-10-15 55.9%
2020-09-15 40.2%
2020-10-15 52.8%
2021-09-15 52.5%
2021-10-15 71.9%
2022-09-15 53.3%
2022-10-15 70.3%
2023-09-15 47.7%
2023-10-15 64.7%
2024-09-15 49.6%
2024-10-15 60.5%
# Reindex to continuous monthly grid so summer gaps appear as line breaks
_t0_iao = min(
v4_mean_iao_available.time.values.min(),
fused_mean_iao_full.time.values.min(),
)
_t1_iao = max(
v4_mean_iao_available.time.values.max(),
fused_mean_iao_full.time.values.max(),
)
_t0_iao, _t1_iao = pd.Timestamp(_t0_iao), pd.Timestamp(_t1_iao)
_grid_start_iao = (
pd.Timestamp(year=_t0_iao.year, month=9, day=1)
if _t0_iao.month >= 9
else pd.Timestamp(year=_t0_iao.year - 1, month=9, day=1)
)
all_months_iao = pd.date_range(_grid_start_iao, _t1_iao, freq='MS') + pd.Timedelta(days=14)
# Freeboard
fused_fb_iao_plot = _growth_season_only(fused_mean_fb_iao_full.reindex(time=all_months_iao))
v4_fb_iao_common_plot = _growth_season_only(v4_mean_fb_iao_common.reindex(time=all_months_iao))
fused_fb_iao_common_plot = _growth_season_only(fused_mean_fb_iao_common.reindex(time=all_months_iao))
# Snow depth
fused_sd_iao_plot = _growth_season_only(fused_mean_sd_iao_full.reindex(time=all_months_iao))
v4_sd_iao_common_plot = _growth_season_only(v4_mean_sd_iao_common.reindex(time=all_months_iao))
fused_sd_iao_common_plot = _growth_season_only(fused_mean_sd_iao_common.reindex(time=all_months_iao))
# Thickness
fused_iao_plot = _growth_season_only(fused_mean_iao_full.reindex(time=all_months_iao))
v4_iao_common_plot = _growth_season_only(v4_mean_iao_common.reindex(time=all_months_iao))
fused_iao_common_plot = _growth_season_only(fused_mean_iao_common.reindex(time=all_months_iao))
fig, axes = plt.subplots(3, 1, figsize=(10.5, 6.8), sharex=True)
# Row 0: freeboard
ax = axes[0]
v4_fb_iao_common_plot.plot(ax=ax, label='IS2SITMOGR4-V4 (common footprint)', color='k', marker='o', markersize=3, lw=1.5, ls='--')
fused_fb_iao_common_plot.plot(ax=ax, label='IS2SMGPSIT-V1 (common footprint)', color='b', marker='o', markersize=3, ls='--')
fused_fb_iao_plot.plot(ax=ax, label='IS2SMGPSIT-V1 (full IAO)', color='c', marker='x', markersize=3, lw=1)
ax.set_ylabel('Total freeboard (m)')
ax.set_xlabel('')
ax.tick_params(axis='x', labelbottom=False)
ax.grid(axis='y')
ax.set_ylim(0, 0.45)
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.legend(loc='upper right', ncol=3, frameon=False)
# Row 1: snow depth
ax = axes[1]
v4_sd_iao_common_plot.plot(ax=ax, label='IS2SITMOGR4-V4 (common footprint)', color='k', marker='o', markersize=3, lw=1.5, ls='--')
fused_sd_iao_common_plot.plot(ax=ax, label='IS2SMGPSIT-V1 (common footprint)', color='b', marker='o', markersize=3, ls='--')
fused_sd_iao_plot.plot(ax=ax, label='IS2SMGPSIT-V1 (full IAO)', color='c', marker='x', markersize=3, lw=1)
ax.set_ylabel('Snow depth (m)')
ax.tick_params(axis='x', labelbottom=False)
ax.grid(axis='y')
ax.set_xlabel('')
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
# Row 2: thickness
ax = axes[2]
v4_iao_common_plot.plot(ax=ax, label='IS2SITMOGR4-V4 (common footprint)', color='k', marker='o', markersize=3, lw=1.5, ls='--')
fused_iao_common_plot.plot(ax=ax, label='IS2SMGPSIT-V1 (common footprint)', color='b', marker='o', markersize=3, ls='--')
fused_iao_plot.plot(ax=ax, label='IS2SMGPSIT-V1 (full IAO)', color='c', marker='x', markersize=3, lw=1)
ax.set_ylabel('Sea ice thickness (m)')
ax.set_xlabel('Date')
ax.grid(axis='y')
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
# Shade summer months (May–Aug) on all rows
times_iao = pd.to_datetime(all_months_iao)
for ax in axes:
for year in np.unique(times_iao.year):
for month in [5, 6, 7, 8]:
start = pd.Timestamp(year=year, month=month, day=1)
end = pd.Timestamp(year=year, month=month + 1, day=1) if month < 12 else pd.Timestamp(year=year + 1, month=1, day=1)
if start >= times_iao.min() and start <= times_iao.max():
ax.axvspan(start, end, color='0.8', alpha=0.2, linewidth=0)
axes[2].xaxis.set_major_locator(mdates.MonthLocator(bymonth=[9, 1, 5], bymonthday=1))
axes[2].xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m'))
for ax in axes:
ax.tick_params(axis='y', labelsize=9)
ax.yaxis.label.set_size(9)
axes[2].tick_params(axis='x', labelsize=9)
axes[2].xaxis.label.set_size(9)
for label in axes[2].get_xticklabels():
label.set_rotation(30)
label.set_ha('right')
label.set_rotation_mode('anchor')
plt.tight_layout()
plt.savefig('../paper/figures/inner_arctic_seasonal_thickness_freeboard_snowdepth_comparison.png', dpi=300)
plt.show()
Summary#
The ICESat-2–SMOS–SMAP Zarr provides daily fields and can be loaded directly from S3 with
read_is2smgpsitv1_zarr().IS2SITMOGR4-V4 remains monthly; comparisons are done by resampling IS2SMGPSIT-V1 or selecting overlapping months.
Record-level thickness and volume time-series diagnostics are continued in notebook 11d.