Sea ice analysis

Work in Progress

This notebook is currently under development. Content and results may change!

[1]:

%load_ext autoreload
%autoreload

[2]:

%matplotlib inline
import warnings, yaml
warnings.filterwarnings("ignore")
from mom6_tools.MOM6grid import MOM6grid
from mom6_tools.m6toolbox import cime_xmlquery
from ncar_jobqueue import NCARCluster
from dask.distributed import Client
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature
import numpy as np
import xarray as xr
from mom6_tools.m6plot import polarplot

Basemap module not found. Some regional plots may not function properly

ERROR 1: PROJ: proj_create_from_database: Open of /glade/work/gmarques/conda-envs/mom6-tools/share/proj failed

[3]:

# Read in the yaml file
diag_config_yml_path = "diag_config.yml"
diag_config_yml = yaml.load(open(diag_config_yml_path,'r'), Loader=yaml.Loader)

[4]:

caseroot = diag_config_yml['Case']['CASEROOT']
casename = cime_xmlquery(caseroot, 'CASE')
DOUT_S = cime_xmlquery(caseroot, 'DOUT_S')
rundir = cime_xmlquery(caseroot, 'RUNDIR')

if DOUT_S:
  OUTDIR = cime_xmlquery(caseroot, 'DOUT_S_ROOT')+'/ocn/hist/'
  OUTDIR_ice = cime_xmlquery(caseroot, 'DOUT_S_ROOT')+'/ice/hist/'
else:
  OUTDIR = cime_xmlquery(caseroot, 'RUNDIR')
  OUTDIR_ice = OUTDIR

print('Rundir directory is:', rundir)
print('Ice history files are at:', OUTDIR_ice)
print('Casename is:', casename)

Rundir directory is: /glade/derecho/scratch/gmarques/g.e30_a03c.GJRAv4.TL319_t232_wgx3_hycom1_N75.2024.079/run
Ice history files are at: /glade/derecho/scratch/gmarques/archive/g.e30_a03c.GJRAv4.TL319_t232_wgx3_hycom1_N75.2024.079/ice/hist/
Casename is: g.e30_a03c.GJRAv4.TL319_t232_wgx3_hycom1_N75.2024.079

[5]:

# create an empty class object
class args:
  pass# create an empty class object

# load avg dates
avg = diag_config_yml['Avg']

args.start_date = avg['start_date']
args.end_date = avg['end_date']
args.casename = casename
args.ice = casename+diag_config_yml['Fnames']['ice']
args.geom =   casename+diag_config_yml['Fnames']['geom']
args.static = casename+diag_config_yml['Fnames']['static']
args.savefigs = False
args.time_series = True
args.nw = 6 # requesting 6 workers

[6]:

variables = ['hi', 'aice']

def preprocess(ds):
    '''Return the dataset with variables'''
    return ds[variables]

%time ds = xr.open_mfdataset(OUTDIR_ice+args.ice, preprocess=preprocess)

CPU times: user 2min 4s, sys: 2min 27s, total: 4min 32s
Wall time: 9min 7s

[7]:

ds

[7]:

<xarray.Dataset> Size: 2GB
Dimensions:  (time: 732, nj: 480, ni: 540)
Coordinates:
  * time     (time) object 6kB 0001-01-16 12:00:00 ... 0061-12-16 12:00:00
    TLON     (nj, ni) float32 1MB dask.array<chunksize=(480, 540), meta=np.ndarray>
    TLAT     (nj, ni) float32 1MB dask.array<chunksize=(480, 540), meta=np.ndarray>
    ULON     (nj, ni) float32 1MB dask.array<chunksize=(480, 540), meta=np.ndarray>
    ULAT     (nj, ni) float32 1MB dask.array<chunksize=(480, 540), meta=np.ndarray>
    NLON     (nj, ni) float32 1MB dask.array<chunksize=(480, 540), meta=np.ndarray>
    NLAT     (nj, ni) float32 1MB dask.array<chunksize=(480, 540), meta=np.ndarray>
    ELON     (nj, ni) float32 1MB dask.array<chunksize=(480, 540), meta=np.ndarray>
    ELAT     (nj, ni) float32 1MB dask.array<chunksize=(480, 540), meta=np.ndarray>
Dimensions without coordinates: nj, ni
Data variables:
    hi       (time, nj, ni) float32 759MB dask.array<chunksize=(1, 480, 540), meta=np.ndarray>
    aice     (time, nj, ni) float32 759MB dask.array<chunksize=(1, 480, 540), meta=np.ndarray>
Attributes:
    title:               g.e30_a03c.GJRAv4.TL319_t232_wgx3_hycom1_N75.2024.079
    contents:            Diagnostic and Prognostic Variables
    source:              CICE Sea Ice Model, unknown_version_name
    comment:             All years have exactly 365 days
    comment2:            File written on model date 00010201
    comment3:            seconds elapsed into model date:      0
    time_period_freq:    month_1
    time_axis_position:  middle
    conventions:         CF-1.0
    history:             This dataset was created on 2024-10-02 at 15:04
    io_flavor:           io_pio2 cdf1

[8]:

print('\n Selecting data between {} and {}...'.format(args.start_date, args.end_date))
%time ds_sel = ds.sel(time=slice(args.start_date, args.end_date))


 Selecting data between 0031-01-01 and 0062-01-01...
CPU times: user 10.4 ms, sys: 0 ns, total: 10.4 ms
Wall time: 44.3 ms

[9]:

print('\n Computing time mean...')
%time ds_sel_mean = ds_sel.mean('time').compute()


 Computing time mean...
CPU times: user 1min 18s, sys: 2min 53s, total: 4min 12s
Wall time: 9min 18s

[10]:

print('Computing monthly climatology...')
%time  ds_monthly = ds_sel.groupby("time.month").mean('time').compute()

Computing monthly climatology...
CPU times: user 1min 29s, sys: 3min 34s, total: 5min 4s
Wall time: 10min 54s

[11]:

# load obs
obs_path = '/glade/work/gmarques/cesm/datasets/Seaice/'
obs_NH = xr.open_dataset(obs_path+'OBS_NSIDC_sat_NH_T2Ms_sic.nc')
obs_SH = xr.open_dataset(obs_path+'OBS_NSIDC_sat_SH_T2Ms_sic.nc')

[12]:

obs_NH_mean = obs_NH.mean('time').compute()
obs_SH_mean = obs_SH.mean('time').compute()

[13]:

obs_NH_monthly = obs_NH.groupby("time.month").mean('time').compute()
obs_SH_monthly = obs_SH.groupby("time.month").mean('time').compute()

Sea ice thickness

Annual mean

[14]:

# Sea ice thickness
aice_mean = np.ma.masked_invalid(ds_monthly['aice'].mean('month').values) * 100.
hi_mean = np.ma.masked_invalid(ds_monthly['hi'].mean('month').values)
hi_mean = np.ma.masked_where(aice_mean < 15.0, hi_mean)
suptitle = ('ANN mean sea ice thickness, ' + str(args.start_date) + ' to ' + str(args.end_date))
%matplotlib inline
# SH
fig = plt.figure(figsize=(10,8))
ax = plt.gca(projection=ccrs.SouthPolarStereo())
polarplot(hi_mean, grd, title=dcase.casename, debug=False, colormap=plt.cm.gist_ncar, clim=(0,2), axis=ax)
plt.suptitle(suptitle)

# NH
fig = plt.figure(figsize=(10,8))
ax = plt.gca(projection=ccrs.NorthPolarStereo())
polarplot(hi_mean, grd, title=dcase.casename, debug=False, colormap=plt.cm.gist_ncar, clim=(0,4), axis=ax, proj='NP')
plt.suptitle(suptitle)

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[14], line 9
      7 # SH
      8 fig = plt.figure(figsize=(10,8))
----> 9 ax = plt.gca(projection=ccrs.SouthPolarStereo())
     10 polarplot(hi_mean, grd, title=dcase.casename, debug=False, colormap=plt.cm.gist_ncar, clim=(0,2), axis=ax)
     11 plt.suptitle(suptitle)

TypeError: gca() got an unexpected keyword argument 'projection'

<Figure size 1000x800 with 0 Axes>

JFM

[ ]:

aice_JFM = np.ma.masked_invalid(ds_monthly['aice'].sel(month=[1,2,3]).mean('month').values) * 100.
hi_JFM = np.ma.masked_invalid(ds_monthly['hi'].sel(month=[1,2,3]).mean('month').values)
hi_JFM = np.ma.masked_where(aice_JFM < 15.0, hi_JFM)
suptitle = ('JFM mean sea ice thickness, ' + str(args.start_date) + ' to ' + str(args.end_date))
%matplotlib inline
# SH
fig = plt.figure(figsize=(10,8))
ax = plt.gca(projection=ccrs.SouthPolarStereo())
polarplot(hi_JFM, grd, title=dcase.casename, debug=False, colormap=plt.cm.gist_ncar, clim=(0,2), axis=ax)
plt.suptitle(suptitle)

# NH
fig = plt.figure(figsize=(10,8))
ax = plt.gca(projection=ccrs.NorthPolarStereo())
polarplot(hi_JFM, grd, title=dcase.casename, debug=False, colormap=plt.cm.gist_ncar, clim=(0,4), axis=ax, proj='NP')
plt.suptitle(suptitle)

AMJ

[ ]:

aice_AMJ = np.ma.masked_invalid(ds_monthly['aice'].sel(month=[4,5,6]).mean('month').values) * 100.
hi_AMJ = np.ma.masked_invalid(ds_monthly['hi'].sel(month=[4,5,6]).mean('month').values)
hi_AMJ = np.ma.masked_where(aice_AMJ < 15.0, hi_AMJ)
suptitle = ('AMJ mean sea ice thickness, ' + str(args.start_date) + ' to ' + str(args.end_date))
%matplotlib inline
# SH
fig = plt.figure(figsize=(10,8))
ax = plt.gca(projection=ccrs.SouthPolarStereo())
polarplot(hi_AMJ, grd, title=dcase.casename, debug=False, colormap=plt.cm.gist_ncar, clim=(0,2), axis=ax)
plt.suptitle(suptitle)

# NH
fig = plt.figure(figsize=(10,8))
ax = plt.gca(projection=ccrs.NorthPolarStereo())
polarplot(hi_AMJ, grd, title=dcase.casename, debug=False, colormap=plt.cm.gist_ncar, clim=(0,4), axis=ax, proj='NP')
plt.suptitle(suptitle)

JAS

[ ]:

aice_JAS = np.ma.masked_invalid(ds_monthly['aice'].sel(month=[7,8,9]).mean('month').values) * 100.
hi_JAS = np.ma.masked_invalid(ds_monthly['hi'].sel(month=[7,8,9]).mean('month').values)
hi_JAS = np.ma.masked_where(aice_JAS < 15.0, hi_JAS)
suptitle = ('JAS mean sea ice thickness, ' + str(args.start_date) + ' to ' + str(args.end_date))
%matplotlib inline
# SH
fig = plt.figure(figsize=(10,8))
ax = plt.gca(projection=ccrs.SouthPolarStereo())
polarplot(hi_JAS, grd, title=dcase.casename, debug=False, colormap=plt.cm.gist_ncar, clim=(0,2), axis=ax)
plt.suptitle(suptitle)

# NH
fig = plt.figure(figsize=(10,8))
ax = plt.gca(projection=ccrs.NorthPolarStereo())
polarplot(hi_JAS, grd, title=dcase.casename, debug=False, colormap=plt.cm.gist_ncar, clim=(0,4), axis=ax, proj='NP')
plt.suptitle(suptitle)

OND

[ ]:

aice_OND = np.ma.masked_invalid(ds_monthly['aice'].sel(month=[10,11,12]).mean('month').values) * 100.
hi_OND = np.ma.masked_invalid(ds_monthly['hi'].sel(month=[10,11,12]).mean('month').values)
hi_OND = np.ma.masked_where(aice_OND < 15.0, hi_OND)
suptitle = ('OND mean sea ice thickness, ' + str(args.start_date) + ' to ' + str(args.end_date))
%matplotlib inline
# SH
fig = plt.figure(figsize=(10,8))
ax = plt.gca(projection=ccrs.SouthPolarStereo())
polarplot(hi_OND, grd, title=dcase.casename, debug=False, colormap=plt.cm.gist_ncar, clim=(0,2), axis=ax)
plt.suptitle(suptitle)

# NH
fig = plt.figure(figsize=(10,8))
ax = plt.gca(projection=ccrs.NorthPolarStereo())
polarplot(hi_OND, grd, title=dcase.casename, debug=False, colormap=plt.cm.gist_ncar, clim=(0,4), axis=ax, proj='NP')
plt.suptitle(suptitle)

[ ]:

(obs_SH_mean.sic * 100.).plot()

[ ]:

tmp = obs_SH_mean['sic'].values * 100.
#tmp1 = np.ma.masked_where(tmp < 15.0, tmp)
fig = plt.figure(figsize=(10,8))
ax = plt.gca(projection=ccrs.SouthPolarStereo())
ax.contour(obs_SH_mean.lon.values, obs_SH_mean.lat.values, tmp, levels=[16.5], transform=ccrs.PlateCarree(), colors='black')

Sea ice concentration

[ ]:

aice_mean = np.ma.masked_where(aice_mean <= 15.0, aice_mean)

suptitle = ('ANN mean ice concentration, ' + str(args.start_date) + ' to ' + str(args.end_date))
%matplotlib inline
# SH
fig = plt.figure(figsize=(10,8))
ax = plt.gca(projection=ccrs.SouthPolarStereo())
polarplot(aice_mean, grd, title=dcase.casename, debug=False, colormap=plt.cm.gist_ncar, clim=(0,100), axis=ax)
ax.contour(obs_SH_mean.lon.values, obs_SH_mean.lat.values, tmp, levels=[16.5], transform=ccrs.PlateCarree(), colors='black',
           linestyles='solid')
plt.suptitle(suptitle)

# NH
fig = plt.figure(figsize=(10,8))
ax = plt.gca(projection=ccrs.NorthPolarStereo())
polarplot(aice_mean, grd, title=dcase.casename, debug=False, colormap=plt.cm.gist_ncar, clim=(0,100), axis=ax, proj='NP')
plt.suptitle(suptitle)

Northern Hemisphere

Maximum sea ice thickness

[ ]:

from mom6_tools.latlon_analysis import create_xarray_dataset, plot_area_ave_stats
from mom6_tools.m6plot import myStats

[ ]:

dtime = seaice.time.values
# variable name
var = 'hi'
# create datasets
ds_sh = create_xarray_dataset(var,'m',dtime)
ds_nh = create_xarray_dataset(var,'m',dtime)

# loop in time
for t in range(0,len(dtime)):
    # northern hemisphere
    tmp = np.ma.masked_invalid(seaice[var].sel(time=dtime[t]).values)
    tmp_nh = np.ma.masked_where(grd.geolat < 0, tmp)
    # get stats
    sMin, sMax, mean, std, rms = myStats(tmp_nh, area=None)
    # update Dataset
    ds_nh[var][1,t]  = sMax
    # southern hemisphere
    tmp_sh = np.ma.masked_where(grd.geolat > 0, tmp)
    # get stats
    sMin, sMax, mean, std, rms = myStats(tmp_sh, area=None)
    # update Dataset
    ds_sh[var][1,t] = sMax

Southern Hemisphere

[ ]:

sh = ds_sh.sel(stats='max')
fig = plt.figure(figsize=(12, 6))
ax  = fig.add_subplot(111)
sh.hi.plot(ax=ax)
plt.title('Maximum sea ice thickness, southern hemisphere', fontsize=14)
plt.xlabel('Time [years]')

Northern Hemisphere

[ ]:

nh = ds_nh.sel(stats='max')
fig = plt.figure(figsize=(12, 6))
ax  = fig.add_subplot(111)
nh.hi.plot(ax=ax)
plt.title('Maximum sea ice thickness, northern hemisphere', fontsize=14)
plt.xlabel('Time [years]')

Sea ice extent

[ ]:

dtime = seaice.time.values
# variable name
var = 'aice'
#
area_sh = []
area_nh = []

# loop in time
for t in range(len(dtime)):
    tmp = np.ma.masked_invalid(seaice[var].sel(time=dtime[t]).values)
    # mask southern hemisphere
    tmp_nh = np.ma.masked_where(grd.geolat < 0, grd.area_t)
    # mask northern hemisphere
    tmp_sh = np.ma.masked_where(grd.geolat > 0, grd.area_t)
    # mask if aggregrated concentration <= 0.15
    tmp_nh = np.ma.masked_where(tmp <= 0.15, tmp_nh*tmp)
    tmp_sh = np.ma.masked_where(tmp <= 0.15, tmp_sh*tmp)
    # append data
    area_sh.append(np.sum(tmp_sh)*1.0e-12) # 1.0e6 km^2
    area_nh.append(np.sum(tmp_nh)*1.0e-12) # 1.0e6 km^2

Southern hemisphere

[ ]:

fig = plt.figure(figsize=(12, 6))
plt.plot(dtime,area_sh)
plt.title('Sea ice area, southern hemisphere', fontsize=14)
plt.xlabel('Time [years]', fontsize=14); plt.ylabel('x 1.0e6 km$^2$', fontsize=14)

Northern hemisphere

[ ]:

fig = plt.figure(figsize=(12, 6))
plt.plot(dtime,area_nh)
plt.title('Sea ice area, northern hemisphere', fontsize=14)
plt.xlabel('Time [years]', fontsize=14); plt.ylabel('x 1.0e6 km$^2$', fontsize=14)

Under development

[ ]:

import matplotlib.path as mpath
fig = plt.figure(figsize=[10, 8])
ice = (seaice.aice[-1,:,:].data)
ax1 = plt.subplot(1, 1, 1, projection=ccrs.NorthPolarStereo())
ax1.set_extent([-180, 180, 50, 90], ccrs.PlateCarree())
ax1.add_feature(cartopy.feature.LAND)
ax1.gridlines()
# Compute a circle in axes coordinates, which we can use as a boundary
# for the map. We can pan/zoom as much as we like - the boundary will be
# permanently circular.
theta = np.linspace(0, 2*np.pi, 100)
center, radius = [0.5, 0.5], 0.5
verts = np.vstack([np.sin(theta), np.cos(theta)]).T
circle = mpath.Path(verts * radius + center)
ax1.set_boundary(circle, transform=ax1.transAxes)
#colormap = 'rainbow'
cs = ax1.pcolormesh(grd.geolon,grd.geolat,ice,transform=ccrs.PlateCarree(),cmap='Blues_r', shading='flat')
fig.colorbar(cs)
# Add Land
ax1.add_feature( cartopy.feature.LAND, zorder=1, edgecolor='none', facecolor='#fae5c9') #fae5c9')
# add Ocean
ax1.add_feature(cartopy.feature.OCEAN)
# Add coastline
ax1.coastlines(color='black')
# Add lat lon rings
ax1.gridlines(alpha='0.1',color='black')
im1 = ax1.contour(grd.geolon,grd.geolat,ice,[15],colors='red', transform=ccrs.PlateCarree())
plt.title('Model, Years ' + str(args.year_start) + ' to ' + str(args.year_end))