<jupyter_start><jupyter_text>Calculating Northern Hemisphere Sea Ice MetricsThis notebook describes how to calculate northern hemisphere sea ice extent and sea ice area from monthly grids of sea ice concentration.Northern hemisphere sea ice extent is the total area for each timestep of northern hemisphere grid cells with greater than 15% sea ice concentration.  Northern hemisphere sea ice area is the sum for each timestep of the products of fractional sea ice concentration and grid cell area.  I follow the ESGF and CMIP naming conventions and call northern hemipshere sea ice extent `siextentn` and nothern hemipshere sea ice area `siarean`.$$sixetentn = \sum{\alpha \times siconc \times areacello}; \alpha = \left\{                                                               \begin{array}{11}                                                              1; siconc \geq 0.15\\                                                              0; siconc < 0.15                                                              \end{array}                                                              \right.$$$$siarean = \sum{siconc \times areacello}$$Two files are required each model; a `siconc` file and an `areacello` file.  `areacello` is the area in $m^2$ for ocean grid cells.  The appropriate grid cell area variable for `siconc` is given in the `cell_measures` attribute for `siconc`.  This can be checked on a Linux machine by "grepping" for the `cell_measures` line in the netcdf header using `ncdump -h`<jupyter_code>%%bash
for f in `ls /home/apbarret/Data/CMIP6/siconc/SImon/*/historical/r1i1p1f1/siconc_SImon_*_historical_r1i1p1f1_gn_185001-*.nc`; 
do 
result=`ncdump -h $f | grep "cell_measures"`; 
echo "$(basename -- $f) $result"; 
done<jupyter_output>siconc_SImon_ACCESS-CM2_historical_r1i1p1f1_gn_185001-201412.nc 		siconc:cell_measures = "area: areacello" ;
siconc_SImon_ACCESS-ESM1-5_historical_r1i1p1f1_gn_185001-201412.nc 		siconc:cell_measures = "area: areacello" ;
siconc_SImon_AWI-ESM-1-1-LR_historical_r1i1p1f1_gn_185001-185012.nc 		siconc:cell_measures = "area: areacello" ;
siconc_SImon_BCC-CSM2-MR_historical_r1i1p1f1_gn_185001-201412.nc 		siconc:cell_measures = "area: areacello" ;
siconc_SImon_BCC-ESM1_historical_r1i1p1f1_gn_185001-201412.nc 		siconc:cell_measures = "area: areacello" ;
siconc_SImon_CAMS-CSM1-0_historical_r1i1p1f1_gn_185001-201412.nc 		siconc:cell_measures = "area: areacello" ;
siconc_SImon_CanESM5_historical_r1i1p1f1_gn_185001-201412.nc 		siconc:cell_measures = "area: areacello" ;
siconc_SImon_CAS-ESM2-0_historical_r1i1p1f1_gn_185001-201412.nc 		siconc:cell_measures = "area: areacello" ;
siconc_SImon_CESM2-FV2_historical_r1i1p1f1_gn_185001-189912.nc 		siconc:cell_measures = "area: areacello" ;
siconc_SImon_CES[...]<jupyter_text>Frustratingly, models do not have a consistent naming convention for dimensions.  Some modelling groups have used `x` and `y` for grid dimensions and some have used `i` and `j`.  Latitude and longitude coordinates can be `lat` and `lon` or `latitude` and `longitude`.  In some cases, dimension names are not consistent between `siconc` and `areacello` files.  So dimension names in one of the files have to be renamed to so that dimensions names match.<jupyter_code>import os
import glob
import datetime as dt

import matplotlib.pyplot as plt

import xarray as xr
import pandas as pd
DATA_PATH = "/home/apbarret/Data/CMIP6"
SICONC_PATH = os.path.join(DATA_PATH, "siconc", "SImon")
AREACELLO_PATH = os.path.join(DATA_PATH, "areacello")<jupyter_output><empty_output><jupyter_text>Matching models with `areacello`<jupyter_code>model_list = [os.path.basename(p) for p in glob.glob(os.path.join(SICONC_PATH, '*'))]
model_list
model = "CESM2"
experiment = "historical"
ensemble = "r1i1p1f1"

EXCLUDE_THESE_MODELS = ['CESM2-WACCM']<jupyter_output><empty_output><jupyter_text>Get the data<jupyter_code>def find_siconc_files(model, experiment, ensemble):
    files = glob.glob(os.path.join(SICONC_PATH, model, experiment, ensemble, "*.nc"))
    if not files:
        raise Exception(f"No siconc files found for {model} {experiment} {ensemble}")
    return files


def get_siconc(filepaths):
    '''Returns xarray.Dataset for a given model run'''
    ds = xr.open_mfdataset(filepaths, combine="by_coords", use_cftime=True)
    #ds['time'] = ds.indexes['time'].to_datetimeindex()
    if 'lat' in ds.variables:
        ds = ds.rename({'lat': 'latitude', 'lon': 'longitude'})
    if 'ni' in ds.dims:
        ds = ds.rename_dims({'ni': 'i', 'nj': 'j'})
    if 'ni' in ds.variables:
        ds = ds.rename({'ni': 'i', 'nj': 'j'})
    if 'nav_lat' in ds.variables:
        ds = ds.rename({'nav_lat': 'latitude', 'nav_lon': 'longitude'})
    if 'x' in ds.variables:
        ds = ds.rename({'x': 'i', 'y': 'j'})
    if 'x' in ds.dims:
        ds = ds.rename({'x': 'i', 'y': 'j'})
    ds['siconc'] = ds.siconc * 1e-2
    return ds


def find_areacello_file(model):
    files = glob.glob(os.path.join(AREACELLO_PATH, f"areacello_*_{model}_*.nc"))
    if not files:
        raise Exception(f"No areacello file for {model}")
    return files[0]


def get_areacello(filepath, verbose=False):
    '''Returns xarray.Dataset containing ocean gridcell area'''
    if verbose: print(f"Getting ocean cell area from {area_file}")
    ds = xr.open_dataset(filepath, use_cftime=True)
    if 'nlat' in ds.variables:
        ds = ds.rename({'nlat': 'j', 'nlon': 'i'})  # rename nlat and nlon dimensions to match siconc
    if 'nav_lat' in ds.variables:
        ds = ds.rename({'nav_lat': 'latitude', 'nav_lon': 'longitude'})
    if 'lat' in ds.variables:
        ds = ds.rename({'lat': 'latitude', 'lon': 'longitude'})
    if 'x' in ds.variables:
        ds = ds.rename({'x': 'i', 'y': 'j'})
    if 'x' in ds.dims:
        ds = ds.rename({'x': 'i', 'y': 'j'})
    if any([substr in filepath for substr in ['FGOALS-f3-L', 'FGOALS-g3']]):
        ds['j'] = list(reversed(ds.indexes['j']))
    ds['areacello'] = ds.areacello * 1e-6  # convert from m^2 to km^2
    return ds


def calc_siarean(siconc, area):
    siarea = siconc.siconc.where(siconc.siconc > 0.) * area.areacello
    return siarea.where(area.latitude > 0.).sum(dim=['j', 'i'])


def calc_siextentn(siconc, area):
    siextent = siconc.siconc.where(siconc.siconc > .15) * area.areacello
    return siextent.where(area.latitude > 0.).sum(dim=['j', 'i'])


def load_data(model):
    '''Loads siconc and area data as xarray.Datasets'''
    try:
        siconc_files = find_siconc_files(model, experiment, ensemble)
    except Exception as error:
        raise Exception(error)
    
    try:
        areacello_file = find_areacello_file(model)
    except Exception as error:
        raise Exception(error)
        
    siconc_ds = get_siconc(siconc_files)
    area_ds = get_areacello(areacello_file)

    assert set(area_ds.areacello.dims).issubset(set(siconc_ds.siconc.dims)), f"{area_ds.areacello.dims} is not subset of {siconc_ds.siconc.dims} for {model}"
    assert 'latitude' in area_ds.variables and 'longitude' in area_ds.variables, f"latitude and longitude are not variables in areacello Dataset of {model}"
    assert 'latitude' in siconc_ds.variables and 'longitude' in siconc_ds.variables, f"latitude and longitude are not variables in siconc Dataset of {model}"
    assert 'i' in area_ds.areacello.dims and 'j' in area_ds.areacello.dims, f"Expect dimensions ('i', 'j') for areacello for {model}, found {area_ds.areacello.dims} instead"
    assert 'i' in siconc_ds.siconc.dims and 'j' in siconc_ds.siconc.dims, f"Expect dimensions ('i', 'j') for siconc for {model}, found {siconc_ds.areacello.dims} instead"

    return siconc_ds, area_ds
siarean_dct = {}
siextentn_dct = {}
for model in model_list:
    
    print(model)
    
    if model in EXCLUDE_THESE_MODELS:
        continue
    
    try:
        siconc_ds, area_ds = load_data(model)
    except Exception as error:
        print(error)
        print("Skipping...\n")
        continue

    #siarean = calc_siarean(siconc_ds, area_ds)
    #siarean['time'] = [dt.datetime(t.year, t.month, 1) for t in siarean.time.values]
    #siarean_dct[model] = siarean.to_series()
    
    siextentn = calc_siextentn(siconc_ds, area_ds)
    siextentn['time'] = [dt.datetime(t.year, t.month, 1) for t in siextentn.time.values]
    siextentn_dct[model] = siextentn.to_series()
    
    siconc_ds.close()
    area_ds.close()

#siarean_df = pd.DataFrame(siarean_dct)
siextentn_df = pd.DataFrame(siextentn_dct)
siextentn_df
siextentn_df.to_csv('cmip6_siextentn_historical_r1i1p1f1.from_siconc.csv')
siextentn_df[siextentn_df.index.month == 9].plot()
id=22
model = 'CIESM'
fig, ax = plt.subplots(figsize=(15,5))
siarean_df[siarean_df.index.month == 9].loc[:,model].plot(title=model, ylim=(0,10e6))
model='CESM2'
siconc_file = find_siconc_files(model, experiment, ensemble)
areacello_file = find_areacello_file(model)

#data = xr.open_mfdataset(siconc_file, combine='by_coords')
#area = xr.open_dataset(areacello_file)
data = get_siconc(siconc_file)
area = get_areacello(areacello_file)

siarean = calc_siarean(data, area)
siarean['time'] = [dt.datetime(t.year, t.month, 1) for t in siarean.time.values]


ts = siarean.to_series()

ts[ts.index.month == 9].plot()
area.areacello.max()
area
area['j'] = list(reversed(area.indexes['j']))
model='CESM2'
siconc_file = find_siconc_files(model, experiment, ensemble)
print(siconc_file)
cesm2 = get_siconc(siconc_file)
cesm2.siconc.max().values
model='CIESM'
siconc_file = find_siconc_files(model, experiment, ensemble)
print(siconc_file)
ciesm = get_siconc(siconc_file)
ciesm.siconc.max().values
diff = cesm2.siconc.isel(time=1779) - ciesm.siconc.isel(time=1779)
diff
diff.plot()
ds
ds.siconc.max().values<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12370286653652522, 'fraction_numerical': 0.02465837871159971, 'mean_word_length': 3.901309164149043, 'pattern_counts': {'":': 0, '<': 17, '<?xml version=': 0, '>': 20, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '662023', 'license': 'MIT', 'n_tokens_mistral': 3858, 'n_tokens_neox': 3669, 'n_words': 825, 'path': 'notebooks/northern_hemisphere_sea_ice_from siconc.ipynb', 'repo_name': 'andypbarrett/cmip6', 'size': '9733'}