- Introduction
- Mass loss from citet:mankoff_2021
- Discharge
- Greenland outline
- Bare ice area
- Albedo
- Melt
- PROMICE In situ / Point obs
GISS contributions to the Arctic Report Card
# 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- ~ 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
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
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
- Using v90 from https://doi.org/10.22008/promice/data/ice_discharge/d/v02
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
df['2013-01-01':].resample('1D').interpolate().mean()
discharge 500.629915 err 46.894040 dtype: float64
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
See ./figs_tmp sub-folder for graphics
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
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')
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>
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.gpkgimport 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)
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"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)
ls TMgrass -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"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)
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']]
