Table of contents

• Introduction
• Mass loss from citet:mankoff_2021
  • Baseline mass loss (calendar years)
  • 1991 – 2020 hydro years
  • 2023 mass loss
• Discharge
  • Baseline discharge (1991 through 2020)
  • Average since 2013
  • 2023 discharge through latest update
  • Trends
    • All GIS
    • By region
    • Publication graphic
• Greenland outline
• Bare ice area
• Albedo
  • Crop to GL
  • Figure
• Melt
  • Figure
• PROMICE In situ / Point obs

Introduction

GISS contributions to the Arctic Report Card

Mass loss from citet:mankoff_2021

# v708 fetched on [2023-11-03 Fri] https://dataverse.geus.dk/api/access/datafile/:persistentId?persistentId=doi:10.22008/FK2/OHI23Z/4KAVFS -O MB_region.nc

Baseline mass loss (calendar years)

• ~ 218 Gt yr-1

%cd '/home/kdm/projects/arctic_report_card'

import xarray as xr
ds = xr.open_dataset("MB_region.nc")

print(ds\
      .sum(dim='region')\
      .sel({'time':slice('2002-01-01','2022-12-31')})\
      .mean()*365)

/home/kdm/projects/arctic_report_card
<xarray.Dataset>
Dimensions:      ()
Data variables: (12/19)
    MB           float64 -218.2
    MB_err       float64 91.57
    MB_ROI       float64 -218.2
    MB_ROI_err   float64 156.8
    SMB          float64 291.5
    SMB_err      float64 26.23
    ...           ...
    BMB_err      float64 5.586
    BMB_ROI      float64 23.95
    BMB_ROI_err  float64 5.497
    MB_HIRHAM    float64 -222.3
    MB_MAR       float64 -210.1
    MB_RACMO     float64 -218.6

1991 – 2020 hydro years

ds.sum(dim='region')\
  .sel({'time':slice('1990-09-01','2020-08-31')})\
  .resample({'time':'AS-SEP'})\
  .sum()\
  .mean()\
  [['MB','MB_err','SMB','SMB_err','BMB','BMB_err','D','D_err']]

<xarray.Dataset>
Dimensions:  ()
Data variables:
    MB       float32 -162.4
    MB_err   float32 87.96
    SMB      float32 325.3
    SMB_err  float32 29.28
    BMB      float32 23.06
    BMB_err  float32 5.501
    D        float32 464.6
    D_err    float32 42.68

2023 mass loss

ds.sum(dim='region')\
  .sel({'time':slice('2022-09-01','2023-08-31')})\
  .sum()\
  [[_ for _ in ds.data_vars if 'err' not in _]]

<xarray.Dataset>
Dimensions:    ()
Data variables:
    MB         float32 -214.1
    MB_ROI     float32 -214.1
    SMB        float32 317.1
    SMB_ROI    float32 317.1
    D          float32 504.1
    D_ROI      float32 504.1
    BMB        float32 27.2
    BMB_ROI    float32 27.2
    MB_HIRHAM  float32 -136.7
    MB_MAR     float32 -291.6
    MB_RACMO   float32 0.0

As per ARC 2023 GRACE for 2023 hydro year [2022-09-01 Thu] - [2023-08-31 Thu] is -156 +- 22 Gt.

• NOTE: This number is already scaled by 0.84 so it represents only main ice sheet.

hydro year SMB from ‘observations’

• [ ] -156 + 504 + 27 = 375
• [ ] -291 + 504 + 27 = 240 # MAR SMB = MB_MAR + D + BMB
• [ ] -136 + 504 + 27 = 395 # HIRHAM SMB = MB_HIRHAM + D + BMB

Discharge

• Using v90 from https://doi.org/10.22008/promice/data/ice_discharge/d/v02

Baseline discharge (1991 through 2020)

import xarray as xr
ds = xr.open_dataset("~/data/Mankoff_2020/ice/v90/region.nc").sum(dim='region')

df = ds[['discharge','err']].to_dataframe()['1990':]
df['1991-01-01':'2020-12-31':].mean()
df['1991-01-01':'2020-12-31':].resample('1D').interpolate().sum()/30/365
df['1991-01-01':'2020-12-31':].resample('1D').interpolate().mean()

discharge    465.048117
err           42.717198
dtype: float64

Average since 2013

df['2013-01-01':].resample('1D').interpolate().mean()

discharge    500.629915
err           46.894040
dtype: float64

2023 discharge through latest update

print("Last timestamp: ", df.index[-1])
df['2023-01-01':'2023-12-31'].resample('1D').interpolate().mean()

Last timestamp:  2023-08-10 00:00:00
discharge    503.397368
err           47.069253
dtype: float64

Trends

See ./figs_tmp sub-folder for graphics

All GIS

df['discharge'].resample('1D').interpolate().resample('YS').mean().plot(drawstyle='steps-post')
df['discharge'].resample('1D').interpolate().resample('YS').mean().tail()

time
2019-01-01    503.009784
2020-01-01    513.816107
2021-01-01    517.558611
2022-01-01    514.149532
2023-01-01    503.431076
Freq: AS-JAN, Name: discharge, dtype: float64

By region

dsR = xr.open_dataset("~/data/Mankoff_2020/ice/v90/region.nc")

# dsR = dsR['discharge'].resample({'time':'1D'}).interpolate().resample({'time':'MS'}).mean()
dsR = dsR['discharge'].resample({'time':'1D'}).interpolate().resample({'time':'YS'}).mean()
_ = dsR.plot.line(x='time', drawstyle='steps-post')

Publication graphic

import matplotlib.gridspec as gridspec
import numpy as np
import pandas as pd
from adjust_spines import adjust_spines as adj
import matplotlib.pyplot as plt
import datetime as dt

from cycler import cycler
plt.rcParams['axes.prop_cycle'] = cycler('color', ['#1f77b4', '#ff7f0e', '#2ca02c', '#d62728', \
                                                   '#9467bd', '#8c564b', '#e377c2', '#bcbd22', '#17becf'])

fig = plt.figure(1, figsize=(9,7)) # w,h
fig.clf()
grid = plt.GridSpec(2, 1, height_ratios=[1,6], hspace=0.1) # h, w

ax_D = fig.add_subplot(grid[1,:])

from adjust_spines import adjust_spines as adj
adj(ax_D, ['left','bottom'])

ROOT="./out/"
ROOT="/home/kdm/data/Mankoff_2020/ice/v90/"
D = pd.read_csv(ROOT+"region_D.csv", index_col=0, parse_dates=True)
err = pd.read_csv(ROOT+"region_err.csv", index_col=0, parse_dates=True)
coverage = pd.read_csv(ROOT+"region_coverage.csv", index_col=0, parse_dates=True)

THRESH = coverage < 0.5
D[THRESH] = np.nan
err[THRESH] = np.nan
coverage[THRESH] = np.nan

# PROMICE drop in SE. Need 200 m data
D = D.iloc[:-5]
err = err.iloc[:-5]
coverage = coverage.iloc[:-5]

def pad_df(df):
    df = pd.concat([pd.DataFrame(index=np.array(['1986-01-01']).astype('datetime64[ns]')), df] )
    idx = str(df.index.year.max())+'-12-31'
    df = pd.concat([df, pd.DataFrame(index=np.array([idx]).astype('datetime64[ns]'))])
    df = df.sort_index()
    return df

D = pad_df(D)
err = pad_df(err)
coverage = pad_df(coverage)

### Take annual average from daily interpolated rather than the existing samples.
D_day_year = D.resample('1D',axis='rows').mean().interpolate(method='time',limit_area='inside').resample('A',axis='rows').mean()
err_day_year=err.resample('1D',axis='rows').mean().interpolate(method='time',limit_area='inside').resample('A',axis='rows').mean()

# No annual average if few sample
num_obs = D.resample('Y').count().values
D_day_year[num_obs<=3] = np.nan
err_day_year[num_obs<=3] = np.nan

MS=4
Z=99
for r in D.columns:
    e = ax_D.errorbar(D[r].index, D[r].values, fmt='o', mfc='none', ms=MS)
    C = e.lines[0].get_color()
    D_day_year[r].plot(drawstyle='steps', linewidth=2, ax=ax_D,
                       color=C,
                       alpha=0.75, zorder=Z)
    for i in np.arange(D.index.size):
        if np.isnan(D.iloc[i][r]): continue
        alpha = coverage.iloc[i][r]
        if alpha < 0: alpha = 0
        if alpha > 1: alpha = 1
        ax_D.errorbar(D.iloc[i].name, D.iloc[i][r],
                      yerr=err.iloc[i][r], ecolor='gray',
                      marker='o', ms=MS,
                      # mfc='k', mec='k',
                      color=C,
                      mfc=C, mec=C,
                      alpha=alpha)

    tx = pd.Timestamp(str(D[r].dropna().index[-1].year) + '-01-01') + dt.timedelta(days=380)
    ty = D_day_year[r].dropna().iloc[-1]
    # if r in ['CE', 'SW']: ty=ty-4
    if r == 'CE': ty=ty-4
    # if r == 'NE': ty=ty+4
    # if r == 'NO': ty=ty-2
    ax_D.text(tx, ty, r, verticalalignment='center', horizontalalignment='left')

import matplotlib.dates as mdates
ax_D.xaxis.set_major_locator(mdates.YearLocator())

# plt.legend()
ax_D.legend("", framealpha=0)
ax_D.set_xlabel('Time [Years]')
ax_D.set_ylabel('Discharge [Gt yr$^{-1}$]')
ax_D.set_xlim(D.index[0], D.index[-1])
ax_D.set_xticklabels(D.index.year.unique())

ax_D.xaxis.set_tick_params(rotation=-90)
for tick in ax_D.xaxis.get_majorticklabels():
    tick.set_horizontalalignment("left")

plt.savefig('./discharge_ts_regions.png', transparent=False, bbox_inches='tight', dpi=300)
plt.savefig('./discharge_ts_regions.svg', transparent=False, bbox_inches='tight', dpi=300)

Err_pct = (err_day_year.values/D_day_year.values*100).round().astype(int).astype(str)
Err_pct[Err_pct.astype(float)<0] = 'NaN'
tbl = (D_day_year.round().fillna(value=0).astype(int).astype(str) + ' ('+Err_pct+')')
tbl.index = tbl.index.year.astype(str)
tbl.columns = [_ + ' (Err %)' for _ in tbl.columns]
tbl

<Figure size 900x700 with 1 Axes>

Greenland outline

grass -c EPSG:3413 G_3413

v.import input=/home/kdm/data.me/GIS/NaturalEarth/ne_10m_admin_0_countries/ne_10m_admin_0_countries.shp output=countries
v.extract input=countries output=greenland where='name = "Greenland"'
v.out.ogr input=greenland output=greenland.gpkg

v.import input=/home/kdm/data/Zwally_2012/sectors/sectors.shp output=zwally_2012
g.region vector=zwally_2012 res=100 -ap
v.to.rast input=zwally_2012 output=z_rast use=val val=1
r.to.vect input=z_rast output=ice_edge type=area
v.out.ogr input=ice_edge output=ice_edge.gpkg

Bare ice area

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import xarray as xr
import datetime

from matplotlib import rc
rc('font', size=11)
rc('text', usetex=False)
# matplotlib.pyplot.xkcd()

# plt.close(1)
fig = plt.figure(1, figsize=(5,4)) # w,h
fig.clf()
fig.set_tight_layout(True)
import matplotlib.gridspec as gridspec

ax = fig.add_subplot(111)
colors = ['purple','k', 'r', 'darkorange', 'b', 'g','lightgreen']

ds = xr.open_mfdataset('./Adrien/SICE_GrIS_bare_ice_area_*.nc')
df = ds.to_dataframe()

this_y = datetime.datetime.now().year

for i,y in enumerate(df.index.year.unique()[::-1]):
    data = df[df.index.year == y]
    data = data.resample('1D').ffill()
    data = data[(data.index.dayofyear > 130) & (data.index.dayofyear < 267)]
    ax.plot(data.index.dayofyear,
            data['bare_ice_area_km2'],
            # drawstyle='steps-post',
            color=colors[i],
            linewidth = (2 if y == this_y else 1),
            label=str(y))

ax.legend(fontsize=9, frameon=True, bbox_to_anchor=(0, 0.9), loc='upper left')

from adjust_spines import adjust_spines as adj
adj(ax, ['left','bottom'])

ax.set_ylabel('Bare ice area [km$^{2}$]')
import matplotlib.dates as mdates

label = data.index[(data.index.day == 1) | (data.index.day == 15)]
ax.set_xticks(label.dayofyear)
ax.set_xticklabels([str(_)[5:10] for _ in label])
ax.set_xticklabels(['May 15','June 1','June 15','July 1','July 15','Aug 1','Aug 15','Sep 1','Sep 15'])
plt.xticks(rotation=45)


# ax.get_yaxis().set_major_formatter(
#     mpl.ticker.FuncFormatter(lambda x, p: format(int(x), ',')))

ax.grid(visible=True, which='major', axis='y', alpha=0.33)
ax.grid(visible=True, which='major', axis='x', alpha=0.33)

plt.savefig('bare_ice.png', transparent=False, bbox_inches='tight', dpi=300)
plt.savefig('bare_ice.svg', transparent=False, bbox_inches='tight', dpi=300)

Albedo

Crop to GL

grass -c ./G_3413/AW

g.region vector=greenland@PERMANENT res=500 -pa

r.import input=Adrien/SICE_2023_JJA_albedo_anomaly_vs_2017_2022.tif output=anom extent=input

# d.mon wx0
# d.rast anom

eval $(g.region -upg raster=anom)

r.mask vector=greenland@PERMANENT
g.region zoom=MASK
r.mapcalc "cropped = anom"


r.mask -r
g.region raster=cropped -pa # set bounds based on crop
g.region e=$e w=$w -pa # expand e/w to original to include Canada
r.mapcalc "cropped_NS = anom"

g.region raster=cropped
r.out.gdal input=cropped output=Adrien/cropped.tif format=GTiff createopt="COMPRESS=DEFLATE"

g.region raster=cropped_NS
r.out.gdal input=cropped_NS output=Adrien/cropped_NS.tif format=GTiff createopt="COMPRESS=DEFLATE"

Figure

import matplotlib
import matplotlib.pyplot as plt
import datetime
import numpy as np
import pandas as pd
import geopandas as gp
import rasterio as rio
import rasterio.mask
from rasterio.plot import plotting_extent
import cmocean
from mpl_toolkits.axes_grid1.inset_locator import inset_axes

from matplotlib import rc
rc('font', size=11)
rc('text', usetex=False)
# matplotlib.pyplot.xkcd()

C_land = "#EAEAEA"
C_ocean = "#D0CFD4"

# plt.close(1)
fig = plt.figure(1, figsize=(8,8)) # w,h
#gcm = get_current_fig_manager()
#gcm.window.move(-1000,0)
#gcm.resize(gcm.window.size().height(), gcm.window.size().width())
# get_current_fig_manager().window.move(0,0)
fig.clf()
# fig.set_tight_layout(True)
import matplotlib.gridspec as gridspec

gs = gridspec.GridSpec(2,2, width_ratios=[1,1], height_ratios=[5,1]) #w,h



ax_albedo_map = plt.subplot(gs[0,0])
ax_albedo_plot = plt.subplot(gs[1,0])

if 'o' not in locals():
    o = gp.read_file('greenland.gpkg')
    
o.plot(color=C_land, ax=ax_albedo_map, facecolor='none', zorder=-1)
# o.plot(facecolor='None', edgecolor='gray', ax=ax_albedo_map, zorder=-1, alpha=1, linewidth=0.5)


r_albedo = rio.open('./Adrien/cropped.tif')
r_albedo_extent = plotting_extent(r_albedo)
r_albedo = r_albedo.read(1)
r_albedo[r_albedo== -999] = np.nan

cmapGr = matplotlib.cm.get_cmap(plt.cm.BrBG_r)
cmapBl = matplotlib.cm.get_cmap(plt.cm.RdBu)
colors = np.vstack(([cmapGr(i) for i in np.arange(128,257)[::-1]], [cmapBl(i) for i in np.arange(128,257)]))
import matplotlib.colors as mcolors
cmap = mcolors.LinearSegmentedColormap.from_list('my_colormap', colors)
# cmap = matplotlib.cm.get_cmap(cmocean.cm.balance_r)

im_albedo = ax_albedo_map.imshow(r_albedo, extent=r_albedo_extent, cmap=cmap, vmin=-0.1, vmax=0.1)

ax_albedo_map.axis('off')

ax_albedo_cb = inset_axes(ax_albedo_map,
                          width="5%",  # width = 5% of parent_bbox width
                          height="25%",  # height : 50%
                          loc='lower right',
                          bbox_to_anchor=(-0.25, 0, 1, 1),
                          bbox_transform=ax_albedo_map.transAxes,
                          borderpad=0)

# cb_albedo = fig.colorbar(im, cax=ax_albedo_cb)
cb_albedo = fig.colorbar(im_albedo, cax=ax_albedo_cb)
cb_albedo.set_label('Albedo anomaly\n[unitless]')

# df = pd.read_csv('JEB/MODIS_S3_JJA.csv',
#                   parse_dates=True, index_col=1)
# df.index = [datetime.datetime(int(_),1,1) for _ in df.index]
# df.loc[df.index[-1] + (df.index[-1]-df.index[-2])] = df.iloc[-1]
# ax_albedo_plot.plot(df.index,
#                     # df.sum(axis='columns').values,
#                     df['JJA_MODIS_S3'].values,
#                     drawstyle='steps-post', color='k')


df = xr.open_dataset('./Adrien/MODIS_GrIS_JJA_mean_albedo.nc').to_dataframe()
df.index = [datetime.datetime(int(_),1,1) for _ in df.index]
df.loc[df.index[-1] + (df.index[-1]-df.index[-2])] = df.iloc[-1]
ax_albedo_plot.plot(df.index,
                    df['JJA_GrIS_mean_albedo_MODIS'].values,
                    drawstyle='steps-post', color='k')

# df = xr.open_dataset('./Adrien/SICE_GrIS_JJA_mean_albedo.nc').to_dataframe()
# df.index = [datetime.datetime(int(_),1,1) for _ in df.index]
# df.loc[df.index[-1] + (df.index[-1]-df.index[-2])] = df.iloc[-1]
# ax_albedo_plot.plot(df.index,
#                     df['JJA_mean_albedo'].values,
#                     drawstyle='steps-post', color='g')



adj(ax_albedo_plot, ['left','bottom'])

ax_albedo_plot.set_ylim(0.76,0.81)
ax_albedo_plot.set_yticks([0.76, 0.77, 0.78, 0.79, 0.80, 0.81])
# ax_albedo_plot.spines['left'].set_bounds(0.74, 0.82)
ax_albedo_plot.set_ylabel('Albedo\n[unitless]')
ax_albedo_plot.set_xticks(ax_albedo_plot.get_xticks()+365*2+1)
# # ax_albedo_plot.xticks(rotation=70)
# # plt.setp(ax_albedo_plot.xaxis.get_majorticklabels(), rotation=70)
import matplotlib.dates as mdates
ax_albedo_plot.xaxis.set_major_formatter(mdates.DateFormatter('%Y'))

ax_albedo_plot.grid(visible=True, which='major', axis='y', alpha=0.33)
ax_albedo_plot.plot(df.index[[0,-1]], [df['JJA_GrIS_mean_albedo_MODIS'].mean()]*2, 'k--', alpha=0.5)

plt.savefig('albedo.png', transparent=False, bbox_inches='tight', dpi=300)
plt.savefig('albedo.svg', transparent=False, bbox_inches='tight', dpi=300)

Melt

ls TM

grass -c ./G_3413/TM

g.region vector=greenland@PERMANENT res=500 -pa

r.import input=TM/greenland_melt_anomaly_colorless_20230401-20230831.tif output=melt extent=input

# d.mon wx0
# d.rast melt

eval $(g.region -upg raster=melt)

r.mask vector=greenland@PERMANENT
g.region zoom=MASK
r.mapcalc "cropped = melt"

g.region raster=cropped
r.out.gdal input=cropped output=TM/cropped.tif format=GTiff createopt="COMPRESS=DEFLATE"

Figure

import numpy as np
import pandas as pd
import geopandas as gp
import rasterio as rio
import rasterio.mask
import matplotlib
import matplotlib.pyplot as plt
from rasterio.plot import plotting_extent
import cmocean
from mpl_toolkits.axes_grid1.inset_locator import inset_axes

from matplotlib import rc
rc('font', size=12)
rc('text', usetex=False)

fig = plt.figure(1, figsize=(8,8)) # w,h
fig.clf()

gs = gridspec.GridSpec(2,2, width_ratios=[1,1], height_ratios=[4,1]) #w,h

ax_melt_map = plt.subplot(gs[0,1])
ax_melt_plot = plt.subplot(gs[1,1])

C_land = "#EAEAEA"
C_ocean = "#D0CFD4"

# ax_melt_map.set_facecolor(C_ocean)

if 'r_melt' not in locals():
    r_melt = rio.open('./TM/cropped.tif')
    r_melt_extent = plotting_extent(r_melt)

    r_melt = r_melt.read(1)
    r_melt[r_melt== -999] = np.nan

if 'o' not in locals():
    o = gp.read_file('greenland.gpkg')
    
o.plot(color=C_land, ax=ax_melt_map, facecolor='none', zorder=-1)

cmap = matplotlib.cm.get_cmap(cmocean.cm.balance)
im_melt = ax_melt_map.imshow(r_melt, extent=r_melt_extent,
                             cmap=cmap,
                             vmin=-40, vmax=40)

ax_melt_map.axis('off')

ax_melt_cb = inset_axes(ax_melt_map,
                        width="5%",  # width = 5% of parent_bbox width
                        height="25%",  # height : 50%
                        loc='lower right',
                        bbox_to_anchor=(-0.25, 0, 1, 1),
                        bbox_transform=ax_melt_map.transAxes,
                        borderpad=0)


cb_melt = fig.colorbar(im_melt, cax=ax_melt_cb)
cb_melt.set_label('Melt anomaly\n[days]')


df0 = pd.read_csv('TM/greenland-daily-melt.csv', parse_dates=True, index_col=0)
df1 = pd.read_csv('TM/greenland-daily-melt-climatology.csv')
df1['date'] = [pd.to_datetime('2023-01-01') + pd.to_timedelta(doy-1, unit='D') for doy in df1['doy']]
df1 = df1.set_index('date')
df = df0.merge(df1, left_index=True, right_index=True)
df[df['qc_flag'] != True] = np.nan

df = df.apply(lambda x: x/df['icesheet_area_km2_x']*100)

ax_melt_plot.plot(df['Median'], color='k', linestyle='--', drawstyle='steps-post', label='Median')
ax_melt_plot.plot(df['melting_area_km2'],
         color=np.array(cmap(185, bytes=True)[0:3])/255,
         drawstyle='steps-post',
         label='2023',
         linewidth=1.0)

ax_melt_plot.fill_between(df.index,
                 df['10'].values.flatten(),
                 df['90'].values.flatten(),
                 color='gray',
                 step='post',
                 label='Interdecile range',
                 alpha=0.25)

ax_melt_plot.fill_between(df.index,
                 df['25'].values.flatten(),
                 df['75'].values.flatten(),
                 color='k',
                 step='post',
                 label='Interquartile range',
                 alpha=0.25)

ax_melt_plot.legend(fontsize=9, frameon=False, bbox_to_anchor=(0, 1.25), loc='upper left', ncol=2)

from adjust_spines import adjust_spines as adj
adj(ax_melt_plot, ['left','bottom'])

ax_melt_plot.set_ylim(0,60)
ax_melt_plot.set_yticks([0,20,40,60])
ax_melt_plot.spines['left'].set_bounds(0,60)
ax_melt_plot.set_ylabel('Melt area\n[%]')
# ax_melt_plot.xticks(rotation=70)
# plt.setp(ax_melt_plot.xaxis.get_majorticklabels(), rotation=70)
import matplotlib.dates as mdates

ax_melt_plot.xaxis.set_major_formatter(mdates.DateFormatter('%b'))

ax_melt_plot.grid(visible=True, which='major', axis='y', alpha=0.33)

plt.savefig('melt.png', transparent=False, bbox_inches='tight', dpi=300)
plt.savefig('melt.svg', transparent=False, bbox_inches='tight', dpi=300)

PROMICE In situ / Point obs

import matplotlib.pyplot as plt
import matplotlib as mpl
import matplotlib.cm as cm
import numpy as np
import pandas as pd
import geopandas as gp
import rasterio as rio
import rasterio.mask
from rasterio.plot import plotting_extent
import cmocean
from mpl_toolkits.axes_grid1.inset_locator import inset_axes

from matplotlib import rc
rc('font', size=10)
rc('text', usetex=False)

fig = plt.figure(1, figsize=(8,8)) # w,h
fig.clf()
# fig.set_tight_layout(True)
import matplotlib.gridspec as gridspec

gs = gridspec.GridSpec(2,2, width_ratios=[1,1], height_ratios=[5,1]) #w,h

ax_map = plt.subplot(gs[0,1])

C_land = "#EAEAEA"
C_ocean = "#D0CFD4"
sub = ['THU_L','KPC_L','UPE_L','SCO_L','KAN_L','NUK_L','TAS_L','QAS_L']

if 'o' not in locals():
    o = gp.read_file('greenland.gpkg')
    
o.plot(color=C_land, ax=ax_map, facecolor='none', zorder=-1)

ice = gp.read_file('ice_edge.gpkg')
ice.boundary.plot(color='k', ax=ax_map, facecolor='None', alpha=0.25, linewidth=0.5, zorder=-1)

ax_map.axis('off')

anom = pd.read_csv('./DVA/PROMICE ablation anomalies (%) (1991-2020 ref).csv',
                  index_col=0, parse_dates=True)
unc = anom.loc['Uncertainty']
anom = anom.loc['2023']

abl = pd.read_csv('./promice_ice_ablation_2023.txt',
                  delim_whitespace=True, index_col=0)
abl = abl.loc[2023]
abl = abl[abl.index.str.contains('|'.join(sub))]
abl.index = [_.split('_')[0] for _ in abl.index]

s = gp.read_file('/home/kdm/data.me/PROMICE/stations.gpkg', index_col=0)\
    .drop(columns=['description','timestamp','begin','end','altitudeMode',
                   'tessellate','visibility','drawOrder','icon',
                   'extrude'])\
    .to_crs('EPSG:3413')

s = s[s['Name'].str.contains('|'.join(sub))]
s['Name'] = [_.split('_')[0] for _ in s['Name']]

s['x'] = s['geometry'].x
s['y'] = s['geometry'].y

s['lon'] = s.to_crs('EPSG:4326')['geometry'].x
s['lat'] = s.to_crs('EPSG:4326')['geometry'].y
s.to_csv('stations.csv')

s = s.merge(anom, left_on='Name', right_index=True)\
     .rename(columns={'2023':'anom'})

s = s.merge(abl, left_on='Name', right_index=True)\
     .rename(columns={2023:'abl'})

s = s.merge(unc, left_on='Name', right_index=True)\
     .rename(columns={'Uncertainty':'unc'})

# ax_map.scatter(s['x'], s['y'], c=s['anom'], s=s['abl']*100, cmap=mpl.cm.RdBu_r)
s['color'] = s['anom'].where(np.abs(s['anom']) > s['unc'])
sc = s.where(~np.isnan(s['color'])).dropna()

# C = sc['color']; C = (C - C.min()) / (C.max()-C.min()); C=(255*C).astype(int)
C = sc['color']; C = ((C + 100)/200 * 255).astype(int)
cmap = mpl.cm.RdBu_r
C = cmap(C)
# C = mpl.cm.RdBu_r(sc['color']/np.max(sc['color'])*255)


im = ax_map.scatter(sc['x'], sc['y'], facecolor=C, s=sc['abl']*100, edgecolor='k', alpha=1, vmin=-100, vmax=100)

sw = s.where(np.isnan(s['color'])).dropna(subset=['Name'])
ax_map.scatter(sw['x'], sw['y'], facecolor='w', s=sw['abl']*100, edgecolor='k')

# ax_map.scatter(-38.4576926,72.579521, facecolor='k')
# summit = gp.GeoDataFrame(geometry=gp.points_from_xy([-38.4576926],[72.579521])).set_crs('EPSG:4326').to_crs('EPSG:3413')
# ax_map.scatter(summit['geometry'].x,summit['geometry'].y, color='k')
# ax_map.annotate('Summit',
#                 xy=(summit['geometry'].x, summit['geometry'].y),
#                 xycoords='data',
#                 xytext=(summit['geometry'].x, summit['geometry'].y-75000),
#                 textcoords='data',
#                 fontsize=12, color='k',
#                 # fontweight='bold',
#                 ha="center", va="center")

def do_text(st, color):
    xoffset = 0 if st['Name'] != 'THU' else -150000
    t0 = ax_map.annotate(st['Name'],
                         xy=(st['x'], st['y']),
                         xycoords='data',
                         xytext=(st['x']+xoffset, st['y']),
                         textcoords='data',
                         fontsize=6, color=color, fontweight='bold',
                         ha="center", va="center")

    plussign = '+' if st["anom"] > 0 else ''
    xoffset = {'KPC':3.0E5,
              'THU':0, # 3.0E5
              'UPE':-3.2E5,
              'SCO':+3.1E5,
              'KAN':-3.2E5,
              'TAS':3.2E5,
              'NUK':-3.2E5,
              'QAS':-3E5}
    yoffset = {'KPC':0,
              'THU':-1.7E5,
              'UPE':0,
              'SCO':0,
              'KAN':0,
              'TAS':0,
              'NUK':0,
              'QAS':0}

    
    t1 = ax_map.annotate(f'{st["abl"]} m \n {plussign}{np.round(st["anom"]).astype(int)} %',
                         xy=(st['x']+xoffset[st['Name']], st['y']+yoffset[st['Name']]),
                         xycoords='data',
                         xytext=(st['x']+xoffset[st['Name']], st['y']+yoffset[st['Name']]),
                         textcoords='data',
                         ha='center', va="center",
                         bbox=dict(boxstyle="round4,pad=0.2",
                                   fc="w", ec="k", lw=2, alpha=0.25),
                         # arrowprops=dict(arrowstyle="->",
                         #                 connectionstyle="arc3"),
                         )


 # ax.text(s['x'].values, s['y'].values, s['Name'].values)
# [ax.text(_['x'].values, _['y'].values, _['Name'].values) for _ in s]
for idx in sc.index:
    st = sc.loc[idx]
    do_text(st, 'k')

for idx in sw.index:
    st = sw.loc[idx]
    do_text(st, 'k')



# REGIONS
region = gp.read_file('/home/kdm/projects/total_mass_balance/tmp/region_interior.gpkg')
region.plot(ax=ax_map, edgecolor='k', facecolor='None', alpha=1)
ax_map.text(-1E5, -1.1E6, 'NO')#, transform=ax_map.TransAxes)
ax_map.text(-2E5, -1.7E6, 'NW')#, transform=ax_map.TransAxes)
ax_map.text(3E5, -1.5E6, 'NE')#, transform=ax_map.TransAxes)
ax_map.text(-1E5, -2.1E6, 'CW')#, transform=ax_map.TransAxes)
ax_map.text(4E5, -2.1E6, 'CE')#, transform=ax_map.TransAxes)
ax_map.text(-1.6E5, -2.7E6, 'SW')#, transform=ax_map.TransAxes)
ax_map.text(1.6E5, -2.5E6, 'SE')#, transform=ax_map.TransAxes)




ax_map_cb = inset_axes(ax_map,
                       width="3%",  # width = 5% of parent_bbox width
                       height="17%",  # height : 50%
                       loc='lower right',
                       bbox_to_anchor=(-0.25, 0.05, 1, 1),
                       axes_kwargs={'yticks':[-100.,100.]},
                       bbox_transform=ax_map.transAxes,
                       borderpad=0)
cb = fig.colorbar(cm.ScalarMappable(norm=None, cmap=cmap),
                  cax=ax_map_cb,
                  label='Ablation Anomaly\n[%]')
cb.ax.set_yticks([0,0.5,1])
cb.ax.set_yticklabels([-100,0,100])
    
plt.savefig('ablation.svg', transparent=False, bbox_inches='tight', dpi=300)
plt.savefig('ablation.png', transparent=False, bbox_inches='tight', dpi=300)
     
# [['Name','anom','abl','unc']]