This IPython notebook produces the error x distance-size ratio figures for the section "Evaluation of the accuracy" of the article.
The data was calculated by the tesseroid_vs_spherical_shell.ipynb notebook and saved to CSV files in the data directory of the main repository.
We'll use the pandas library to load and manipulate the data. Plots will be made with matplotlib.
In [1]:
%matplotlib inline
from __future__ import division
from IPython.display import set_matplotlib_formats
set_matplotlib_formats('svg')
import pandas as pd
import numpy as np
from matplotlib import pyplot as plt
In [2]:
plt.rcParams['axes.labelsize'] = 9.0 # fontsize of the x any y labels
plt.rcParams['xtick.labelsize'] = 9.0 # fontsize of the tick labels
plt.rcParams['ytick.labelsize'] = 9.0 # fontsize of the tick labels
plt.rcParams['legend.fontsize'] = 9.0
plt.rcParams['font.family'] = 'serif'
plt.rcParams['font.serif'] = 'Computer Modern Roman'
plt.rcParams['text.usetex'] = True # use latex for all text handling
plt.rcParams['text.color'] = '3a3a3a'
plt.rcParams['figure.facecolor'] = 'white'
plt.rcParams['axes.linewidth'] = 1
plt.rcParams['axes.edgecolor'] = '3a3a3a'
plt.rcParams['axes.facecolor'] = 'white'
plt.rcParams['lines.linewidth'] = 1
plt.rcParams['lines.markersize'] = 4
plt.rcParams['xtick.major.size'] = 2
plt.rcParams['ytick.major.size'] = 2
And load the computed shell effect per height so that we can calculate the relative difference.
In [3]:
shell = pd.read_csv('../data/shell-per-height.csv', index_col=0)
In [4]:
shell
Out[4]:
Load the computed differences between the shell and tesseroid effects for each experiment.
We'll use the groupby method of the pandas DataFrame to compute the maximum absolute difference per distance-size ratio and keep that.
In [5]:
def load_max_diff(fname):
data = pd.read_csv(fname)
return data.groupby(['size', 'ratio']).max()
In [6]:
pole = load_max_diff('../data/difference-size-1-height-2000-pole.csv')
equator = load_max_diff('../data/difference-size-1-height-2000-equator.csv')
goce = load_max_diff('../data/difference-size-1-height-260000-pole.csv')
big = load_max_diff('../data/difference-size-30-height-2000-pole.csv')
This is what the data looks like after taking the maximum difference per ratio.
In [7]:
pole.head()
Out[7]:
Make a plot of the maximum relative difference as a function of the distance-size ratio used in the computations.
In [8]:
ratio = pole.index.levels[1] # Get the unique values of the distance-size ratio used
In [9]:
plotargs = dict(color='#3a3a3a', markeredgewidth=0.5, markeredgecolor='#3a3a3a')
styles = "-s -^ -v -o".split()
fig, subplots = plt.subplots(3, 1, figsize=(3.33, 5), sharex='col')
axes = subplots.ravel()
fields = ['pot', 'gz', 'gzz']
titles = [r'V', r'gz', r'gzz']
subfigure = ['(a)', '(b)', '(c)']
for ax, f, title, sub in zip(axes, fields, titles, subfigure):
ax.text(-0.21, 0.9, sub, fontsize=12, fontdict={'weight': 'bold'},
transform=ax.transAxes)
ax.text(0.5, 0.9, title, fontsize=11,
horizontalalignment='center', verticalalignment='center',
bbox={'facecolor': 'w',
'edgecolor': '#9b9b9b',
'linewidth': 0.5, 'pad': 8},
transform=ax.transAxes)
shell_low = np.abs(shell[shell.height == 2000][f].values)
shell_high = np.abs(shell[shell.height == 260000][f].values)
ax.plot(ratio, 100*pole[f]/shell_low, styles[0], label='pole',
**plotargs)
ax.plot(ratio, 100*equator[f]/shell_low, styles[1], label='equator',
**plotargs)
ax.plot(ratio, 100*goce[f]/shell_high, styles[2], label='260 km',
**plotargs)
ax.plot(ratio, 100*big[f]/shell_low, styles[3], label=r'$30^\circ$ size',
**plotargs)
ax.hlines(0.1, ratio.min(), ratio.max(), colors=['#3a3a3a'], linewidth=1.5)
ax.set_xlim(ratio.min(), ratio.max())
ax.set_yscale('log')
ax.set_xticks(range(11))
ax.set_yticks(ax.get_yticks()[2:-2])
ax.set_ylabel('Difference (\\%)')
ax.grid(True, linewidth=0.5, color='#aeaeae')
ax.set_axisbelow(True)
ax.minorticks_off()
ax = axes[-1]
ax.set_xlabel('Distance-size ratio')
ax.legend(borderpad=0.5, numpoints=1, bbox_to_anchor=(1, 1),
fancybox=True, shadow=False, fontsize=9, )
plt.tight_layout(pad=0.3, h_pad=0, w_pad=0)
plt.subplots_adjust(hspace=0, wspace=0)
plt.savefig('../figs/distance-size-curves.eps') # Save the figure to an EPS
Load the data for the gzz difference per height. I'll use the pre-loaded spherical shell data to calculate the relative difference.
In [10]:
gzz = pd.read_csv('../data/difference-gzz-per-height-size-1-pole.csv',
index_col=0).groupby(['height', 'ratio']).max()
This is what the first few lines of the data look like.
In [11]:
gzz.head()
Out[11]:
In [12]:
heights = gzz.index.levels[0].values # Get the unique values of the computation height
In [13]:
heights
Out[13]:
Now we'll make the difference x ratio plot with a curve for each height.
In [14]:
color = '#3a3a3a'
styles = "-s -^ -v -o -*".split()
fig = plt.figure(figsize=(3.33, 3))
ax = plt.subplot(111)
for h, sty in zip(heights, styles):
shell_value = np.abs(shell[shell.height == h]['gzz'].values)
diff = 100*gzz.loc[int(h)]/shell_value
markersize = 4
if sty[-1] == '*':
markersize = 7
ax.plot(ratio, diff, sty, label='{:.0f} km'.format(h/1000), color=color,
markeredgewidth=0.5, markeredgecolor=color, markersize=markersize)
ax.hlines(0.1, ratio.min(), ratio.max(), colors=['#3a3a3a'], linewidth=1.5)
ax.set_xlim(ratio.min(), ratio.max())
ax.set_yscale('log')
ax.set_xticks(range(11))
ax.set_yticks(ax.get_yticks())
ax.minorticks_off()
ax.grid(True, linewidth=0.5, color='#aeaeae')
ax.set_axisbelow(True)
ax.set_ylabel('Difference (\\%)')
ax.set_xlabel('Distance-size ratio')
ax.legend(borderpad=0.5, numpoints=1, bbox_to_anchor=(1, 1),
fancybox=True, shadow=False, fontsize=9, )
plt.tight_layout(pad=0.25, h_pad=0, w_pad=0)
plt.subplots_adjust(hspace=0, wspace=0)
plt.savefig('../figs/gzz-with-height.eps') # Save the figure to an EPS