Commit 3845ad94 authored by Kennedy, Joseph H's avatar Kennedy, Joseph H
Browse files

Add the SMB analysis developed by M.M. Forrester (less data)

parent 5cdd6c3d
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Fri Apr 7 05:53:20 2017
@author: maryforrester
"""
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Tue Mar 28 14:31:32 2017
@author: maryforrester
"""
#%% #%% Modules and directories
import os
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from livvkit.util import elements as EL
describe_basins = """
Histograms of differences between modeled and observed annual surface mass
balance in m w.e. yr-1, at altimetry locations in the Greenland Ice Sheet.
Colors correspond to the eight major drainage basins denoted by Zwally et al.,
2012. High frequencies near zero imply greater model agreement with altimetry
data. Observation data were compiled from the IceBridge dataset (Lewis et al.,
2017) in the accumulation zones. The light blue line highlights 0 difference
between model and observed; values above (below) this line indicate that the
model overestimates (underestimates) SMB in comparison to altimetry
observations. Axes have been normalized for the eight drainage histograms in
order to show differences between the number of available comparison points per
basin. Note that because this plot compares nearest neighbor values without any
spatial interpolation, one coarse model cell may be compared to multiple (if
not many) observed SMB estimates, as the isochrones dataset is of very high
spatial resolution.
"""
describe_all = """
Histogram of differences between modeled and observed annual surface mass
balance in m w.e. yr-1, at altimetry locations in the Greenland Ice Sheet. High
frequencies near zero imply greater model agreement with altimetry data.
Observation data were compiled from the IceBridge dataset (Lewis et al., 2017)
in the accumulation zones. The light blue line highlights 0 difference between
model and observed; values above (below) this line indicate that the model
overestimates (underestimates) SMB in comparison to altimetry observations.
Note that because this plot compares nearest neighbor values without any
spatial interpolation, one coarse model cell may be compared to multiple (if
not many) observed SMB estimates, as the isochrones dataset is of very high
spatial resolution.
"""
def main(args):
img_list = []
plt.style.use = 'default'
IB_df = pd.read_csv(args.iceBridge, index_col=False)
#SMBdf_avg = SMBdf_avg[SMBdf_avg.thk_flag==1]
#Read in all data
#%% model-obs differences: histograms, colored by region
working = IB_df[IB_df.thk_flag==1]
diff_data = pd.DataFrame(data=working.mod_b - working.b,index = working.index,columns=['difference'])
diff_data['zwallyBasin'] = working.zwallyBasin
diff_data['majorbasins'] = np.floor(working.zwallyBasin).astype('str')
#%%
diff_drain = diff_data.groupby('zwallyBasin')
#histplt = diff_drain.hist(bins=20)
colors = {'0.0':'moccasin','1.0':'steelblue','2.0':'darkturquoise',
'3.0':'green','4.0':'lightsalmon','5.0':'mediumpurple',
'6.0':'grey','7.0':'purple','8.0':'firebrick'}
#Exclude cells with unidentified basin and the southernmost drainage (smallest basin)
list_regions = ['1.0','2.0','3.0','4.0','5.0','6.0','7.0','8.0']
pltcrds = [(0,0),(0,1),(1,0),(1,1),(2,0),(2,1),(3,0),(3,1)]
fig, axes = plt.subplots(nrows=4, ncols=2, figsize=(10, 20))
for i in range(0,len(list_regions)):
regname = list_regions[i]
inde = pltcrds[i]
reg_data = diff_data[diff_data.majorbasins==regname]
subs = reg_data['difference'].hist(ax=axes[inde],color=colors[regname])
subs.set_xlabel('Model - observed SMB difference\n(kg m$^{-2}$ a$^{-1}$)')
subs.set_ylabel('Cell frequency')
subs.set_title('Region ' + regname[0],size=12)
subs.set_facecolor('white')
subs.axvline(x=0.0,c="lightskyblue",linewidth=2)
subs.set_xlim(-0.4,0.4)
subs.set_ylim(0,1000)
img_file = os.path.join(args.out, 'IB_diffhist_basins.png')
plt.savefig(img_file)
plt.close()
img_link = os.path.join(os.path.basename(args.out), os.path.basename(img_file))
img_elem = EL.image('IceBridge difference histogram by basin',
' '.join(describe_basins.split()),
img_link)
img_elem['Height'] = args.img_height
img_list.append(img_elem)
#%%
difftest = pd.DataFrame(data = IB_df.mod_b - IB_df.b,index = IB_df.index,columns=['difference'])
fig = plt.figure(figsize = (8,8))
hist1 = difftest['difference'].hist(bins=20)
hist1.set_xlabel('Model - observed SMB difference\n(kg m$^{-2}$ a$^{-1}$)')
hist1.set_ylabel('Cell frequency')
hist1.set_xlim(-0.4,0.4)
hist1.axvline(x=0.0,c="lightskyblue",linewidth=2)
img_file = os.path.join(args.out, 'IB_diffhist_all.png')
plt.savefig(img_file)
plt.close()
img_link = os.path.join(os.path.basename(args.out), os.path.basename(img_file))
img_elem = EL.image('IceBridge difference histogram',
' '.join(describe_all.split()),
img_link)
img_elem['Height'] = args.img_height
img_list.append(img_elem)
return img_list
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Fri Apr 7 05:41:16 2017
@author: maryforrester
"""
#%% #%% Modules and directories
import os
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from livvkit.util import elements as EL
describe = """
Modeled annual surface mass balance versus estimates derived from IceBridge
altimetry data detailed in Lewis et al. (2017). SMB estimates from these were
taken in the accumulation zone (SMB>0), and comparisons are given in kg m-2
yr-1. Colors correspond to the eight major drainage basins denoted by Zwally et
al., 2012 A 1:1 line is drawn in blue. Note that because this plot compares
nearest neighbor values without any spatial interpolation, one coarse model
cell may be compared to multiple (if not many) observed SMB estimates, as the
isochrones dataset is of very high spatial resolution.
"""
def main(args):
img_list = []
plt.style.use = 'default'
IB_df = pd.read_csv(args.iceBridge, index_col=False)
majorbasins = np.floor(IB_df.zwallyBasin).astype('str')
#%% Model vs. observation scatters: 2-panel - accumulation, ablation
accplot = IB_df[IB_df.thk_flag==1] #Remove data outside of modeled ice sheet for comparisons
colors = {'0.0':'moccasin','1.0':'steelblue','2.0':'darkturquoise',
'3.0':'green','4.0':'lightsalmon','5.0':'mediumpurple',
'6.0':'grey','7.0':'purple','8.0':'firebrick'}
#Accumulation
testplot = accplot.plot(kind='scatter',x='b',y='mod_b',s=30,
c=majorbasins.apply(lambda x: colors[x]), \
legend=True,alpha=0.1,figsize=(8,8))
testplot.set_xlabel('Radar SMB estimate (m w.e.)')
testplot.set_ylabel('Modeled SMB (m w.e.)')
#testplot.set_axis_bgcolor('white')
testplot.set_ylim(-.2,1.2)
testplot.set_xlim(-.2,1.2)
testplot.plot([0, 1], [0, 1], transform=testplot.transAxes,zorder=0)
testplot.set_title('Accumulation')
img_file = os.path.join(args.out, 'IB_AccCompare.png')
plt.savefig(img_file)
plt.close()
img_link = os.path.join(os.path.basename(args.out), os.path.basename(img_file))
img_elem = EL.image('Modeled SMB vs. IceBridge',
' '.join(describe.split()),
img_link)
img_elem['Height'] = args.img_height
img_list.append(img_elem)
return img_list
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Thu Apr 6 15:29:38 2017
@author: maryforrester
"""
#%% Modules and directories
import os
import numpy as np
import numpy.ma as ma
import matplotlib.pyplot as plt
import pandas as pd
from mpl_toolkits.basemap import Basemap
from netCDF4 import Dataset
from livvkit.util import elements as EL
describe_spatial = """
Spatial plot of annual surface mass balance along IceBridge transects in the
Greenland accumulation zone. SMB is given in m w.e., and was estimated from
altimetry data.
"""
describe_diff = """
Spatial plot of differences between modeled and observed annual surface mass
balance along IceBridge altimetry transects in the Greenland ice sheet.
Differences are given in m w.e. Modeled SMB values are compared to accumulation
zone data compiled from Lewis et al. (2017). Blue (red) indicates locations at
which the model overestimates (underestimates) surface mass balance when
compared to annual IceBridge estimates.
"""
def main(args):
img_list = []
IB_df = pd.read_csv(args.iceBridge)
#%% Read in model data
nc1 = Dataset(args.smb, mode='r')
nc2 = Dataset(args.latlon, mode='r')
nc3 = Dataset(args.elevation, mode='r')
lats_model = nc2.variables['lat'][:]
lons_model = nc2.variables['lon'][:]
elev_model = nc3.variables['usrf'][:]
SMB_model = nc1.variables['acab'][:] * 1000 #Convert to kg/m2/yr
thk_model = nc1.variables['thk'][:]
nc1.close()
nc2.close()
nc3.close()
#%%
plt.figure(figsize = (12,12))
#This will be the center of our map
lon_0 = lons_model.mean()
lat_0 = lats_model.mean()
#Create the basemap
m = Basemap(width=2000000,height=3000000,
rsphere=(6378137.00,6356752.3142),\
resolution='l',projection='lcc',\
lat_0=lat_0,lon_0=lon_0)
m.drawcoastlines()
m.drawcountries()
m.fillcontinents(color='gainsboro',lake_color='aqua',zorder=1)
m.drawparallels(np.arange(-80., 81., 5.), labels=[1,0,0,0], fontsize=10, color='lightgrey',zorder=1)
m.drawmeridians(np.arange(-180., 181., 10.), labels=[0,0,0,1], fontsize=10, color='lightgrey',zorder=1)
#Format observation data to add
lat_obs = IB_df['Y'].values
lon_obs = IB_df['X'].values
xobs, yobs = m(lon_obs,lat_obs)
smbobs = IB_df['b'].values
obs = m.scatter(xobs,yobs,c=smbobs,marker='o',cmap='viridis',zorder=3,lw=0)
cbar = m.colorbar(obs, location='bottom', pad="5%")
cbar.set_label('Surface mass balance (m w.e. a$^{-1}$)')
img_file = os.path.join(args.out, 'IB_spatial.png')
plt.savefig(img_file)
plt.close()
img_link = os.path.join(os.path.basename(args.out), os.path.basename(img_file))
img_elem = EL.image('IceBridge transects',
' '.join(describe_spatial.split()),
img_link)
img_elem['Height'] = args.img_height
img_list.append(img_elem)
#%% Same thing, but plotting the differences now between altimetry estimate and model
plt.figure(figsize = (12,12))
#Create the basemap
m = Basemap(width=2000000,height=3000000,
rsphere=(6378137.00,6356752.3142),\
resolution='l',projection='lcc',\
lat_0=lat_0,lon_0=lon_0)
m.drawcoastlines()
m.drawcountries()
m.fillcontinents(color='gainsboro',lake_color='aqua',zorder=1)
m.drawparallels(np.arange(-80., 81., 5.), labels=[1,0,0,0], fontsize=10, color='lightgrey',zorder=1)
m.drawmeridians(np.arange(-180., 181., 10.), labels=[0,0,0,1], fontsize=10, color='lightgrey',zorder=1)
#Format observation data to add
lat_obs = IB_df['Y'].values
lon_obs = IB_df['X'].values
xobs, yobs = m(lon_obs,lat_obs)
smbobs = IB_df['mod_b'].values - IB_df['b'].values
obs = m.scatter(xobs,yobs,c=smbobs,marker='o',cmap='RdBu',zorder=3,vmin = -0.2,vmax = 0.2, lw=0)
cbar = m.colorbar(obs, location='bottom', pad="5%")
cbar.set_label('Difference in surface mass balance \n (Model - IceBridge) (m w.e. a$^{-1}$)')
img_file = os.path.join(args.out, 'IB_spatial_difference.png')
plt.savefig(img_file)
plt.close()
img_link = os.path.join(os.path.basename(args.out), os.path.basename(img_file))
img_elem = EL.image('IceBridge transect differences',
' '.join(describe_diff.split()),
img_link)
img_elem['Height'] = args.img_height
img_list.append(img_elem)
return img_list
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Tue Mar 28 14:31:32 2017
@author: maryforrester
"""
#%% #%% Modules and directories
import os
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from livvkit.util import elements as EL
describe_hist = """
Histogram of differences between modeled and observed annual surface mass
balance in kg m-2 yr-1, at all core and firn locations in the Greenland Ice
Sheet. High frequencies near zero imply greater model agreement with core data.
Data are compiled from Cogley et al. (2004), Bales et al. (2009), and the
PROMICE database (Machguth et al., 2016), from both ablation and accumulation
zones. The light blue line highlights 0 difference between model and observed;
values above (below) this line indicate that the model overestimates
(underestimates) SMB in comparison to core observations.
"""
describe_basins = """
Histogram of differences between modeled and observed annual surface mass
balance in kg m-2 yr-1, at core and firn locations in the Greenland Ice Sheet.
Colors correspond to the eight major drainage basins denoted by Zwally et al.,
2012 High frequencies near zero imply greater model agreement with core data.
Data were compiled from Cogley et al. (2004), Bales et al., 2009, and the
PROMICE database (Machguth et al., 2016), from both ablation and accumulation
zones. The light blue line highlights 0 difference between model and observed;
values above (below) this line indicate that the model overestimates
(underestimates) SMB in comparison to core observations. Axes have been
normalized for the eight drainage histograms in order to show differences
between number of available comparison points per basin.
"""
def main(args):
img_list = []
plt.style.use = 'default'
SMBdf_avg = pd.read_csv(args.core)
#SMBdf_avg = SMBdf_avg[SMBdf_avg.thk_flag==1]
#Read in all data
#%% model-obs differences: histograms, colored by region
working = SMBdf_avg[SMBdf_avg.thk_flag==1]
diff_data = pd.DataFrame(data=working.mod_b - working.b,index = working.index,columns=['difference'])
diff_data['zwallyBasin'] = working.zwallyBasin
diff_data['majorbasins'] = np.floor(working.zwallyBasin).astype('str')
#%%
diff_drain = diff_data.groupby('zwallyBasin')
#histplt = diff_drain.hist(bins=20)
colors = {'0.0':'moccasin','1.0':'steelblue','2.0':'darkturquoise',
'3.0':'green','4.0':'lightsalmon','5.0':'mediumpurple',
'6.0':'grey','7.0':'purple','8.0':'firebrick'}
#Exclude cells with unidentified basin and the southernmost drainage (smallest basin)
list_regions = ['1.0','2.0','3.0','4.0','5.0','6.0','7.0','8.0']
pltcrds = [(0,0),(0,1),(1,0),(1,1),(2,0),(2,1),(3,0),(3,1)]
fig, axes = plt.subplots(nrows=4, ncols=2, figsize=(10, 20))
for i in range(0,len(list_regions)):
regname = list_regions[i]
inde = pltcrds[i]
reg_data = diff_data[diff_data.majorbasins==regname]
subs = reg_data['difference'].hist(ax=axes[inde],color=colors[regname])
subs.set_xlabel('Model - observed SMB difference\n(kg m$^{-2}$ a$^{-1}$)')
subs.set_ylabel('Cell frequency')
subs.set_title('Region ' + regname[0],size=12)
subs.set_facecolor('white')
subs.axvline(x=0.0,c="lightskyblue",linewidth=2)
subs.set_xlim(-2000,2000)
subs.set_ylim(0,65)
img_file = os.path.join(args.out, 'core_diffhist_basins.png')
plt.savefig(img_file)
plt.close()
img_link = os.path.join(os.path.basename(args.out), os.path.basename(img_file))
img_elem = EL.image('Core difference histogram by basin',
' '.join(describe_basins.split()),
img_link)
img_elem['Height'] = args.img_height
img_list.append(img_elem)
difftest = pd.DataFrame(data = SMBdf_avg.mod_b - SMBdf_avg.b,index = SMBdf_avg.index,columns=['difference'])
fig = plt.figure(figsize = (8,8))
hist1 = difftest['difference'].hist(bins=30)
hist1.set_xlabel('Model - observed SMB difference\n(kg m$^{-2}$ a$^{-1}$)')
hist1.set_ylabel('Cell frequency')
hist1.set_xlim(-4000,4000)
hist1.axvline(x=0.0,c="lightskyblue",linewidth=2)
img_file = os.path.join(args.out, 'core_diffhist_all.png')
plt.savefig(img_file)
plt.close()
img_link = os.path.join(os.path.basename(args.out), os.path.basename(img_file))
img_elem = EL.image('Core difference histogram',
' '.join(describe_hist.split()),
img_link)
img_elem['Height'] = args.img_height
img_list.append(img_elem)
return img_list[::-1]
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Thu Apr 6 15:29:38 2017
@author: maryforrester
"""
#%% Modules and directories
import os
import numpy as np
import numpy.ma as ma
import matplotlib.pyplot as plt
import pandas as pd
from mpl_toolkits.basemap import Basemap
from netCDF4 import Dataset
from livvkit.util import elements as EL
describe = """
Filled contour of modeled annual surface mass balance of the Greenland ice
sheet, with firn and core field estimates overlaid as filled circles. Data were
compiled from Cogley et al. (2004), Bales et al., 2009, and the PROMICE
database (Machguth et al., 2016), from both ablation and accumulation zones.
"""
def main(args):
img_list = []
SMBdf_avg = pd.read_csv(args.core)
#%% Read in model data
nc1 = Dataset(args.smb, mode='r')
nc2 = Dataset(args.latlon, mode='r')
nc3 = Dataset(args.elevation, mode='r')
lats_model = nc2.variables['lat'][:]
lons_model = nc2.variables['lon'][:]
elev_model = nc3.variables['usrf'][:]
SMB_model = nc1.variables['acab'][:] * 1000 #Convert to kg/m2/yr
thk_model = nc1.variables['thk'][:]
nc1.close()
nc2.close()
nc3.close()
#%%
plt.figure(figsize = (12,12))
#This will be the center of our map
lon_0 = lons_model.mean()
lat_0 = lats_model.mean()
#Create the basemap
m = Basemap(width=2000000,height=3000000,
rsphere=(6378137.00,6356752.3142),\
resolution='l',projection='lcc',\
lat_0=lat_0,lon_0=lon_0)
m.drawcoastlines()
m.drawcountries()
m.fillcontinents(color='gainsboro',lake_color='aqua',zorder=1)
m.drawparallels(np.arange(-80., 81., 5.), labels=[1,0,0,0], fontsize=10, color='lightgrey',zorder=1)
m.drawmeridians(np.arange(-180., 181., 10.), labels=[0,0,0,1], fontsize=10, color='lightgrey',zorder=1)
#Add mask for thk<0.0001
mask = (thk_model.flatten()<0.0001)
SMBflat = SMB_model.flatten()
SMBflat[np.where(mask)] = np.nan
SMBflat.shape = SMB_model.shape
msmb = ma.masked_invalid(SMBflat) #Make sure to convert this to a masked array
#Plotting the pseudocolor
xi, yi = m(lons_model.squeeze(),lats_model.squeeze() )
smb = m.pcolormesh(xi,yi,msmb.squeeze(),cmap='Spectral',zorder=2,vmin=-1500,vmax=1500)
cbar = m.colorbar(smb, location='bottom', pad="5%")
cbar.set_label('Surface mass balance (kg m$^{-2}$ a$^{-1}$)')
#Format observation data to add
lat_obs = SMBdf_avg['Y'].values
lon_obs = SMBdf_avg['X'].values
xobs, yobs = m(lon_obs,lat_obs)
smbobs = SMBdf_avg.b
obs = m.scatter(xobs,yobs,c=smbobs,marker='o',edgecolor='black',
cmap='Spectral',zorder=3,vmin=-1500,vmax=1500)
img_file = os.path.join(args.out, 'core_spatial.png')
plt.savefig(img_file)
plt.close()
img_link = os.path.join(os.path.basename(args.out), os.path.basename(img_file))
img_elem = EL.image('Modeled SMB with field overlay',
' '.join(describe.split()),
img_link)
img_elem['Height'] = args.img_height
img_list.append(img_elem)
return img_list
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Tue Mar 28 14:28:29 2017
@author: maryforrester
"""
#%% #%% Modules and directories
import os
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from livvkit.util import elements as EL
describe_felev = """
Field estimates of annual surface mass balance in kg m-2 yr-1 as a function of
observed elevation in meters. Data were compiled from Cogley et al. (2004),
Bales et al. (2009), and the PROMICE database (Machguth et al., 2016), from
both ablation and accumulation zones. Colors correspond to the eight major
drainage basins denoted by Zwally et al., 2012 Size of the point represents the
number of years in the field record, with larger points containing
comparatively more temporal information than smaller points. Each point
represents an annual surface mass balance estimate averaged across at least one
year of data.