Commit 0e09c069 authored by Kennedy, Joseph H's avatar Kennedy, Joseph H
Browse files

Make basin #ing match Zwally, make units consistent and SI, some PEP8 cleanup

parent d47dc988
......@@ -6,14 +6,6 @@ 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
......@@ -24,32 +16,32 @@ 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.
balance, at altimetry locations on 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
2017) in the accumulation zones. The light blue line highlights zero 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
not many) observed SMB estimates, as the isochrones dataset has 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
balance, 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
in the accumulation zones. The light blue line highlights zero 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
not many) observed SMB estimates, as the isochrones dataset has very high
spatial resolution.
"""
......@@ -77,14 +69,14 @@ def main(args):
#histplt = diff_drain.hist(bins=20)
colors = {'0.0':'moccasin','1.0':'steelblue','2.0':'darkturquoise',
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)]
zwally_regions = [2, 3, 4, 1, 5, 6, 7, 8]
pltcrds = [(0,1), (1,0), (1,1), (0,0), (2,0), (2,1), (3,0), (3,1)]
fig, axes = plt.subplots(nrows=4, ncols=2, figsize=(10, 20))
......@@ -95,12 +87,13 @@ def main(args):
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_title('Region {}'.format(zwally_regions[i]),size=12)
subs.set_facecolor('white')
subs.axvline(x=0.0,c="lightskyblue",linewidth=2)
subs.set_xlim(-0.4,0.4)
subs.set_xlim(-400,400)
subs.set_ylim(0,1000)
plt.tight_layout()
img_file = os.path.join(args.out, 'IB_diffhist_basins.png')
plt.savefig(img_file)
plt.close()
......@@ -119,7 +112,7 @@ def main(args):
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.set_xlim(-400,400)
hist1.axvline(x=0.0,c="lightskyblue",linewidth=2)
img_file = os.path.join(args.out, 'IB_diffhist_all.png')
......
......@@ -47,11 +47,11 @@ def main(args):
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_xlabel('Radar SMB estimate (kg m$^{-2}$ a$^{-1}$)')
testplot.set_ylabel('Modeled SMB (kg m$^{-2}$ a$^{-1}$)')
#testplot.set_axis_bgcolor('white')
testplot.set_ylim(-.2,1.2)
testplot.set_xlim(-.2,1.2)
testplot.set_ylim(-200,1200)
testplot.set_xlim(-200,1200)
testplot.plot([0, 1], [0, 1], transform=testplot.transAxes,zorder=0)
testplot.set_title('Accumulation')
......
......@@ -19,14 +19,14 @@ 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
Greenland accumulation zone. SMB is given in kg m^-2 a^-1, 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
Differences are given in m^-2 a^-1 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.
......@@ -80,7 +80,7 @@ def main(args):
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}$)')
cbar.set_label('Surface mass balance (kg m$^-2$ a$^{-1}$)')
img_file = os.path.join(args.out, 'IB_spatial.png')
plt.savefig(img_file)
......@@ -112,10 +112,10 @@ def main(args):
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)
obs = m.scatter(xobs,yobs,c=smbobs,marker='o',cmap='RdBu',zorder=3,vmin = -200,vmax = 200, 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}$)')
cbar.set_label('Difference in surface mass balance \n (Model - IceBridge) (kg m$^-2$ a$^{-1}$)')
img_file = os.path.join(args.out, 'IB_spatial_difference.png')
plt.savefig(img_file)
......
......@@ -71,8 +71,8 @@ def main(args):
#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)]
zwally_regions = [2, 3, 4, 1, 5, 6, 7, 8]
pltcrds = [(0,1), (1,0), (1,1), (0,0), (2,0), (2,1), (3,0), (3,1)]
fig, axes = plt.subplots(nrows=4, ncols=2, figsize=(10, 20))
......@@ -83,12 +83,13 @@ def main(args):
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_title('Region {}'.format(zwally_regions[i]),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)
plt.tight_layout()
img_file = os.path.join(args.out, 'core_diffhist_basins.png')
plt.savefig(img_file)
plt.close()
......
......@@ -38,7 +38,7 @@ Annual surface mass balance in kg m-2 yr-1 as a function of observed elevation
in meters at four different transect locations in the ablation zone. Red dots
are modeled SMB and elevation, while blue dots are observed SMB and elevation
from firn and core field estimates. Observation data were compiled from the
PROMICE database (Machguth et al., 2009), from ablation zones.
PROMICE database (Machguth et al., 2016).
"""
describe_facc = """
......@@ -85,10 +85,11 @@ def main(args):
#testplot = SMBdf_avg.plot(kind='scatter',x='Z',y='b',s=size*20,c='steelblue')
testplot = SMBdf_avg.plot(kind='scatter',x='Z',y='b',s=size*10+1,c=majorbasins.apply(lambda x: colors[x]),legend=True,alpha=0.6)
testplot.set_ylabel('Modeled SMB (kg m$^{-2}$ a$^{-1}$) ')
testplot.set_ylabel('Field SMB estimate (kg m$^{-2}$ a$^{-1}$) ')
testplot.set_xlabel('Elevation (m)')
#testplot.set_facecolor('white')
testplot.axhline(y=0.0,c="darkgrey",linewidth=1,zorder=0)
testplot.set_ylim(-5500,1500)
img_file = os.path.join(args.out, 'core_smbelev_obs.png')
plt.savefig(img_file)
......@@ -103,10 +104,11 @@ def main(args):
#Model with added transparency
testplot2 = SMBdf_avg.plot(kind='scatter',x='mod_Z',y='mod_b',s=30,c=majorbasins.apply(lambda x: colors[x]),legend=True,alpha=0.6)
testplot2.set_ylabel('Field SMB estimate (kg m$^{-2}$ a$^{-1}$)')
testplot2.set_ylabel('Modeled SMB (kg m$^{-2}$ a$^{-1}$)')
testplot2.set_xlabel('Elevation (m)')
#testplot2.set_facecolor('white')
testplot2.axhline(y=0.0,c="darkgrey",linewidth=1,zorder=0)
testplot2.set_ylim(-5500, 1500)
img_file = os.path.join(args.out, 'core_smbelev_mod.png')
plt.savefig(img_file)
......@@ -122,7 +124,7 @@ def main(args):
#%% Elevation vs. SMB scatters: 4-panel - 4 glacier transects
working = SMBdf_avg[SMBdf_avg.thk_flag==1]
glacier_list = ['414','180','215','454']
glac_longnames = ['Qammanarssup Sermia','Nioghalvjerds-fjorden','Storstommen','K-transect']
glac_longnames = ['Qamanârssêp Sermia','Nioghalfvjerdsfjorden','Storstmmen','K-transect']
#Function to grab indices in SMBdf (observed and modeled) given a particular glacier
#These are the same transects used in Noel et al., 2016
......@@ -133,7 +135,7 @@ def main(args):
glac_name = glacier_list[i]
glac_data = working[working.glacier_ID==glac_name]
glacplot.scatter(x=glac_data.Z,y=glac_data.b,color='slateblue',label='Observed')
glacplot.scatter(x=glac_data.mod_Z,y=glac_data.mod_b,color='lightcoral',label='Modeled')
glacplot.scatter(x=glac_data.mod_Z,y=glac_data.mod_b,marker='d',color='lightcoral',label='Modeled')
if i==0:
glacplot.legend(loc='upper left')
glacplot.set_facecolor('white')
......@@ -151,7 +153,7 @@ def main(args):
plt.close()
img_link = os.path.join(os.path.basename(args.out), os.path.basename(img_file))
img_elem = EL.image('Transect ablation vs. obs. elevation',
img_elem = EL.image('Ablation transect SMB vs. obs. elevation',
' '.join(describe_transelev.split()),
img_link)
img_elem['Height'] = args.img_height
......
......@@ -9,80 +9,45 @@ Created on Thu Apr 6 15:29:38 2017
#%% Modules and directories
import os
import numpy as np
import pandas as pd
import numpy.ma as ma
import matplotlib.pyplot as plt
import pandas as pd
from mpl_toolkits.basemap import Basemap
from matplotlib import colors as c
from netCDF4 import Dataset
import argparse
this_script_loc = os.path.dirname(os.path.abspath(__file__))
base_path = os.path.dirname(this_script_loc)
def rel_base_path(file):
file = os.path.relpath(file, base_path)
return file
def parse_args(args=None):
parser = argparse.ArgumentParser(description=__doc__,
formatter_class=argparse.ArgumentDefaultsHelpFormatter)
parser.add_argument('--elevation',
default = '../Greenland_5km_v1.1_SacksRev_c110629.nc',
type=rel_base_path,
help='File name of elevation model.')
parser.add_argument('--latlon',
default = '../Greenland_5km_v1.1_SacksRev_c110629.nc',
type=rel_base_path,
help='File name of latlon model.')
parser.add_argument('--smb',
default = '../b.e10.BG20TRCN.f09_g16.002_ANN_196001_200512_climo.nc',
type=rel_base_path,
help='File name of smb model.')
parser.add_argument('--core',
default = '../outputs/SMB_Obs_Model.csv',
type=rel_base_path,
help='The core output file.')
parser.add_argument('--zwally',
default = '../outputs/model_zwally_basins.csv',
type=rel_base_path,
help='The zwally output file.')
parser.add_argument('--iceBridge',
default = '../outputs/IceBridge_modelXY.csv',
type=rel_base_path,
help='The Ice Bridge output file.')
parser.add_argument('--out',
default = '../outputs',
type=rel_base_path,
help='The output directory for the figure.')
return parser.parse_args(args)
from livvkit.util import elements as EL
describe_metadata = """
Map describing the datasets used in model surface mass balance validation,
including core/firn locations, altimetry transects, and drainage basin
delineations in the Greenland Ice Sheet. Drainage basins were denoted by Zwally
et al. (2012), and their colors correspond to those used in other histograms
and scatterplots that compare modeled and observed SMB. Altimetry data (white
lines) were obtained from the IceBridge project (Lewis et al., 2017) and SMB
values along transects represent the annual average from the temporal domain of
each transect (the oldest resolvable internal reflecting horizon varies for
each airborne radar location, but can reach as far back as the early 1700s).
Core/firn estimates are colored based on their source; dark blue circles are
locations compiled by Cogley et al. (2004) and updated by Bales et al. (2009)
PARCA cores, while yellow triangles are estimates supplied by PROMICE (Machguth
et al., 2016). Core data points are sized by the length of their temporal
record, with larger points indicating annual estimates were taken from a
greater number of yearly SMB values.
"""
def main(args):
img_list = []
IB_df = pd.read_csv(args.iceBridge)
SMBdf_avg = pd.read_csv(args.core)
zwally_data = pd.read_csv(args.zwally)
#os.chdir('/Users/maryforrester/Desktop/testing')
#print(os.getcwd())
SMBdf_avg = pd.read_csv(os.path.join(base_path, args.core))
zwally_data = pd.read_csv(os.path.join(base_path,args.zwally))
#%% Read in model data
nc1 = Dataset(os.path.join(base_path, args.smb), mode='r')
nc2 = Dataset(os.path.join(base_path, args.latlon), mode='r')
nc3 = Dataset(os.path.join(base_path, args.elevation), mode='r')
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'][:]
......@@ -95,7 +60,6 @@ def main(args):
nc2.close()
nc3.close()
#%%
plt.figure(figsize = (12,12))
forcolors = {'0.0':'moccasin','1.0':'steelblue','2.0':'darkturquoise',
......@@ -134,37 +98,52 @@ def main(args):
#cbar = m.colorbar(smb, location='bottom', pad="5%")
#cbar.set_label('Zwally drainage basin')
basin_labels = {'1': (-48, 80.5),
'2': (-29, 76),
'3': (-32, 70),
'4': (-40.5, 66),
'5': (-47.5, 62),
'6': (-49, 67.5),
'7': (-50, 71),
'8': (-55, 75),}
for lbl, loc in basin_labels.items():
xi, yi = m(loc[0], loc[1])
plt.text(xi, yi, lbl, fontsize=16, fontweight='bold', color='k')
#Adding in firn/core measurement locations. Size by the number of years in the record
lat_obs = SMBdf_avg['Y'].values
lon_obs = SMBdf_avg['X'].values
xobs, yobs = m(lon_obs,lat_obs)
forsize = SMBdf_avg.nyears
forsize.loc[forsize>50] = 50
forsize.loc[forsize<5] = 5
obs1 = m.scatter(xobs,yobs,marker='o',edgecolor='black',s = 4*forsize,zorder=4)
SMBdf_avg.loc[SMBdf_avg['nyears'] > 50] = 50
SMBdf_avg.loc[SMBdf_avg['nyears'] < 5] = 5
obs1 = m.scatter(xobs,yobs,marker='o',edgecolor='black',s = 4*SMBdf_avg['nyears'],zorder=4)
#Color the ablation zone (PROMICE) measurements yellow
SMBdf_promice = SMBdf_avg[SMBdf_avg.source=='promice']
SMBdf_promice = SMBdf_avg[SMBdf_avg.source=='promice'].copy()
lat_obs = SMBdf_promice['Y'].values
lon_obs = SMBdf_promice['X'].values
xobs, yobs = m(lon_obs,lat_obs)
forsize = SMBdf_promice.nyears
forsize.loc[forsize>50] = 50
forsize.loc[forsize<5] = 5
obs2 = m.scatter(xobs,yobs,marker='^',color='yellow',edgecolor='black',s = 4*forsize,zorder=4)
SMBdf_promice.loc[SMBdf_promice['nyears'] > 50] = 50
SMBdf_promice.loc[SMBdf_promice['nyears'] < 5] = 5
obs1 = m.scatter(xobs,yobs,marker='^',color='yellow',edgecolor='black',s = 4*SMBdf_promice['nyears'],zorder=4)
#Add IceBridge transects
IB_df = pd.read_csv(os.path.join(base_path, args.iceBridge))
lat_obs = IB_df['Y'].values
lon_obs = IB_df['X'].values
print(len(lon_obs))
#print(len(lon_obs))
xobs, yobs = m(lon_obs,lat_obs)
obs3 = m.scatter(xobs,yobs,marker='o',lw=0,zorder=3,s=3,color='white')
plt.savefig(os.path.join(base_path, args.out, 'metadata_spatial.png'))
#plt.show()
#plt.close()
img_file = os.path.join(args.out, 'plot_meta.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('Data locations',
' '.join(describe_metadata.split()),
img_link)
img_elem['Height'] = args.img_height
img_list.append(img_elem)
if __name__ == '__main__':
main(parse_args())
return img_list
......@@ -91,7 +91,7 @@ def main(args):
for i in range(0,len(IB_df)):
avgsmb = np.nanmean(IB_df.iloc[i,2:19])
IB_df.loc[i,'b'] = avgsmb
IB_df.loc[i,'b'] = avgsmb * 1000 # convert m w.e. a^-1 to kg m^2 a^-1
#%% Read in model data
......@@ -102,7 +102,7 @@ def main(args):
lats_model = nc2.variables['lat'][:]
lons_model = nc2.variables['lon'][:]
elev_model = nc3.variables['usrf'][:]
SMB_model = nc1.variables['acab'][:]
SMB_model = nc1.variables['acab'][:] * 1000 # convert m w.e. a^-1 to kg m^2 a^-1
thk_model = nc1.variables['thk'][:]
nc1.close()
......
......@@ -47,12 +47,8 @@ et al. (2012).
from __future__ import absolute_import, division, print_function, unicode_literals
import os
import sys
import argparse
import numpy as np
from netCDF4 import Dataset
import livvkit
from livvkit.util import elements as EL
......@@ -72,6 +68,8 @@ with FN.temp_sys_path(os.path.dirname(__file__)):
import smb.plot_IB_hist as IB_hist
import smb.plot_IB_spatial as IB_spatial
import smb.plot_IB_scatter as IB_scatter
import smb.plot_metadata_spatial as meta
base_path = os.path.dirname(os.path.abspath(__file__))
......@@ -98,7 +96,7 @@ def run(name, config):
args = parse_args(config_arg_list)
#PREPROCESS
# PREPROCESS
if 'zwally' in config['preprocess']:
zwally.main(zwally.parse_args())
if 'icebridge' in config['preprocess']:
......@@ -106,9 +104,11 @@ def run(name, config):
if 'core' in config['preprocess']:
core.main(core.parse_args())
#PLOTS
# PLOTS
img_list = []
img_list.extend(meta.main(args))
img_list.extend(c_spatial.main(args))
img_list.extend(IB_spatial.main(args))
......@@ -121,11 +121,14 @@ def run(name, config):
img_list.extend(transects[:3])
element_list = []
element_list.append(EL.gallery('Figures', img_list))
element_list = [EL.gallery('Figures', img_list)]
ref_bib = utils.bib2html(config['References'])
element_list.append(EL.html('<div class="references"><h3>References</h3> If you use this LIVVkit extension for any part of your modeling or analyses, please cite:' + ref_bib + '</div>'))
element_list.append(EL.html(' '.join(['<div class="references"><h3>References</h3>',
'If you use this LIVVkit extension for any part of your',
'modeling or analyses, please cite:' + ref_bib + '</div>',
])
))
return EL.page(name, __doc__, element_list)
......@@ -159,27 +162,27 @@ def rel_base_path(file):
def parse_args(args=None):
parser = argparse.ArgumentParser(description=__doc__,
formatter_class=argparse.ArgumentDefaultsHelpFormatter)
formatter_class=argparse.ArgumentDefaultsHelpFormatter)
parser.add_argument('--elevation',
default = 'data/Greenland_5km_v1.1_SacksRev_c110629.nc',
type=rel_base_path,
help='File name of elevation model.')
default='data/Greenland_5km_v1.1_SacksRev_c110629.nc',
type=rel_base_path,
help='File name of elevation model.')
parser.add_argument('--latlon',
default = 'data/Greenland_5km_v1.1_SacksRev_c110629.nc',
type=rel_base_path,
help='File name of latlon model.')
default='data/Greenland_5km_v1.1_SacksRev_c110629.nc',
type=rel_base_path,
help='File name of latlon model.')
parser.add_argument('--smb',
default = 'data/b.e10.BG20TRCN.f09_g16.002_ANN_196001_200512_climo.nc',
type=rel_base_path,
help='File name of smb model.')
default='data/b.e10.BG20TRCN.f09_g16.002_ANN_196001_200512_climo.nc',
type=rel_base_path,
help='File name of smb model.')
parser.add_argument('--out',
default = 'outputs',
type=rel_base_path,
help='The output directory for the figures.')
default='outputs',
type=rel_base_path,
help='The output directory for the figures.')
args, _ = parser.parse_known_args(args)
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment