9B. Sea Ice Supplementary

Explore the Sea Ice Data

Lawrence Livermore National Laboratory United States Department of Energy Earth System Grid Federation 2 - US Consortium

Summary

In this notebook, we are going to explore the dataset used for the PMP Sea Ice demo notebook. Let’s explore the sea ice data for fun!

Notebook Authors: Jiwoo Lee, Ana Ordonez, Paul Durack, Peter Gleckler (PCMDI, Lawrence Livermore National Laboratory)

Table of Contents

Note: Links to the sections work best when viewing this notebook via nbviewer. - 1. Environment setup - 2. Model Data * 2.1 Load data - 2.1.1 Open dataset - 2.1.2 Visualize the data * 2.2 Sea ice extent - 3. Reference Data * 3.1 Load data - 3.1.1 Open Reference Dataset for Arctic - 3.1.2 Open Reference Dataset for Antartica * 3.2 Sea ice extent - 4. Diagnostics: Climatology Annual Cycle - 5. Evaluation Metrics * 5.1 Mean Square Error (Annual Mean) * 5.2 Temporal Mean Square Error (Annual Cycle)

1. Environment setup

We will use multiple libraries for this analysis.

xCDAT: an open source Python tool built to make climate data analysis easy. xCDAT is an extension of xarray for data analysis on structured grids.
numpy: a dependency required to manage n-dimensional array data.
matplotlib and cartopy: required for data visualization.

[1]:

import xcdat as xc
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs

Go back to Top

2. Model Data

2.1 Load data

2.1.1 Open dataset

This demo uses one of the numerous CMIP6 models – E3SM-1-0. The Sea-Ice Area Percentage (Ocean Grid; ‘siconc’) and Grid-Cell Area for Ocean Variables (‘areacello’) variables are needed and can be found in the directories listed below. In addition, six other models are available that can augment the analyses in this demo:

/p/user_pub/pmp/demo/sea-ice/links_siconc
/p/user_pub/pmp/demo/sea-ice/links_area

This data can be found from ESGF (https://esgf-node.llnl.gov/). - e.g., Search query: https://aims2.llnl.gov/search?project=CMIP6&activeFacets={“experiment_id”:”historical”,”variable_id”:”siconc”}

Go back to Top

[2]:

!ls /p/user_pub/pmp/demo/sea-ice/links_siconc/E3SM-1-0/historical/r1i2p2f1/siconc

siconc_SImon_E3SM-1-0_historical_r1i2p2f1_gr_185001-185912.nc
siconc_SImon_E3SM-1-0_historical_r1i2p2f1_gr_186001-186912.nc
siconc_SImon_E3SM-1-0_historical_r1i2p2f1_gr_187001-187912.nc
siconc_SImon_E3SM-1-0_historical_r1i2p2f1_gr_188001-188912.nc
siconc_SImon_E3SM-1-0_historical_r1i2p2f1_gr_189001-189912.nc
siconc_SImon_E3SM-1-0_historical_r1i2p2f1_gr_190001-190912.nc
siconc_SImon_E3SM-1-0_historical_r1i2p2f1_gr_191001-191912.nc
siconc_SImon_E3SM-1-0_historical_r1i2p2f1_gr_192001-192912.nc
siconc_SImon_E3SM-1-0_historical_r1i2p2f1_gr_193001-193912.nc
siconc_SImon_E3SM-1-0_historical_r1i2p2f1_gr_194001-194912.nc
siconc_SImon_E3SM-1-0_historical_r1i2p2f1_gr_195001-195912.nc
siconc_SImon_E3SM-1-0_historical_r1i2p2f1_gr_196001-196912.nc
siconc_SImon_E3SM-1-0_historical_r1i2p2f1_gr_197001-197912.nc
siconc_SImon_E3SM-1-0_historical_r1i2p2f1_gr_198001-198912.nc
siconc_SImon_E3SM-1-0_historical_r1i2p2f1_gr_199001-199912.nc
siconc_SImon_E3SM-1-0_historical_r1i2p2f1_gr_200001-200912.nc
siconc_SImon_E3SM-1-0_historical_r1i2p2f1_gr_201001-201112.nc

[3]:

!ls /p/user_pub/pmp/demo/sea-ice/links_area/E3SM-1-0/

areacello_Ofx_E3SM-1-0_historical_r1i1p1f1_gr.nc

[4]:

# Load model dataset
ds = xc.open_mfdataset(
    "/p/user_pub/pmp/demo/sea-ice/links_siconc/E3SM-1-0/historical/r1i2p2f1/siconc/siconc_SImon_E3SM-1-0_historical_r1i2p2f1_*_*.nc"
)
area = xc.open_dataset(
    "/p/user_pub/pmp/demo/sea-ice/links_area/E3SM-1-0/areacello_Ofx_E3SM-1-0_historical_r1i1p1f1_gr.nc"
)

2024-02-09 12:08:07,978 [WARNING]: dataset.py(open_dataset:109) >> "No time coordinates were found in this dataset to decode. If time coordinates were expected to exist, make sure they are detectable by setting the CF 'axis' or 'standard_name' attribute (e.g., ds['time'].attrs['axis'] = 'T' or ds['time'].attrs['standard_name'] = 'time'). Afterwards, try decoding again with `xcdat.decode_time`."

See what is in the dataset:

[5]:

ds

[5]:

<xarray.Dataset>
Dimensions:    (time: 1944, bnds: 2, lat: 180, lon: 360)
Coordinates:
  * lat        (lat) float64 -89.5 -88.5 -87.5 -86.5 ... 86.5 87.5 88.5 89.5
  * lon        (lon) float64 0.5 1.5 2.5 3.5 4.5 ... 356.5 357.5 358.5 359.5
    type       |S7 b'sea_ice'
  * time       (time) object 1850-01-16 12:00:00 ... 2011-12-16 12:00:00
Dimensions without coordinates: bnds
Data variables:
    time_bnds  (time, bnds) object dask.array<chunksize=(1, 2), meta=np.ndarray>
    lat_bnds   (lat, bnds) float64 dask.array<chunksize=(180, 2), meta=np.ndarray>
    lon_bnds   (lon, bnds) float64 dask.array<chunksize=(360, 2), meta=np.ndarray>
    siconc     (time, lat, lon) float32 dask.array<chunksize=(1, 180, 360), meta=np.ndarray>
Attributes: (12/48)
    Conventions:            CF-1.7 CMIP-6.2
    activity_id:            CMIP
    branch_method:          standard
    branch_time_in_child:   0.0
    creation_date:          2023-01-30T21:36:30Z
    data contact:           Sam Stevenson (sstevenson@ucsb.edu)
    ...                     ...
    variant_label:          r1i2p2f1
    license:                CMIP6 model data produced by E3SM-Project is lice...
    cmor_version:           3.7.0
    references:             Stevenson, S., Huang, X., Zhao, Y., Di Lorenzo, E...
    version:                v20230811
    branch_time_in_parent:  3560.0

[6]:

area

Go back to Top

2.1.2 Visualize the data

Go back to Top

a. Quick look

Let’s take a quick look at a single timestep of the data:

[7]:

ds["siconc"].isel(time=0).plot()

[7]:

<matplotlib.collections.QuadMesh at 0x7fdd797dd7b0>

../_images/examples_Demo_9b_seaIce_data_explore_15_1.png

[8]:

ds["siconc"].isel(time=0).to_numpy()

[8]:

array([[      nan,       nan,       nan, ...,       nan,       nan,
              nan],
       [      nan,       nan,       nan, ...,       nan,       nan,
              nan],
       [      nan,       nan,       nan, ...,       nan,       nan,
              nan],
       ...,
       [99.87347 , 99.88154 , 99.88288 , ..., 99.86003 , 99.861336,
        99.86454 ],
       [99.656136, 99.65681 , 99.657486, ..., 99.655624, 99.65451 ,
        99.65535 ],
       [99.54839 , 99.55058 , 99.554344, ..., 99.54922 , 99.54881 ,
        99.54852 ]], dtype=float32)

b. Snapshot on map

Viewing the same data from another angle: - Example mapping scripts are adapted and revised from https://docs.xarray.dev/en/latest/user-guide/plotting.html#maps.

Go back to Top

[9]:

from cartopy.feature import LAND as cartopy_land
from cartopy.feature import OCEAN as cartopy_ocean

[10]:

ax = plt.axes(projection=ccrs.Orthographic(-80, 35))

ax.set_global()

# Color over the ocean
ax.add_feature(cartopy_ocean, zorder=1, edgecolor="none", facecolor="lightblue")

# Data
ds["siconc"].isel(time=0).plot.pcolormesh(ax=ax, transform=ccrs.PlateCarree(), cmap="Blues_r", zorder=3)

# Color over the land
ax.add_feature(cartopy_land, zorder=2, edgecolor="k", facecolor="lightgrey")

# Add coastlines
ax.coastlines(color="black", zorder=4)

ax.set_title("sea ice concentration over the Arctic region")

[10]:

Text(0.5, 1.0, 'sea ice concentration over the Arctic region')

../_images/examples_Demo_9b_seaIce_data_explore_19_1.png

Let’s plot the data using another angle, zooming in on the Arctic region:

Additional matplotlib colorshemes can be found here.

[11]:

# Generate a map
proj = ccrs.NorthPolarStereo()
ax = plt.subplot(111, projection=proj)
ax.set_global()

# Color over the ocean
ax.add_feature(cartopy_ocean, zorder=1, edgecolor="none", facecolor="lightblue")

# Color over the land
ax.add_feature(cartopy_land, zorder=2, edgecolor="k", facecolor="lightgrey")

# Data
ds["siconc"].isel(time=0).plot.pcolormesh(
    ax=ax,
    x="lon",
    y="lat",
    transform=ccrs.PlateCarree(),
    cmap="Blues_r",
    cbar_kwargs={"label": "ice concentration (%)"},
    zorder=3,
)

# Add coastlines
ax.coastlines(color="black", zorder=4)

ax.set_extent([-180, 180, 43, 90], ccrs.PlateCarree())
ax.set_title("sea ice concentration over the Arctic region")

[11]:

Text(0.5, 1.0, 'sea ice concentration over the Arctic region')

../_images/examples_Demo_9b_seaIce_data_explore_21_1.png

c. See the data evolving in time

Arctic:

Go back to Top

[12]:

proj_plot = ccrs.Orthographic(0, 90)

years_to_show = [2000]  # You can edit this to check another year(s). It could include multiple years, e.g., [2000, 2001, 2002]


p = ds['siconc'].sel(time = ds.time.dt.year.isin(years_to_show)).plot(x='lon', y='lat',
                                                                 transform=ccrs.PlateCarree(),
                                                                 aspect=ds.sizes["lon"] / ds.sizes["lat"],  # for a sensible figsize
                                                                 subplot_kws={"projection": proj_plot},
                                                                 col='time', col_wrap=6, robust=True, cmap='Blues_r', zorder=3)

# We have to set the map's options on all four axes
for ax,i in zip(p.axs.flat,  ds.time.sel(time = ds.time.dt.year.isin(years_to_show)).values):
    ax.set_title(i.strftime("%B %Y"), fontsize=18)
    # Color over the ocean
    ax.add_feature(cartopy_ocean, zorder=1, edgecolor="none", facecolor="lightblue")
    # Color over the land
    ax.add_feature(cartopy_land, zorder=2, edgecolor="k", facecolor="lightgrey")
    # Draw coastlines
    ax.coastlines(zorder=4)

../_images/examples_Demo_9b_seaIce_data_explore_23_0.png

Go back to Top

Antartic:

[13]:

proj_plot = ccrs.Orthographic(0, -90)

#years_to_show = [2000]  # You can uncomment and edit this to check another year(s). It could include multiple years, e.g., [2000, 2001, 2002]

p = ds['siconc'].sel(time = ds.time.dt.year.isin(years_to_show)).plot(x='lon', y='lat',
                                                                 transform=ccrs.PlateCarree(),
                                                                 aspect=ds.sizes["lon"] / ds.sizes["lat"],  # for a sensible figsize
                                                                 subplot_kws={"projection": proj_plot},
                                                                 col='time', col_wrap=6, robust=True, cmap='Blues_r', zorder=3)

# We have to set the map's options on all four axes
for ax,i in zip(p.axs.flat,  ds.time.sel(time = ds.time.dt.year.isin(years_to_show)).values):
    ax.set_title(i.strftime("%B %Y"), fontsize=18)
    # Color over the ocean
    ax.add_feature(cartopy_ocean, zorder=1, edgecolor="none", facecolor="lightblue")
    # Color over the land
    ax.add_feature(cartopy_land, zorder=2, edgecolor="k", facecolor="lightgrey")
    # Draw coastlines
    ax.coastlines(zorder=4)

../_images/examples_Demo_9b_seaIce_data_explore_26_0.png

2.2 Sea ice extent

Go back to Top

Below, we will calculate the sea ice extent from sea ice concentration.

[14]:

sea_ice_extent = (ds.siconc.where(ds.lat > 0) * 1e-2 * area.areacello * 1e-6).where(ds.siconc > 15).sum(("lat", "lon"))

Note for the above line - where siconc > 15: to consider sea ice extent instead of total sea ice area - multiply 1e-2: to convert percentage (%) to fraction

To learn more about the sea ice extent, see this article from NSIDC.

Go back to Top

Visualize the sea ice extent time series:

[15]:

sea_ice_extent.sel({"time": slice("1981-01-01", "2010-12-31")}).plot()

[15]:

[<matplotlib.lines.Line2D at 0x7fd9e254e5c0>]

../_images/examples_Demo_9b_seaIce_data_explore_33_1.png

Let’s add the sea_ice_extent dataArray to the dataset (ds) for later use:

[16]:

ds["extent"] = sea_ice_extent

[17]:

ds

[17]:

<xarray.Dataset>
Dimensions:    (time: 1944, bnds: 2, lat: 180, lon: 360)
Coordinates:
  * lat        (lat) float64 -89.5 -88.5 -87.5 -86.5 ... 86.5 87.5 88.5 89.5
  * lon        (lon) float64 0.5 1.5 2.5 3.5 4.5 ... 356.5 357.5 358.5 359.5
    type       |S7 b'sea_ice'
  * time       (time) object 1850-01-16 12:00:00 ... 2011-12-16 12:00:00
Dimensions without coordinates: bnds
Data variables:
    time_bnds  (time, bnds) object dask.array<chunksize=(1, 2), meta=np.ndarray>
    lat_bnds   (lat, bnds) float64 dask.array<chunksize=(180, 2), meta=np.ndarray>
    lon_bnds   (lon, bnds) float64 dask.array<chunksize=(360, 2), meta=np.ndarray>
    siconc     (time, lat, lon) float32 dask.array<chunksize=(1, 180, 360), meta=np.ndarray>
    extent     (time) float32 dask.array<chunksize=(1,), meta=np.ndarray>
Attributes: (12/48)
    Conventions:            CF-1.7 CMIP-6.2
    activity_id:            CMIP
    branch_method:          standard
    branch_time_in_child:   0.0
    creation_date:          2023-01-30T21:36:30Z
    data contact:           Sam Stevenson (sstevenson@ucsb.edu)
    ...                     ...
    variant_label:          r1i2p2f1
    license:                CMIP6 model data produced by E3SM-Project is lice...
    cmor_version:           3.7.0
    references:             Stevenson, S., Huang, X., Zhao, Y., Di Lorenzo, E...
    version:                v20230811
    branch_time_in_parent:  3560.0

Go back to Top

4. Diagnostics: Climatology Annual Cycle

Let’s analyze the climatology of Sea-Ice extent, focusing on its annual cycle.

We will use xCDAT’s `temporal.climatology <https://xcdat.readthedocs.io/en/latest/generated/xarray.Dataset.temporal.climatology.html>`__ function to calculate the monthly climatology field.

Go back to Top

[26]:

obs_ds_clim = obs_ds.temporal.climatology(
    "extent", freq="month", reference_period=("1981-01-01", "2010-12-31"))

[27]:

ds_clim = ds.temporal.climatology(
    "extent", freq="month", reference_period=("1981-01-01", "2010-12-31"))

[28]:

fig, ax = plt.subplots()

# Unify the time axis to plot lines on a same figure
ds_clim["time"] = [x for x in range(1, 13)]
obs_ds_clim["time"] = [x for x in range(1, 13)]

# Add lines
ds_clim["extent"].plot(ax=ax, label="Model (E3SM-1-0)")
obs_ds_clim["extent"].plot(ax=ax, label="OBS (OSI-SAF)")

# Attach legend
plt.title("Arctic climatological sea ice extent\n1981-2010")
plt.xlabel("month")
plt.ylabel("Extent (km${^2}$)")
plt.legend(loc="upper right", fontsize=9)

[28]:

<matplotlib.legend.Legend at 0x7fda0079d960>

../_images/examples_Demo_9b_seaIce_data_explore_56_1.png

Go back to Top

5. Evaluation Metrics

The term “metric” indicates a score or statistics that can represent a model’s performance in reproducing observed features. A “metric” can be defined in many different ways, but the primary purpose of metrics is to summarize the performances of many different models and provide a framework for benchmarking model realism.

In this notebook, we define Mean Square Error (MSE) and Temporal MSE as metrics.

5.1 Mean Square Error (Annual Mean)

Go back to Top

[29]:

mse = (ds_clim["extent"].data.mean() - obs_ds_clim["extent"].data.mean()) ** 2 * 1e-12

[30]:

print(f"Mean Square Error: {float(mse)} 10^12 km^4")

Mean Square Error: 1.9092291476097731 10^12 km^4

5.2 Temporal Mean Square Error (Annual Cycle)

Go back to Top

[31]:

tmse = np.sum(((ds_clim["extent"].data - obs_ds_clim["extent"].data) ** 2)) / len(ds_clim["extent"]) * 1e-12

[32]:

print(f"Temporal Mean Square Error: {float(tmse)} 10^12 km^4")

Temporal Mean Square Error: 3.808760315117328 10^12 km^4

Go back to Top

9B. Sea Ice Supplementary

Explore the Sea Ice Data

1. Environment setup

2. Model Data

2.1 Load data

2.1.1 Open dataset

2.1.2 Visualize the data

a. Quick look

b. Snapshot on map

c. See the data evolving in time

2.2 Sea ice extent

3. Reference Data

3.1 Load data

3.1.1 Open Reference Dataset for Arctic

3.1.2 Open Reference Dataset for Antartica

3.2 Sea ice extent

4. Diagnostics: Climatology Annual Cycle

5. Evaluation Metrics

5.1 Mean Square Error (Annual Mean)

5.2 Temporal Mean Square Error (Annual Cycle)