Gravimetric Mass Balance for Ice Sheet

logo

Gravimetric Mass Balance for Ice Sheet#

By Sebastian B. Simonsen and Natalia Havelund, DTU Space, Comments can be posted to S. Simonsen at ssim@dtu.dk

About#

This tutorial will demonstrate how to plot the Gravimetric Mass Balance (GMB) data for selected regions across the Greenland or Antarctic ice sheets, using data from the Copernicus Climate Change Service (C3S).

It will show you how to download data from the C3S Climate Data Store (CDS), calculate average change rates over a selectable period, and display the results.

Prepare your environment#

To download data from the CDS, you must complete three steps;

  1. Register for free with the CDS at https://cds.climate.copernicus.eu/#!/home

  2. Set up the CDS Application Program Interface (API), as described at https://cds.climate.copernicus.eu/api-how-to

  3. Agree to the dataset terms of use - if you have not already done so you will be directed to the correct webpage the first time you try to download data. Accept the terms and conditions, and re-run the downloading code.

Import libraries#

The data is stored in netCDF4 files. To work with them, several libraries are needed. These are used to retrieve the data, unpack it, run calculations and plot the results.

%load_ext autoreload
%autoreload 2

# CDS API
import cdsapi

# File handling
import shutil
import glob
import os

# Calculations
import numpy as np
import xarray as xr

# Mapping
import cartopy.crs as ccrs
import cartopy.feature as cfeature

# Plotting
import matplotlib.pyplot as plt

# Ignore distracting warnings
import warnings
warnings.filterwarnings('ignore')

Download data#

Select the ice sheet of interest and retrieve the dataset. An example of how to retrieve the data is given on the CDS, under ‘Show API request’ on the dataset’s webpage at https://cds.climate.copernicus.eu/cdsapp#!/dataset/satellite-ice-sheet-elevation-change?tab=form

The dataset comes in a zip file, and may include several files. Only the most recent is necessary - this contains all of the data in previous files, plus the latest update.

# Set ice_sheet to the desired ice sheet, either 'AntIS' or 'GrIS'

# This will be used to select data, and later to set up plot defaults

#ice_sheet = 'AntIS' # 

ice_sheet = 'GrIS' # 

basin_num = 0 # if set 0, then mass loss for all is shown. For GrIS 0-8, for AntIS 0-27 basins

# Download data from the CDS

import cdsapi
# Download data
# We will use the [Climate Data Store (CDS) API](https://cds.climate.copernicus.eu/api-how-to) to download the [glacier mass change](https://cds.climate.copernicus.eu/cdsapp#!/dataset/derived-gridded-glacier-mass-change) dataset. 
# _NOTE: To use the CDS API, you first need to [register](https://cds.climate.copernicus.eu/user/register) (if not already), find your `UID` and `API key` on your [acount page](https://cds.climate.copernicus.eu/user), and fill them in below._
# Fill in your UID and API key (separated by a colon :)

KEY = ':'  
c = cdsapi.Client(key=KEY)

c.retrieve(
    'satellite-ice-sheet-mass-balance',
    {
        'variable': 'all',
        'format': 'zip',
    },
    'download.zip')

# Unpack the zip file and remove the zipped data

shutil.unpack_archive('download.zip', '.')
os.remove('download.zip')

# # List the files contained and select the latest netCDF (.nc) file. Print its filename

list_of_files = glob.glob('*.nc') # * means all if need specific format then *.nc
latest_file = max(list_of_files, key=os.path.getctime)
print('In the following you are seeing the data from: \n'+latest_file)

Examine the data#

Open the data file and list its contents. These include descriptive comments.

# Open the NetCDF file

ds = xr.open_dataset(latest_file)

# Print contents list
print('The data file contains the following:')
print(ds)

# Load data
time_in_hours = ds['time']

# Convert xarray DataArray to NumPy array
time_in_hours_array = time_in_hours.values

if 'GrIS' in ice_sheet:
    drainage_basin = ds['GrISBasin']


if 'AntIS' in ice_sheet:
    drainage_basin = ds['AntBasin']    

if basin_num != 0:        
    filtered_basins = drainage_basin.where((drainage_basin[:, 2] - basin_num).astype(int) == 0, drop=True)
    basin_numbers_str = filtered_basins[:, 2].astype(str)
    filtered_basins_2x = filtered_basins.where(basin_numbers_str.str.startswith(str(basin_num)+'.'), drop=True)
    
    # Extracting latitude, longitude, and basin number
    
    lat, lon, basin_number = filtered_basins_2x[:,0].values,filtered_basins_2x[:,1].values,filtered_basins_2x[:,2].values
 
    gmb = ds[str(ice_sheet)+'_'+str(round(basin_num))]
    gmb_er = ds[str(ice_sheet)+'_'+str(round(basin_num))+'_er']


else:
    
    lat, lon = drainage_basin[:,0].values,drainage_basin[:,1].values
    gmb = ds[str(ice_sheet)+'_total']
    gmb_er = ds[str(ice_sheet)+'_total_er']

# Create a figure and axis with a specific projection (e.g., Orthographic projection centered on Greenland)


fig, (ax1, ax2) = plt.subplots(1, 2)
if basin_num != 0:
    fig.suptitle('Mass change for basin number: '+str(basin_num))
else:
    fig.suptitle('Total mass change')
# Remove x-axis labels from the first subplot
ax1.set_xticks([])
ax1.set_yticks([])

# Remove outline around the subplot
ax1.spines['top'].set_visible(False)
ax1.spines['bottom'].set_visible(False)
ax1.spines['right'].set_visible(False)
ax1.spines['left'].set_visible(False)


# Subplot 1 - Map

if 'GrIS' in ice_sheet:
    # lon = ds['lon']
    # lat = ds['lat']
    # dhdt_ave = np.mean(ds['dhdt'][:,:,Tbool], axis=2)
    ext= [-60, -25, 57, 84]

    # ax1 = plt.axes(projection=ccrs.Orthographic(-45, 71))
    ax1 = plt.subplot(1, 2, 1, projection=ccrs.Orthographic(-45, 71))


if 'AntIS' in ice_sheet:
    # lon = ds['longitude']
    # lat = ds['latitude']
    # dhdt_ave = np.mean(ds['sec'][:,:,Tbool], axis=2)
    ext = [-180, 180, -90, -60]
    
    # ax1 = plt.axes(projection=ccrs.Orthographic(0, -90)) # ?
    ax1 = plt.subplot(1, 2, 1, projection=ccrs.Orthographic(0, -90))


# ax1 = plt.subplot(1, 2, 1, projection=ccrs.Orthographic(-45, 71))
ax1.add_feature(cfeature.COASTLINE.with_scale('10m'), linewidth=0.5)
ax1.set_extent(ext, crs=ccrs.PlateCarree())
ax1.add_feature(cfeature.LAND, facecolor='white')
ax1.add_feature(cfeature.OCEAN, facecolor='gray')
# Add latitude and longitude gridlines
gl = ax1.gridlines(draw_labels=True, linewidth=0.5, color='gray', alpha=0.5, linestyle='--')

# Remove latitude labels at the top
gl.top_labels = False

# Set the label spacing to display labels only once on either side
gl.xlabels_top = False
gl.xlabels_bottom = True
gl.ylabels_left = True
gl.ylabels_right = False
gl.xlabel_style = {'size': 10}
gl.ylabel_style = {'size': 10}

# Sort coordinates by longitude
sorted_indices = np.argsort(lon)
lon_sorted = lon[sorted_indices]
lat_sorted = lat[sorted_indices]
 # Plot the sorted coordinates
ax1.plot(lon_sorted, lat_sorted, 'r', linestyle='dashed', transform=ccrs.PlateCarree())
ax1.fill(lon_sorted, lat_sorted, color='red', transform=ccrs.PlateCarree())#, alpha=0.5)

# Calculate the midpoint
midpoint_lat = np.mean(lat)
midpoint_lon = np.mean(lon)
# Add text to the middle of the lat-lon set
if basin_num != 0:
    ax1.text(midpoint_lon-1, midpoint_lat, str(basin_num), color='white', transform=ccrs.PlateCarree())
else:
    ax1.text(midpoint_lon-1, midpoint_lat, 'All', color='white', transform=ccrs.PlateCarree())
    ax1.add_feature(cfeature.LAND, facecolor='red')




# Subplot 2 - Time series


ax2 = plt.subplot(1, 2, 2)
ax2.plot(time_in_hours_array, gmb, '-', color='grey', label='Mass change')
ax2.fill_between(time_in_hours_array, gmb - gmb_er, gmb + gmb_er, alpha=0.3, color='grey', label='Error')
ax2.set_xlabel('Time')
ax2.set_ylabel('Mass change [GT]')


# Set the figure size to a specific value
desired_fig_size = (12, 5)
fig.set_size_inches(desired_fig_size)


# Display the plot

plt.show()
Contents