Sea Ice Classification and local Leads/floes PDF analysis
This tutorial provides a step-by-step guide to analyzing polar region using SWOT data. It focuses on visualizing sea ice classification, binning ice concentration, and estimating ice freeboard (snow + ice thickness) using 2D histograms (Probability Density Functions, or PDFs).
Tutorial Objectives
Subset Swot data with cycle/passes numbers, and a geographical area selection, using
altimetry_downloader_avisoOpen data using
altimetry.ioPerform coordinates reprojection to ESPG:6931 using
pyprojVisualise sea ice classification using
matplotlib+cartopyPerform 2D binning on ice concentration, using
pyinterpEstimate ice freeboard using 2D histograms using
numpy+matplotlib
Note
Required environment to run this notebook:
xarray+numpypyprojpyinterpmatplotlib<3.11+cartopyaltimetry_downloader_aviso: see documentation.altimetry.io: available here.
Import + code
import altimetry_downloader_aviso as dl_aviso
from altimetry.io import AltimetryData, FileCollectionSource
import numpy as np
import xarray as xr
from pathlib import Path
import pyinterp
from pyproj import Proj, Transformer, transform
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from mpl_toolkits.axes_grid1 import make_axes_locatable
import matplotlib.axes as maxes
import logging
logging.basicConfig(level=logging.INFO)
def reproj_coords(ds, crs):
"""
Transforms coordinates to a new crs
"""
lat = ds["latitude"].values
idx = np.where(lat >= LAT_BOUND)
(_id, _) = idx
sel = {'num_lines': slice(_id[0], _id[-1] + 1)}
ds_crop = ds.isel(**sel)
transformer = Transformer.from_crs(4326, crs)
x_proj, y_proj = transformer.transform(
ds_crop.latitude.data,
ds_crop.longitude.data
)
ds_crop = ds_crop.assign(
x_proj=(ds_crop.latitude.dims, x_proj),
y_proj=(ds_crop.latitude.dims, y_proj)
)
return ds_crop
def set_carto(extent=None, nrows=1, ncols=1, projection=ccrs.NorthPolarStereo(), central_lon=0, figsize=(14, 12)):
fig, axs = plt.subplots(
nrows=nrows,
ncols=ncols,
subplot_kw={'projection':projection},
figsize=figsize,
dpi=80
)
gl = axs.gridlines(draw_labels=True)
gl.top_labels = gl.right_labels = False
axs.coastlines(color="black", lw=0.5)
divider = make_axes_locatable(axs)
ax_cb = divider.new_horizontal(size="5%", pad=0.1, axes_class=plt.Axes)
fig = plt.gcf()
fig.add_axes(ax_cb)
if extent:
axs.set_extent(extent, crs=ccrs.PlateCarree())
return fig, axs, ax_cb
def get_bins(lat_bound: float, resolution: float):
"""
Generates a regular grid (xc_bins, yc_bins) in North polar stereographical projection
"""
proj_polar = Proj(proj='stere', lat_0=90, lat_ts=70, lon_0=-45, datum='WGS84', units='m')
x_min, y_min = proj_polar(0, lat_bound)
x_max, y_max = 0, 0
dist_max = np.abs(y_min)
extent = dist_max
num_bins = int(np.ceil((2 * extent) / resolution)) + 1
xc_bins = np.linspace(-extent, extent, num_bins)
yc_bins = np.linspace(-extent, extent, num_bins)
return xc_bins, yc_bins
def plot_ice_conc(Xgeo, Ygeo, bin2d, crs, extent, central_latitude, fig_title, figsize=(10,15)):
fig = plt.figure(figsize=figsize)
ax = plt.axes(projection=ccrs.LambertAzimuthalEqualArea(central_latitude=central_latitude))
ax.coastlines()
ax.set_extent(extent, crs=ccrs.PlateCarree())
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
linewidth=2, color='gray', alpha=0.5, linestyle='--')
im = ax.pcolormesh(
Xgeo,
Ygeo,
(1-bin2d.mean()).T*100,
transform=ccrs.epsg(crs),
vmin=0,
vmax=100,
cmap="Spectral_r"
)
divider = make_axes_locatable(ax)
cax = divider.append_axes("bottom", "3%", pad='7%', axes_class=maxes.Axes)
plt.colorbar(im, label='Sea ice concentration [%]', cax=cax, location='bottom')
ax.set_title(f"{fig_title} \n Mean sure leads (3) to sure leads and floes (3 and 0)\n1 and 2 values are ignored")
Parameters
output_dir= Path.home() / "TMP_DATA" / "subsets"
cycle_number=[13]
half_orbits = [153, 155, 157, 181, 183, 209, 211, 237, 239, 265, 267, 293, 295, 319]
box=(-137, 71, -122, 78)
# coverage and resolution
LAT_BOUND = 60 # min absolute latitude
RESOLUTION = 12500. # pixels size in meters
Download data using altimetry_downloader_aviso
dl_aviso.subset(
'SWOT_L3_LR_SSH_Unsmoothed',
output_dir=output_dir,
cycle_number=cycle_number,
pass_number=half_orbits,
selected_variables=["time", "latitude", "longitude", "quality_flag", "ssha_unedited"],
box=box
)
INFO:altimetry_downloader_aviso.catalog_client.client:Fetching products from Aviso's catalog...
INFO:altimetry_downloader_aviso.catalog_client.granule_discoverer:Filtering SWOT_L3_LR_SSH_Unsmoothed product with filters {'cycle_number': [13], 'pass_number': [153, 155, 157, 181, 183, 209, 211, 237, 239, 265, 267, 293, 295, 319]}...
INFO:altimetry_downloader_aviso.core:Subsetting 0 file(s)...
['/home/atonneau/TMP_DATA/subsets/SWOT_L3_LR_SSH_Unsmoothed_013_153_20240402T005441_20240402T014608_v3.0.nc',
'/home/atonneau/TMP_DATA/subsets/SWOT_L3_LR_SSH_Unsmoothed_013_155_20240402T023735_20240402T032902_v3.0.nc',
'/home/atonneau/TMP_DATA/subsets/SWOT_L3_LR_SSH_Unsmoothed_013_157_20240402T042028_20240402T051153_v3.0.nc',
'/home/atonneau/TMP_DATA/subsets/SWOT_L3_LR_SSH_Unsmoothed_013_181_20240403T005512_20240403T014639_v3.0.nc',
'/home/atonneau/TMP_DATA/subsets/SWOT_L3_LR_SSH_Unsmoothed_013_183_20240403T023806_20240403T032932_v3.0.nc',
'/home/atonneau/TMP_DATA/subsets/SWOT_L3_LR_SSH_Unsmoothed_013_209_20240404T005543_20240404T014710_v3.0.nc',
'/home/atonneau/TMP_DATA/subsets/SWOT_L3_LR_SSH_Unsmoothed_013_211_20240404T023836_20240404T033001_v3.0.nc',
'/home/atonneau/TMP_DATA/subsets/SWOT_L3_LR_SSH_Unsmoothed_013_237_20240405T005614_20240405T014741_v3.0.nc',
'/home/atonneau/TMP_DATA/subsets/SWOT_L3_LR_SSH_Unsmoothed_013_239_20240405T023907_20240405T033032_v3.0.nc',
'/home/atonneau/TMP_DATA/subsets/SWOT_L3_LR_SSH_Unsmoothed_013_265_20240406T005645_20240406T014812_v3.0.nc',
'/home/atonneau/TMP_DATA/subsets/SWOT_L3_LR_SSH_Unsmoothed_013_267_20240406T023938_20240406T033105_v3.0.nc',
'/home/atonneau/TMP_DATA/subsets/SWOT_L3_LR_SSH_Unsmoothed_013_293_20240407T005716_20240407T014843_v3.0.nc',
'/home/atonneau/TMP_DATA/subsets/SWOT_L3_LR_SSH_Unsmoothed_013_295_20240407T024010_20240407T033136_v3.0.nc',
'/home/atonneau/TMP_DATA/subsets/SWOT_L3_LR_SSH_Unsmoothed_013_319_20240407T231454_20240408T000524_v3.0.nc']
Open data using altimetry.io
Open data source
alti_data = AltimetryData(
source=FileCollectionSource(
path=output_dir,
ftype="SWOT_L3_LR_SSH",
subset="Unsmoothed",
),
)
Query data
ds = alti_data.query_orbit(
cycle_number=cycle_number,
pass_number=half_orbits,
variables=["time", "latitude", "longitude", "quality_flag", "ssha_unedited"],
)
ds
INFO:fcollections.core._filesdb:Picked subset {'version': '3.0', 'subset': <ProductSubset.Unsmoothed: 4>}
INFO:fcollections.core._readers:Files to read: 14
<xarray.Dataset>
Dimensions: (num_lines: 18802, num_pixels: 519)
Coordinates:
time (num_lines) datetime64[ns] dask.array<chunksize=(2682,), meta=np.ndarray>
latitude (num_lines, num_pixels) float64 dask.array<chunksize=(2682, 519), meta=np.ndarray>
longitude (num_lines, num_pixels) float64 dask.array<chunksize=(2682, 519), meta=np.ndarray>
Dimensions without coordinates: num_lines, num_pixels
Data variables:
quality_flag (num_lines, num_pixels) uint8 dask.array<chunksize=(2682, 519), meta=np.ndarray>
ssha_unedited (num_lines, num_pixels) float64 dask.array<chunksize=(2682, 519), meta=np.ndarray>
cycle_number (num_lines) uint8 13 13 13 13 13 13 13 ... 13 13 13 13 13 13
pass_number (num_lines) uint16 153 153 153 153 153 ... 319 319 319 319
Attributes: (12/43)
Conventions: CF-1.9
Metadata_Conventions: Unidata Dataset Discovery v1.0
cdm_data_type: Swath
comment: Sea Surface Height measured by Altimetry
geospatial_lat_units: degrees_north
geospatial_lon_units: degrees_east
... ...
geospatial_lon_max: 275.949791
date_modified: 2026-08-13T18:44:42Z
history: 2026-08-13T18:44:42Z: Created by DUACS K...
date_created: 2026-08-13T18:44:42Z
date_issued: 2026-08-13T18:44:42Z
temporality: reproc1. Coordinates reprojection to ESPG:6931
# Arctic
crs = 6931
ds_proj = reproj_coords(ds, crs)
ds_proj
<xarray.Dataset>
Dimensions: (num_lines: 18802, num_pixels: 519)
Coordinates:
time (num_lines) datetime64[ns] dask.array<chunksize=(2682,), meta=np.ndarray>
latitude (num_lines, num_pixels) float64 dask.array<chunksize=(2682, 519), meta=np.ndarray>
longitude (num_lines, num_pixels) float64 dask.array<chunksize=(2682, 519), meta=np.ndarray>
Dimensions without coordinates: num_lines, num_pixels
Data variables:
quality_flag (num_lines, num_pixels) uint8 dask.array<chunksize=(2682, 519), meta=np.ndarray>
ssha_unedited (num_lines, num_pixels) float64 dask.array<chunksize=(2682, 519), meta=np.ndarray>
cycle_number (num_lines) uint8 13 13 13 13 13 13 13 ... 13 13 13 13 13 13
pass_number (num_lines) uint16 153 153 153 153 153 ... 319 319 319 319
x_proj (num_lines, num_pixels) float64 -1.472e+06 ... -1.815e+06
y_proj (num_lines, num_pixels) float64 1.514e+06 ... 1.016e+06
Attributes: (12/43)
Conventions: CF-1.9
Metadata_Conventions: Unidata Dataset Discovery v1.0
cdm_data_type: Swath
comment: Sea Surface Height measured by Altimetry
geospatial_lat_units: degrees_north
geospatial_lon_units: degrees_east
... ...
geospatial_lon_max: 275.949791
date_modified: 2026-08-13T18:44:42Z
history: 2026-08-13T18:44:42Z: Created by DUACS K...
date_created: 2026-08-13T18:44:42Z
date_issued: 2026-08-13T18:44:42Z
temporality: reproc2. Visualization of Sea Ice Classification with Flag Application
Learn how to apply quality flags to filter and visualize sea ice classification data.
Understand how to identify and exclude artifacts for accurate analysis.
extent = [-109, -151, 71, 73]
fig, ax, cb = set_carto(extent)
cb.remove()
# 19 = "ice_unsure"
# 20 = "ice"
# 101 = "not_on_sea"
# 102 = "no_data"
flag_sea_ice = xr.where(
(ds_proj.quality_flag==19) | (ds_proj.quality_flag==20),
0,
xr.where((ds_proj.quality_flag==101) | (ds_proj.quality_flag==102),
np.nan, # not_on_sea, no_data -> nan
1
)
)
ax.pcolormesh(
ds_proj.longitude,
ds_proj.latitude,
flag_sea_ice,
transform=ccrs.PlateCarree(),
cmap="Blues",
vmin=0,
vmax=2,
)
fig = plt.gcf()
3. Binning and Visualization of Sea Ice Concentration
Perform spatial binning to visualize sea ice concentration.
Generate maps to highlight areas of interest for further analysis.
flag_sea_ice = xr.where(
(ds_proj.quality_flag==20),
0, # ice -> 0
xr.where((ds_proj.quality_flag==101) | (ds_proj.quality_flag==102),
np.nan, # not_on_sea, no_data -> nan
1 # ice_unsure -> 1
)
)
xc_bins, yc_bins = get_bins(LAT_BOUND, RESOLUTION)
x = ds_proj.x_proj.data.ravel()
y = ds_proj.y_proj.data.ravel()
z = flag_sea_ice.data.ravel()
x_axis = pyinterp.Axis(np.array(xc_bins, dtype='float64'))
y_axis = pyinterp.Axis(np.array(yc_bins, dtype='float64'))
bin2d_flag_ice_conc = pyinterp.Binning2D(x_axis, y_axis, dtype=np.dtype('float64'))
bin2d_flag_ice_conc.push(x, y, np.asarray(z), simple=False)
Plot result
central_latitude = 74
Xgeo, Ygeo = np.meshgrid(xc_bins, yc_bins)
fig_title = f'Sea Ice concentration, NH'
plot_ice_conc(Xgeo, Ygeo, bin2d_flag_ice_conc, crs, extent, central_latitude, fig_title)
4. 2D Histogram (PDF) Analysis for Leads and Floes
Calculate 2D histograms (PDFs) for sea ice leads (cracks) and floes (ice sheets).
Compare the PDFs of leads and floes to derive differences in elevation.
Use the difference between the two PDFs to estimate the freeboard (snow + ice thickness) from SWOT KaRIn, with an expected order of magnitude of ~20 cm.
ssha_leads = xr.where(flag_sea_ice==1, ds_proj.ssha_unedited, np.nan)
ssha_floes = xr.where(flag_sea_ice==0, ds_proj.ssha_unedited, np.nan)
hist, bins = np.histogram(ssha_leads, bins=np.arange(-0.2,0.5,0.01), density=True)
hist_floes, bins = np.histogram(ssha_floes, bins=np.arange(-0.2,0.5,0.01), density=True)
plt.figure()
plt.plot(bins[1:],hist, label="leads")
plt.plot(bins[1:],hist_floes, label="floes")
plt.legend()
<matplotlib.legend.Legend at 0x761c60090680>
