CryoSat-2 data wrangling

CryoSat-2 data wrangling#

Summary This notebook loads the different CryoSat-2 datasets used in our ICESat-2/CryoSat-2 comparison analysis into a single NETCDF4 file, with descriptive attributes maintained for each dataset. Each dataset is regridded to the ICESat2 grid shape [304, 448] (x,y). The datasets used in this notebook are listed below.

Version history: Version 1 (01/01/2022)

Details on the CryoSat-2 datasets#

NASA GSFC CryoSat-2 monthly mean winter Arctic sea ice thickness#

Download link:
- https://nsidc.org/data/rdeft4/
Reference: Kurtz, N. and J. Harbeck. (2017). CryoSat-2 Level-4 Sea Ice Elevation, Freeboard, and Thickness, Version 1 [Data Set]. Boulder, Colorado USA. NASA National Snow and Ice Data Center Distributed Active Archive Center. https://doi.org/10.5067/96JO0KIFDAS8

CPOM CryoSat-2 monthly mean winter Arctic sea ice thickness#

Download link:
- http://www.cpom.ucl.ac.uk/csopr/seaice.php
Reference: Laxon, S. W. et al. CryoSat-2 estimates of Arctic sea ice thickness and volume. Geophysical Research Letters 40, 732-737 (2013).
Notes: data are posted on a higher-res grid than the other datasets.

AWI CS-2/SMOS CryoSat-2 monthly mean winter Arctic sea ice thickness#

Download link:
- https://spaces.awi.de/pages/viewpage.action?pageId=291898639
Reference: Ricker, R., Hendricks, S., Kaleschke, L., Tian-Kunze, X., King, J., and Haas, C.: A weekly Arctic sea-ice thickness data record from merged CryoSat-2 and SMOS satellite data, The Cryosphere, 11, 1607-1623, https://doi.org/10.5194/tc-11-1607-2017, 2017.
Notes: using the merged SMOS product here as it seems to be the prefered choice of this group.

University of Bristol CryoSat-2 monthly mean all-season Arctic sea ice thickness#

Download link:
- https://data.bas.ac.uk/full-record.php?id=GB/NERC/BAS/PDC/01613
Reference: Landy, J.C., Dawson, G.J., Tsamados, M. et al. A year-round satellite sea-ice thickness record from CryoSat-2. Nature 609, 517–522 (2022). https://doi.org/10.1038/s41586-022-05058-5
Notes: Just using the data made available through the publication, more data has been generated since. Only dataset producing estimates through summer!

JPL ICESat-2/CryoSat-2 monthly mean winter Arctic sea ice thickness#

Download link:
- https://icesat-2.gsfc.nasa.gov/sea-ice-data/kacimi-kwok-2022
Reference: Kacimi, S., Kwok, R. (2022), Arctic snow depth, ice thickness and volume from ICESat-2 and CryoSat-2: 2018-2021, Geophysical Research Letters, doi: 10.1029/2021GL097448.
Notes: no projection information but was told it’s on the NSIDC EPSG 3411 projection

Note

Although you’ll see an option to run this notebook in Binder, this notebook is NOT configured to run in Binder. If you want to wrangle the data yourself, each dataset used to compile the final data product can be downloaded from the links above. The final data product produced by this notebook can be downloaded from the google storage bucket associated with this jupyter book.

Import notebook dependencies#

import xarray as xr # For working with gridded climate data 
import pandas as pd
import numpy as np
import itertools
import pyproj 
from netCDF4 import Dataset
import scipy.interpolate 
from utils.read_data_utils import read_book_data # Helper function for reading the data from the bucket
from utils.plotting_utils import compute_gridcell_winter_means, interactiveArcticMaps, interactive_winter_mean_maps, interactive_winter_comparison_lineplot # Plotting

# Plotting dependencies
import cartopy.crs as ccrs
from textwrap import wrap
import hvplot.xarray
import holoviews as hv
import matplotlib.pyplot as plt
from matplotlib.axes import Axes
from cartopy.mpl.geoaxes import GeoAxes
GeoAxes._pcolormesh_patched = Axes.pcolormesh # Helps avoid some weird issues with the polar projection 
%config InlineBackend.figure_format = 'retina'
import matplotlib as mpl
mpl.rcParams['figure.dpi'] = 150 # Sets figure size in the notebook

# Remove warnings to improve display
import warnings 
warnings.filterwarnings('ignore') 

def getCS2gsfc(dataPathCS2, yearStr, mStr):
    """ Read in GSFC CryoSat-2 sea ice thickness data
        
        Just use the date from the middle of the month (15th for simplicity) 
        as I believe that should be essentially the same as our monthly mean.
        
        Downloaded from the NSIDC: https://n5eil01u.ecs.nsidc.org/ICEBRIDGE/RDEFT4.001/
    
    Args:
        yearStr (str): year
        mStr (str, e.g. '01' or '12'): month
        
    Returns
        xptsT (2d numpy array): x coordinates on our map projection
        yptsT (2d numpy array): y coordinates on our map projection
        thicknessCS (2d numpy array): monthly sea ice thickness estimates
        

    """

    #f = xr.open_dataset(dataPathCS2+'/GSFC/'+yearStr+'/RDEFT4_'+yearStr+mStr+'15.nc')
    f = Dataset(dataPathCS2+'/GSFC/'+yearStr+'/RDEFT4_'+yearStr+mStr+'15.nc', 'r')
    thicknessCS = f.variables['sea_ice_thickness'][:]
    #thicknessCS=ma.masked_where(thicknessCS<0, thicknessCS)
    #thicknessCS=ma.masked_where(np.isnan(thicknessCS), thicknessCS)
    
    thicknessCS[np.where(thicknessCS<0.)]=np.nan
    

    latsCS = f.variables['lat'][:]
    lonsCS = f.variables['lon'][:]

    xptsT, yptsT = mapProj(lonsCS, latsCS)

    return xptsT, yptsT, thicknessCS

def getCS2cpom(dataPathCS2, yearStr, mStr, res=1):
    """ Read in the CPOM-UCL CryoSat-2 sea ice thickness data
        
        Data is quite high-res (5 km) so can be easier to just coarsen before using. 
        Think it's fine to just take every Nth x/y point as a 25 km smoother is applied, which is really the effective resolution of each point.
        
        Downloaded from the CPOM portal: http://www.cpom.ucl.ac.uk/csopr/seaice.php
    
    Args:
        dataPathCS2 (str): location of data
        yearStr (str): year
        mStr (str, e.g. '1' or '12'): month (NB THIS USED TO BE DOULBE STRINGED MONTHS)
        res (int): pick every res data point in x/y.
        
    Returns
        xptsT (2d numpy array): x coordinates on our map projection
        yptsT (2d numpy array): y coordinates on our map projection
        thicknessCS (2d numpy array): monthly sea ice thickness estimates
        

    """
    
    #print(dataPathCS2+'/CPOM/thk_'+yearStr+'_'+mStr+'.map.nc')
    f = Dataset(dataPathCS2+'/CPOM/thk_'+yearStr+'_'+mStr+'.map.nc', 'r')
    latsCS = f.variables['latitude'][::res]
    lonsCS = f.variables['longitude'][::res]
    
    thicknessCS = f.variables['thickness'][::res]

    xptsT, yptsT = mapProj(lonsCS, latsCS)

    #files = glob(dataPath+ystr+mstr+'*')
    return xptsT, yptsT, thicknessCS

def getCS2awismos(dataPathCS2, yearStr, mStr):
    """ Read in the AWI CryoSat-2/SMOW merged sea ice thickness data

        Downloaded from the AWI portal: https://spaces.awi.de/pages/viewpage.action?pageId=291898639
    
    Args:
        dataPathCS2 (str): location of data
        yearStr (str): year
        mStr (str, e.g. '01' or '12'): month
        
    Returns
        xptsT (2d numpy array): x coordinates on our map projection
        yptsT (2d numpy array): y coordinates on our map projection
        thicknessCS (2d numpy array): monthly sea ice thickness estimates
        

    """
    #print(dataPathCS2+'/AWI_SMOS/'+yearStr+'/'+mStr+'/*SMOS*'+yearStr+mStr+'*'+yearStr+mStr+'*.nc')
    f = xr.open_mfdataset(dataPathCS2+'/AWI_SMOS/'+yearStr+'/'+mStr+'/*SMOS*'+yearStr+mStr+'*'+yearStr+mStr+'*.nc').mean(dim="time")

    thicknessCS = f['weighted_mean_sea_ice_thickness'].values

    xptsT, yptsT = mapProj(f.lon, f.lat)

    #files = glob(dataPath+ystr+mstr+'*')
    return xptsT, yptsT, thicknessCS

def getCS2kk(dataPathCS2, yearStr, mStr):
    """ Read in the Kacimi and Kwok CryoSat-2 sea ice thickness data
    
    Assuming it's on EPSG 3411, not given.
    
    Args:
        dataPathCS2 (str): location of data
        yearStr (str): year
        mStr (str, e.g. '1' or '12'): month (NB THIS USED TO BE DOULBE STRINGED MONTHS)
        
    Returns
        xptsT (2d numpy array): x coordinates on our map projection
        yptsT (2d numpy array): y coordinates on our map projection
        thicknessCS (2d numpy array): monthly sea ice thickness estimates
        

    """
    
    print(dataPathCS2+'/KacimiKwok/thk_'+yearStr+mStr+'.txt')
    k = pd.read_csv(dataPathCS2+'/KacimiKwok/thk_'+yearStr+mStr+'.txt')

    thicknessCS = k.data.values
    xptsT = k.X.values*1000.
    yptsT = k.Y.values*1000.

    #files = glob(dataPath+ystr+mstr+'*')
    return xptsT, yptsT, thicknessCS

def getCS2ubris(dataPathCS2):
    """ Read in the University of Bristol CryoSat-2 sea ice thickness data

    
    Args:
        dataPathCS2 (str): location of data
        
    Returns
        xptsT (2d numpy array): x coordinates on our map projection
        yptsT (2d numpy array): y coordinates on our map projection
        thicknessCS (2d numpy array): monthly sea ice thickness estimates
        

    """
    ubris_f = xr.open_dataset(dataPathCS2+'/UBRIS/ubristol_cryosat2_seaicethickness_nh_80km_v1p7.nc', decode_times=False)

    # Issue with time starting from year 0!
    # Re-set it to start from some other year
    ubris_f = ubris_f.rename({'Time':'time'})
    ubris_f['time'] = ubris_f['time']-679352
    ubris_f.time.attrs["units"] = "days since 1860-01-01"
    decoded_time = xr.decode_cf(ubris_f)

    ubris_f['time']=decoded_time.time
    ubris_f = ubris_f.swap_dims({'t': 'time'})

    # Resample to monthly, note that the S just makes the index start on the 1st of the month
    thicknessCS = ubris_f.resample(time="MS").mean()
    xptsT, yptsT = mapProj(thicknessCS.isel(time=0).Longitude, thicknessCS.isel(time=0).Latitude)
    
    return xptsT, yptsT, thicknessCS

def regridToICESat2(dataArrayNEW, xptsNEW, yptsNEW, xptsIS2, yptsIS2):  
    """ Regrid new data to ICESat-2 grid 
    
    Args: 
        dataArrayNEW (xarray DataArray): Numpy variable array to be gridded to ICESat-2 grid 
        xptsNEW (numpy array): x-values of dataArrayNEW projected to ICESat-2 map projection 
        yptsNEW (numpy array): y-values of dataArrayNEW projected to ICESat-2 map projection 
        xptsIS2 (numpy array): ICESat-2 longitude projected to ICESat-2 map projection
        yptsIS2 (numpy array): ICESat-2 latitude projected to ICESat-2 map projection
    
    Returns: 
        gridded (numpy array): data regridded to ICESat-2 map projection
    
    """
    #gridded = []
    #for i in range(len(dataArrayNEW.values)): 
    try:
        gridded = scipy.interpolate.griddata((xptsNEW.flatten(),yptsNEW.flatten()), dataArrayNEW.flatten(), (xptsIS2, yptsIS2), method = 'nearest')
    except:
        try:
            gridded = scipy.interpolate.griddata((xptsNEW,yptsNEW), dataArrayNEW, (xptsIS2, yptsIS2), method = 'nearest')
        except:
            print('Error interpolating..')
    
    return gridded

# Read/download the data 
book_ds = read_book_data()
book_ds['region_mask'] = book_ds['region_mask'].isel(time=0, drop=True)

dataPathCS2='/Users/aapetty/Data/CS2'

# Initialize map projection and project data to it
out_proj = 'EPSG:3411'
out_lons = book_ds.longitude.values
out_lats = book_ds.latitude.values

mapProj = pyproj.Proj("+init=" + out_proj)
xptsIS2, yptsIS2 = mapProj(out_lons, out_lats)

start_date = "Nov 2018"
end_date = "Apr 2021"
date_range = pd.date_range(start=start_date, end=end_date, freq='MS') # MS indicates a time frequency of start of the month
#date_range = date_range[((date_range.month <5) | (date_range.month > 8))]
date_range

DatetimeIndex(['2018-11-01', '2018-12-01', '2019-01-01', '2019-02-01',
               '2019-03-01', '2019-04-01', '2019-05-01', '2019-06-01',
               '2019-07-01', '2019-08-01', '2019-09-01', '2019-10-01',
               '2019-11-01', '2019-12-01', '2020-01-01', '2020-02-01',
               '2020-03-01', '2020-04-01', '2020-05-01', '2020-06-01',
               '2020-07-01', '2020-08-01', '2020-09-01', '2020-10-01',
               '2020-11-01', '2020-12-01', '2021-01-01', '2021-02-01',
               '2021-03-01', '2021-04-01'],
              dtype='datetime64[ns]', freq='MS')

cs2_gsfc = []
valid_dates=[]
for date in date_range:
    #print(date)
    try:
        xpts_cs2_gsfc, ypts_cs2_gsfc, cs2_gsfc_temp = getCS2gsfc(dataPathCS2, str(date.year), '%02d' %(date.month))
    except:
        print(date)
        print('no GSFC CS-2 data so skipping...')
        continue
    valid_dates.append(date)
    
    cs2_gsfc_temp_is2grid = regridToICESat2(cs2_gsfc_temp, xpts_cs2_gsfc, ypts_cs2_gsfc, xptsIS2, yptsIS2) 
    cs2_gsfc_temp_is2grid_xr = xr.DataArray(data = cs2_gsfc_temp_is2grid, 
                            dims = ['y', 'x'], 
                            coords = {'latitude': (('y','x'), out_lats), 'longitude': (('y','x'), out_lons)}, 
                            name = 'CS-2_GSFC')

    cs2_gsfc.append(cs2_gsfc_temp_is2grid_xr)
                         
cs2_gsfc = xr.concat(cs2_gsfc, 'time')
cs2_gsfc = cs2_gsfc.assign_coords(time=valid_dates)
cs2_gsfc_attrs = {'units': 'meters', 'long_name': 'GSFC CryoSat-2 monthly mean Arctic sea ice thickness', 'data_download': 'https://nsidc.org/data/rdeft4/', 
        'download_date': '09-2022', 'citation': 'Kurtz, N. and J. Harbeck. (2017). CryoSat-2 Level-4 Sea Ice Elevation, Freeboard, and Thickness, Version 1 [Data Set]. Boulder, Colorado USA. NASA National Snow and Ice Data Center Distributed Active Archive Center. https://doi.org/10.5067/96JO0KIFDAS8'}

cs2_gsfc = cs2_gsfc.assign_attrs(cs2_gsfc_attrs)  
# Add to book
book_ds['cs2_sea_ice_thickness_GSFC']=cs2_gsfc

2018-11-01 00:00:00
no GSFC CS-2 data so skipping...
2018-12-01 00:00:00
no GSFC CS-2 data so skipping...
2019-01-01 00:00:00
no GSFC CS-2 data so skipping...
2019-02-01 00:00:00
no GSFC CS-2 data so skipping...
2019-03-01 00:00:00
no GSFC CS-2 data so skipping...
2019-04-01 00:00:00
no GSFC CS-2 data so skipping...
2019-05-01 00:00:00
no GSFC CS-2 data so skipping...
2019-06-01 00:00:00
no GSFC CS-2 data so skipping...
2019-07-01 00:00:00
no GSFC CS-2 data so skipping...
2019-08-01 00:00:00
no GSFC CS-2 data so skipping...
2019-09-01 00:00:00
no GSFC CS-2 data so skipping...
2019-10-01 00:00:00
no GSFC CS-2 data so skipping...
2019-11-01 00:00:00
no GSFC CS-2 data so skipping...
2019-12-01 00:00:00
no GSFC CS-2 data so skipping...
2020-01-01 00:00:00
no GSFC CS-2 data so skipping...
2020-02-01 00:00:00
no GSFC CS-2 data so skipping...
2020-03-01 00:00:00
no GSFC CS-2 data so skipping...
2020-04-01 00:00:00
no GSFC CS-2 data so skipping...
2020-05-01 00:00:00
no GSFC CS-2 data so skipping...
2020-06-01 00:00:00
no GSFC CS-2 data so skipping...
2020-07-01 00:00:00
no GSFC CS-2 data so skipping...
2020-08-01 00:00:00
no GSFC CS-2 data so skipping...
2020-09-01 00:00:00
no GSFC CS-2 data so skipping...
2020-10-01 00:00:00
no GSFC CS-2 data so skipping...
2020-11-01 00:00:00
no GSFC CS-2 data so skipping...
2020-12-01 00:00:00
no GSFC CS-2 data so skipping...
2021-01-01 00:00:00
no GSFC CS-2 data so skipping...
2021-02-01 00:00:00
no GSFC CS-2 data so skipping...
2021-03-01 00:00:00
no GSFC CS-2 data so skipping...
2021-04-01 00:00:00
no GSFC CS-2 data so skipping...

---------------------------------------------------------------------------
StopIteration                             Traceback (most recent call last)
File ~/miniconda3/envs/is2book_p39_env/lib/python3.9/site-packages/xarray/core/concat.py:226, in concat(objs, dim, data_vars, coords, compat, positions, fill_value, join, combine_attrs)
    225 try:
--> 226     first_obj, objs = utils.peek_at(objs)
    227 except StopIteration:

File ~/miniconda3/envs/is2book_p39_env/lib/python3.9/site-packages/xarray/core/utils.py:186, in peek_at(iterable)
    185 gen = iter(iterable)
--> 186 peek = next(gen)
    187 return peek, itertools.chain([peek], gen)

StopIteration: 

During handling of the above exception, another exception occurred:

ValueError                                Traceback (most recent call last)
Cell In[12], line 21
     14     cs2_gsfc_temp_is2grid_xr = xr.DataArray(data = cs2_gsfc_temp_is2grid, 
     15                             dims = ['y', 'x'], 
     16                             coords = {'latitude': (('y','x'), out_lats), 'longitude': (('y','x'), out_lons)}, 
     17                             name = 'CS-2_GSFC')
     19     cs2_gsfc.append(cs2_gsfc_temp_is2grid_xr)
---> 21 cs2_gsfc = xr.concat(cs2_gsfc, 'time')
     22 cs2_gsfc = cs2_gsfc.assign_coords(time=valid_dates)
     23 cs2_gsfc_attrs = {'units': 'meters', 'long_name': 'GSFC CryoSat-2 monthly mean Arctic sea ice thickness', 'data_download': 'https://nsidc.org/data/rdeft4/', 
     24         'download_date': '09-2022', 'citation': 'Kurtz, N. and J. Harbeck. (2017). CryoSat-2 Level-4 Sea Ice Elevation, Freeboard, and Thickness, Version 1 [Data Set]. Boulder, Colorado USA. NASA National Snow and Ice Data Center Distributed Active Archive Center. https://doi.org/10.5067/96JO0KIFDAS8'}

File ~/miniconda3/envs/is2book_p39_env/lib/python3.9/site-packages/xarray/core/concat.py:228, in concat(objs, dim, data_vars, coords, compat, positions, fill_value, join, combine_attrs)
    226     first_obj, objs = utils.peek_at(objs)
    227 except StopIteration:
--> 228     raise ValueError("must supply at least one object to concatenate")
    230 if compat not in _VALID_COMPAT:
    231     raise ValueError(
    232         f"compat={compat!r} invalid: must be 'broadcast_equals', 'equals', 'identical', 'no_conflicts' or 'override'"
    233     )

ValueError: must supply at least one object to concatenate

cs2_cpom = []
valid_dates=[]
for date in date_range:
    
    try:
        xpts_cs2_cpom, ypts_cs2_cpom, cs2_cpom_temp = getCS2cpom(dataPathCS2, str(date.year), '%01d' %(date.month), res=1)
    except:
        print(date)
        print('no CS-2 data so skipping...')
        continue
    valid_dates.append(date)
    
    cs2_cpom_temp_is2grid = regridToICESat2(cs2_cpom_temp, xpts_cs2_cpom, ypts_cs2_cpom, xptsIS2, yptsIS2) 
    cs2_cpom_temp_is2grid_xr = xr.DataArray(data = cs2_cpom_temp_is2grid, 
                            dims = ['y', 'x'], 
                            coords = {'latitude': (('y','x'), out_lats), 'longitude': (('y','x'), out_lons)}, 
                            name = 'CS-2_CPOM')

    cs2_cpom.append(cs2_cpom_temp_is2grid_xr)
                         
cs2_cpom = xr.concat(cs2_cpom, 'time')
cs2_cpom = cs2_cpom.assign_coords(time=valid_dates)
cs2_cpom_attrs = {'units': 'meters', 'long_name': 'CPOM CryoSat-2 monthly mean Arctic sea ice thickness', 'data_download': 'http://www.cpom.ucl.ac.uk/csopr/seaice.php', 
        'download_date': '09-2022', 'citation': 'Laxon, S. W. et al. CryoSat-2 estimates of Arctic sea ice thickness and volume. Geophysical Research Letters 40, 732-737 (2013).'}

cs2_cpom = cs2_cpom.assign_attrs(cs2_cpom_attrs)  
# Add to book
book_ds['cs2_sea_ice_thickness_CPOM']=cs2_cpom

2019-05-01 00:00:00
no CS-2 data so skipping...
2019-06-01 00:00:00
no CS-2 data so skipping...
2019-07-01 00:00:00
no CS-2 data so skipping...
2019-08-01 00:00:00
no CS-2 data so skipping...
2019-09-01 00:00:00
no CS-2 data so skipping...
2020-05-01 00:00:00
no CS-2 data so skipping...
2020-06-01 00:00:00
no CS-2 data so skipping...
2020-07-01 00:00:00
no CS-2 data so skipping...
2020-08-01 00:00:00
no CS-2 data so skipping...
2020-09-01 00:00:00
no CS-2 data so skipping...

cs2_awi = []
valid_dates=[]
for date in date_range:
    #print(date)
    try:
        xpts_cs2_awi, ypts_cs2_awi, cs2_awi_temp = getCS2awismos(dataPathCS2, str(date.year), '%02d' %(date.month))
    except:
        print('no CS-2 data so skipping...')
        continue
    valid_dates.append(date)
    
    cs2_awi_temp_is2grid = regridToICESat2(cs2_awi_temp, xpts_cs2_awi, ypts_cs2_awi, xptsIS2, yptsIS2) 
    cs2_awi_temp_is2grid_xr = xr.DataArray(data = cs2_awi_temp_is2grid, 
                            dims = ['y', 'x'], 
                            coords = {'latitude': (('y','x'), out_lats), 'longitude': (('y','x'), out_lons)}, 
                            name = 'CS-2_AWISMOS')

    cs2_awi.append(cs2_awi_temp_is2grid_xr)
                         
cs2_awi = xr.concat(cs2_awi, 'time')
cs2_awi = cs2_awi.assign_coords(time=valid_dates)
cs2_awi_attrs = {'units': 'meters', 'long_name': 'AWI SMOS & CryoSat-2 monthly mean Arctic sea ice thickness', 'data_download': 'https://spaces.awi.de/pages/viewpage.action?pageId=291898639', 
        'download_date': '09-2022', 'citation': 'Ricker, R., Hendricks, S., Kaleschke, L., Tian-Kunze, X., King, J., and Haas, C.: A weekly Arctic sea-ice thickness data record from merged CryoSat-2 and SMOS satellite data, The Cryosphere, 11, 1607-1623, https://doi.org/10.5194/tc-11-1607-2017, 2017.'}

cs2_awi = cs2_awi.assign_attrs(cs2_awi_attrs)  
# Add to book
book_ds['cs2_sea_ice_thickness_AWISMOS']=cs2_awi

no CS-2 data so skipping...
no CS-2 data so skipping...
no CS-2 data so skipping...
no CS-2 data so skipping...
no CS-2 data so skipping...
no CS-2 data so skipping...
no CS-2 data so skipping...
no CS-2 data so skipping...
no CS-2 data so skipping...
no CS-2 data so skipping...

cs2_ubris = []
valid_dates=[]

xptsT_ubris, yptsT_ubris, cs2_ubris_raw = getCS2ubris(dataPathCS2)

for date in date_range:
     
    try:
        cs2_ubris_temp_is2grid = regridToICESat2(cs2_ubris_raw.Sea_Ice_Thickness.sel(time=date).values, xptsT_ubris, yptsT_ubris, xptsIS2, yptsIS2) 
         
    except:
        print(date)
        print('no CS-2 data so skipping...')
        continue
    valid_dates.append(date)
    
    
    cs2_ubris_temp_is2grid_xr = xr.DataArray(data = cs2_ubris_temp_is2grid, 
                            dims = ['y', 'x'], 
                            coords = {'latitude': (('y','x'), out_lats), 'longitude': (('y','x'), out_lons)}, 
                            name = 'CS-2_UBRIS')
    
    cs2_ubris.append(cs2_ubris_temp_is2grid_xr)
                         
cs2_ubris = xr.concat(cs2_ubris, 'time')
#cs2_ubris = cs2_ubris.assign_coords(time=valid_dates)
cs2_ubris_attrs = {'units': 'meters', 'long_name': 'University of Bristol CryoSat-2 Arctic sea ice thickness', 'data_download': 'https://data.bas.ac.uk/full-record.php?id=GB/NERC/BAS/PDC/01613', 
        'download_date': '09-2022', 'citation': 'Landy, J.C., Dawson, G.J., Tsamados, M. et al. A year-round satellite sea-ice thickness record from CryoSat-2. Nature 609, 517–522 (2022). https://doi.org/10.1038/s41586-022-05058-5'} 
cs2_ubris = cs2_ubris.assign_coords(time=valid_dates)
cs2_ubris = cs2_ubris.assign_attrs(cs2_ubris_attrs)  
# Add to book
book_ds['cs2_sea_ice_thickness_UBRIS']=cs2_ubris

2020-08-01 00:00:00
no CS-2 data so skipping...
2020-09-01 00:00:00
no CS-2 data so skipping...
2020-10-01 00:00:00
no CS-2 data so skipping...
2020-11-01 00:00:00
no CS-2 data so skipping...
2020-12-01 00:00:00
no CS-2 data so skipping...
2021-01-01 00:00:00
no CS-2 data so skipping...
2021-02-01 00:00:00
no CS-2 data so skipping...
2021-03-01 00:00:00
no CS-2 data so skipping...
2021-04-01 00:00:00
no CS-2 data so skipping...

cs2_kk = []
valid_dates=[]
for date in date_range:
    
    try:
        xpts_cs2_kk, ypts_cs2_kk, cs2_kk_temp = getCS2kk(dataPathCS2, str(date.year)[-2:], '%02d' %(date.month))
    except:
        print(date)
        print('no CS-2 data so skipping...')
        continue
    valid_dates.append(date)
    
    cs2_kk_temp_is2grid = regridToICESat2(cs2_kk_temp, xpts_cs2_kk, ypts_cs2_kk, xptsIS2, yptsIS2) 
    cs2_kk_temp_is2grid_xr = xr.DataArray(data = cs2_kk_temp_is2grid, 
                            dims = ['y', 'x'], 
                            coords = {'latitude': (('y','x'), out_lats), 'longitude': (('y','x'), out_lons)}, 
                            name = 'CS-2_KK')

    cs2_kk.append(cs2_kk_temp_is2grid_xr)
                         
cs2_kk = xr.concat(cs2_kk, 'time')
cs2_kk = cs2_kk.assign_coords(time=valid_dates)
cs2_kk_attrs = {'units': 'meters', 'long_name': 'Kacimi and Kwok ICESat-2/CryoSat-2 monthly mean Arctic sea ice thickness', 'data_download': 'https://icesat-2.gsfc.nasa.gov/sea-ice-data/kacimi-kwok-2022', 
        'download_date': '09-2022', 'citation': 'Kacimi, S., Kwok, R. (2022), Arctic snow depth, ice thickness and volume from ICESat-2 and CryoSat-2: 2018-2021, Geophysical Research Letters, doi: 10.1029/2021GL097448.'}

cs2_kk = cs2_kk.assign_attrs(cs2_kk_attrs)  
# Add to book
book_ds['cs2_sea_ice_thickness_KK']=cs2_kk

/Users/aapetty/Data/CS2/KacimiKwok/thk_1811.txt
/Users/aapetty/Data/CS2/KacimiKwok/thk_1812.txt
/Users/aapetty/Data/CS2/KacimiKwok/thk_1901.txt
/Users/aapetty/Data/CS2/KacimiKwok/thk_1902.txt
/Users/aapetty/Data/CS2/KacimiKwok/thk_1903.txt
/Users/aapetty/Data/CS2/KacimiKwok/thk_1904.txt
/Users/aapetty/Data/CS2/KacimiKwok/thk_1905.txt
2019-05-01 00:00:00
no CS-2 data so skipping...
/Users/aapetty/Data/CS2/KacimiKwok/thk_1906.txt
2019-06-01 00:00:00
no CS-2 data so skipping...
/Users/aapetty/Data/CS2/KacimiKwok/thk_1907.txt
2019-07-01 00:00:00
no CS-2 data so skipping...
/Users/aapetty/Data/CS2/KacimiKwok/thk_1908.txt
2019-08-01 00:00:00
no CS-2 data so skipping...
/Users/aapetty/Data/CS2/KacimiKwok/thk_1909.txt
2019-09-01 00:00:00
no CS-2 data so skipping...
/Users/aapetty/Data/CS2/KacimiKwok/thk_1910.txt
/Users/aapetty/Data/CS2/KacimiKwok/thk_1911.txt
/Users/aapetty/Data/CS2/KacimiKwok/thk_1912.txt
/Users/aapetty/Data/CS2/KacimiKwok/thk_2001.txt
/Users/aapetty/Data/CS2/KacimiKwok/thk_2002.txt
/Users/aapetty/Data/CS2/KacimiKwok/thk_2003.txt
/Users/aapetty/Data/CS2/KacimiKwok/thk_2004.txt
/Users/aapetty/Data/CS2/KacimiKwok/thk_2005.txt
2020-05-01 00:00:00
no CS-2 data so skipping...
/Users/aapetty/Data/CS2/KacimiKwok/thk_2006.txt
2020-06-01 00:00:00
no CS-2 data so skipping...
/Users/aapetty/Data/CS2/KacimiKwok/thk_2007.txt
2020-07-01 00:00:00
no CS-2 data so skipping...
/Users/aapetty/Data/CS2/KacimiKwok/thk_2008.txt
2020-08-01 00:00:00
no CS-2 data so skipping...
/Users/aapetty/Data/CS2/KacimiKwok/thk_2009.txt
2020-09-01 00:00:00
no CS-2 data so skipping...
/Users/aapetty/Data/CS2/KacimiKwok/thk_2010.txt
/Users/aapetty/Data/CS2/KacimiKwok/thk_2011.txt
/Users/aapetty/Data/CS2/KacimiKwok/thk_2012.txt
/Users/aapetty/Data/CS2/KacimiKwok/thk_2101.txt
/Users/aapetty/Data/CS2/KacimiKwok/thk_2102.txt
/Users/aapetty/Data/CS2/KacimiKwok/thk_2103.txt
/Users/aapetty/Data/CS2/KacimiKwok/thk_2104.txt

Save data to local machine as a netcdf4 file#

# Could add an option to decide if you want to overwrite what's already there or not...

filename = './data/icesat2-cryosat-2-book-data.nc'
save_file = True

if (save_file == True):
    try: 
        book_ds.to_netcdf(path=filename, format='NETCDF4', mode='w')
        print('File ' + '"%s"' % filename + ' saved to directory ' + '"%s"' % os.getcwd())
    except: 
        print("Cannot save file because file by same name already exists")

else: 
    pass

Cannot save file because file by same name already exists