- Antarctic ice shelf calving is from Davison 2023
- Antarctic non-shelf calving is the difference between Davison 2023 and Rignot 2019
- Should be updated to Davison ESSD 2023 submitted
xr.open_dataset('./dat/AQ_calving.nc')
<xarray.Dataset> Size: 8kB
Dimensions: (region: 18, time: 25)
Coordinates:
* time (time) datetime64[ns] 200B 1997-07-01 1998-07-01 ... 2021-07-01
* region (region) int32 72B 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
Data variables:
calving (region, time) float64 4kB ...
uncertainty (region, time) float64 4kB ...
region_name (region) <U5 360B ...
Attributes:
description: Antarctic region ice shelf calving rate
date_created: 20241101T153743Z
title: Calving per region
history: Processed for Schmidt (YYYY; in prep); by Ken Mankoff
source: doi:10.5281/ZENODO.8052519
Conventions: CF-1.8
DOI: https://doi.org/10.5281/zenodo.14020895
import xarray as xr
ds = xr.open_dataset('./dat/AQ_calving.nc')
df = ds.sum(dim='region')['calving'].to_dataframe()
ax = df['calving'].plot(drawstyle='steps-post')
unc = (((ds['uncertainty']**2).sum(dim='region'))**0.5).to_dataframe()
unc.plot(drawstyle='steps-post', ax=ax)
_ = ax.set_ylabel('Calving [Gt yr$^{-1}$]')
![./fig/AQ_calving.png](/NASA-GISS/freshwater-forcing-workshop/raw/main/fig/AQ_calving.png)
import xarray as xr
import matplotlib.pyplot as plt
ds = xr.open_dataset('dat/AQ_calving.nc')
ds = ds.assign_coords(region=[_[0] + ' ['+str(_[1])+']' for _ in zip(ds['region_name'].values,ds['region'].values)])
df = ds['calving'].to_dataframe()
df = df['calving'].unstack().T
df['AQ'] = df.sum(axis='columns')
df.loc['Mean'] = df.mean(axis='rows')
df.index = [_[0:10] for _ in df.index.astype(str)]
df.round()
| A-Ap [1] | Ap-B [2] | B-C [3] | C-Cp [4] | Cp-D [5] | D-Dp [6] | Dp-E [7] | E-Ep [8] | Ep-F [9] | F-G [10] | G-H [11] | H-Hp [12] | Hp-I [13] | I-Ipp [14] | Ipp-J [15] | J-Jpp [16] | Jpp-K [17] | K-A [18] | AQ |
1997-07-01 | 56 | 37 | 42 | 93 | 141 | 114 | 26 | 43 | 108 | 83 | 200 | 37 | 48 | 45 | 10 | 139 | 93 | 48 | 1363 |
1998-07-01 | 56 | 37 | 42 | 93 | 141 | 114 | 26 | 43 | 108 | 83 | 200 | 37 | 48 | 45 | 10 | 1995 | 93 | 48 | 3219 |
1999-07-01 | 56 | 37 | 42 | 93 | 141 | 114 | 26 | 43 | 108 | 83 | 200 | 37 | 48 | 45 | 10 | 139 | 93 | 48 | 1363 |
2000-07-01 | 1 | 13 | 0 | 11 | 53 | 304 | 28 | 0 | 2944 | 19 | 129 | 22 | 215 | 245 | 68 | 2064 | 0 | 11 | 6125 |
2001-07-01 | 105 | 227 | 11 | 179 | 296 | 170 | 114 | 129 | 134 | 302 | 475 | 67 | 138 | 75 | 37 | 153 | 102 | 61 | 2776 |
2002-07-01 | 30 | 17 | 0 | 90 | 62 | 34 | 10 | 13 | 1 | 34 | 230 | 21 | 43 | 436 | 1 | 0 | 8 | 2 | 1033 |
2003-07-01 | 174 | 20 | 2 | 27 | 91 | 42 | 10 | 1058 | 0 | 51 | 176 | 23 | 33 | 27 | 1 | 0 | 10 | 3 | 1750 |
2004-07-01 | 104 | 55 | 64 | 155 | 178 | 135 | 42 | 68 | 108 | 130 | 289 | 48 | 83 | 64 | 11 | 15 | 48 | 40 | 1636 |
2005-07-01 | 15 | 21 | 2 | 24 | 109 | 28 | 18 | 0 | 25 | 52 | 60 | 34 | 46 | 347 | 2 | 0 | 37 | 26 | 850 |
2006-07-01 | 69 | 28 | 0 | 49 | 155 | 72 | 36 | 1 | 9 | 94 | 102 | 27 | 38 | 162 | 5 | 41 | 7 | 17 | 910 |
2007-07-01 | 19 | 27 | 11 | 30 | 173 | 80 | 34 | 9 | 38 | 85 | 103 | 40 | 42 | 12 | 3 | 0 | 10 | 24 | 739 |
2008-07-01 | 21 | 38 | 0 | 33 | 152 | 37 | 12 | 0 | 1 | 17 | 58 | 47 | 286 | 33 | 4 | 0 | 0 | 4 | 744 |
2009-07-01 | 61 | 51 | 31 | 94 | 204 | 136 | 48 | 17 | 76 | 147 | 1670 | 71 | 93 | 58 | 8 | 85 | 69 | 55 | 2974 |
2010-07-01 | 12 | 46 | 6 | 111 | 128 | 614 | 21 | 6 | 32 | 79 | 230 | 81 | 255 | 36 | 2 | 0 | 17 | 45 | 1722 |
2011-07-01 | 45 | 22 | 0 | 48 | 191 | 64 | 20 | 0 | 9 | 76 | 92 | 64 | 71 | 11 | 3 | 0 | 0 | 25 | 741 |
2012-07-01 | 3 | 20 | 12 | 48 | 121 | 74 | 21 | 0 | 20 | 97 | 180 | 39 | 54 | 4 | 4 | 0 | 14 | 10 | 722 |
2013-07-01 | 6 | 20 | 0 | 57 | 121 | 38 | 10 | 0 | 2 | 34 | 698 | 61 | 67 | 3 | 3 | 0 | 4 | 5 | 1130 |
2014-07-01 | 26 | 37 | 4 | 163 | 163 | 68 | 21 | 5 | 20 | 67 | 488 | 93 | 80 | 52 | 13 | 14 | 33 | 26 | 1374 |
2015-07-01 | 78 | 26 | 0 | 88 | 56 | 43 | 67 | 0 | 3 | 148 | 220 | 71 | 107 | 13 | 3 | 5 | 0 | 57 | 986 |
2016-07-01 | 37 | 20 | 0 | 42 | 116 | 29 | 19 | 7 | 10 | 51 | 302 | 34 | 49 | 14 | 2 | 26 | 9 | 25 | 791 |
2017-07-01 | 66 | 82 | 0 | 54 | 176 | 77 | 45 | 0 | 9 | 152 | 307 | 34 | 49 | 14 | 3 | 9 | 4 | 22 | 1102 |
2018-07-01 | 31 | 48 | 16 | 65 | 216 | 52 | 21 | 0 | 10 | 107 | 207 | 35 | 50 | 1325 | 3 | 11 | 2 | 20 | 2219 |
2019-07-01 | 25 | 36 | 0 | 60 | 126 | 43 | 20 | 6 | 3 | 89 | 361 | 36 | 43 | 32 | 4 | 6 | 7 | 21 | 917 |
2020-07-01 | 19 | 36 | 367 | 56 | 153 | 111 | 52 | 0 | 8 | 116 | 210 | 28 | 43 | 50 | 4 | 0 | 0 | 15 | 1269 |
2021-07-01 | 49 | 39 | 2 | 143 | 253 | 500 | 27 | 4 | 27 | 127 | 292 | 31 | 66 | 23 | 109 | 1019 | 7 | 176 | 2895 |
Mean | 47 | 42 | 26 | 76 | 149 | 124 | 31 | 58 | 153 | 93 | 299 | 45 | 84 | 127 | 13 | 229 | 27 | 33 | 1654 |
import numpy as np
import pandas as pd
import geopandas as gpd
import xarray as xr
import datetime
# shelf name with longitude and latitude
df = pd.read_excel("~/data/Davison_2023/adi0186_table_s2.xlsx",
sheet_name = 'Total mass changes',
usecols = (1,2,3), index_col = 0, skiprows = 4)
df = df.dropna()
shelf = gpd.GeoDataFrame(
geometry=gpd.points_from_xy(df.longitude, df.latitude, crs="EPSG:4326"), data=df)
shelf = shelf.to_crs('EPSG:3031')
# region name
region = gpd.read_file("~/data//IMBIE/Rignot/ANT_Basins_IMBIE2_v1.6.shp")
region = region[region['Regions'] != 'Islands']
# find regions nearest each shelf
shelf_region = gpd.sjoin_nearest(shelf,region)
shelf_region = shelf_region.drop(columns=['index_right','latitude','longitude','Regions'])
# load calving time series per shelf
calving = pd.read_excel("~/data/Davison_2023/adi0186_table_s2.xlsx",
sheet_name='Calving', index_col=1, skiprows=3, header=(0,1))
calving = calving.T.dropna().drop(columns=['Antarctic Ice Shelves'])
# WARNING: Calving < 0 implies error is baseline rate. This happens fairly often (small values) and occasionally (large values)
calving[calving < 0] = 0
obs = calving.xs('observed', level='Ice shelf')
obs.index.name = 'Date'
obs.index = pd.to_datetime(obs.index.astype(int).astype(str)+'-07-01', format="%Y-%m-%d")
# unc = calving.drop('observed', level=1, axis=0).reset_index().set_index('level_0').drop(columns=['Ice shelf'])
unc = calving.xs('uncertainty', level='Ice shelf')
unc.index.name = 'Date'
unc.index = obs.index
da_obs = xr.DataArray(data = obs.values,
dims = ['date','shelf'],
coords = {'date':obs.index.values, 'shelf':obs.columns})
ds = xr.Dataset({'calving': da_obs})
ds['uncertainty'] = (('date','shelf'), unc)
ds = ds.where(ds['shelf'] != 'Antarctic Ice Shelves', drop=True)
ds['region'] = (('shelf'), shelf_region['Subregion'])
# ds = ds.groupby('region').sum() # Want to agg() with different functions per column...
# uncertainty is sqrt of sum of squares. Not sure how to do this in-place in Xarray.
ds['unc2'] = ds['uncertainty']**2
ds2 = xr.merge([
ds[['calving','region']].groupby('region').sum(),
ds[['unc2','region']].groupby('region').sum(),
])
ds2['uncertainty'] = np.sqrt(ds2['unc2'])
ds2 = ds2.drop_vars('unc2')
# uncertainty for all of AQ as (sum(u**2))**0.5 matches Davison 2023 row 168 "Antarctic Ice Shelves"
# need to calculate AQ-wide uncertainty at shelf resolution because step-aggregating is not commutative
# ds2['uncertainty_AQ'] = np.sqrt(ds['unc2'].sum(dim='shelf'))
ds = ds2
ds = ds.rename({'date':'time'})
ds['region'] = np.arange(18).astype(np.int32) + 1
ds['region_name'] = (('region'), ['A-Ap', 'Ap-B', 'B-C', 'C-Cp', 'Cp-D',
'D-Dp', 'Dp-E', 'E-Ep', 'Ep-F', 'F-G',
'G-H', 'H-Hp', 'Hp-I', 'I-Ipp', 'Ipp-J',
'J-Jpp', 'Jpp-K', 'K-A'])
ds.attrs['description'] = 'Antarctic region ice shelf calving rate'
ds['calving'].attrs['units'] = 'Gt yr-1'
ds['calving'].attrs['long_name'] = 'Shelf calving'
# ds['calving'].attrs['standard_name'] = 'water_flux_into_sea_water_from_land_ice'
# https://github.com/orgs/cf-convention/discussions/388
ds['calving'].attrs['standard_name'] = 'ice_transport_across_line'
ds['uncertainty'].attrs['long_name'] = 'Uncertainty of shelf calving'
ds['time'].attrs['standard_name'] = 'time'
ds['region'].attrs['long_name'] = 'IMBIE region'
ds.attrs['date_created'] = datetime.datetime.now(datetime.timezone.utc).strftime("%Y%m%dT%H%M%SZ")
ds.attrs['title'] = 'Calving per region'
ds.attrs['history'] = 'Processed for Schmidt (YYYY; in prep); by Ken Mankoff'
ds.attrs['source'] = 'doi:10.5281/ZENODO.8052519'
ds.attrs['Conventions'] = 'CF-1.8'
ds.attrs['DOI'] = 'https://doi.org/10.5281/zenodo.14020895'
comp = dict(zlib=True, complevel=5)
encoding = {}
encoding['time'] = {'dtype': 'i4'}
!rm ./dat/AQ_calving.nc
ds.to_netcdf('./dat/AQ_calving.nc', encoding=encoding)
!ncdump -h ./dat/AQ_calving.nc
netcdf AQ_calving {
dimensions:
region = 18 ;
time = 25 ;
variables:
double calving(region, time) ;
calving:_FillValue = NaN ;
calving:units = "Gt yr-1" ;
calving:long_name = "Shelf calving" ;
calving:standard_name = "ice_transport_across_line" ;
int time(time) ;
time:standard_name = "time" ;
time:units = "days since 1997-07-01 00:00:00" ;
time:calendar = "proleptic_gregorian" ;
int region(region) ;
region:long_name = "IMBIE region" ;
double uncertainty(region, time) ;
uncertainty:_FillValue = NaN ;
uncertainty:long_name = "Uncertainty of shelf calving" ;
string region_name(region) ;
// global attributes:
:description = "Antarctic region ice shelf calving rate" ;
:date_created = "20241101T153743Z" ;
:title = "Calving per region" ;
:history = "Processed for Schmidt (YYYY; in prep); by Ken Mankoff" ;
:source = "doi:10.5281/ZENODO.8052519" ;
:Conventions = "CF-1.8" ;
:DOI = "https://doi.org/10.5281/zenodo.14020895" ;
}