Code for interpreting the hydrogen content across profiles in NAMs using FTIR spectra.
For those new to python, I recommend installing it through anaconda. pynams works with either python 2.7 or 3.
Once you have python, download pynams and install it by going into the Pynams folder on your command line and typing 'python setup.py install'.
After that you should be able to follow along with the code below anywhere you can run a python script. I like Jupyter notebooks (this format) for quick things that I want to show people and spyder for more involved projects. Both come with anaconda.
The first and seconds lines make plots show up nicely in a jupyter notebook.
The third line will let us use some special plotting commands for fiddling with figures later on.
The fourth line makes python 2.7 code compatible with python 3.
In [1]:
%matplotlib inline
%config InlineBackend.figure_formats=["svg"]
import matplotlib.pyplot as plt
from __future__ import print_function, division
In [2]:
import pynams
from pynams import Spectrum
### Here I'm pointing to example data within the pynams folder
import os # package to figure out where pynams is on your computer
FTIR_folder = os.path.join(pynams.__path__[0], 'example_FTIR_spectra\\')
spectrum_augite1 = Spectrum(fname='augite1', folder=FTIR_folder, thickness_microns=876.4)
### Here is an example of how to pass in data from a personal computer
### Mac people note there is no C at the beginning of your path
#FTIR_file_location = 'C:\\Users\\Ferriss\\Documents\\Code\\Pynams\\example_FTIR_spectra\\'
#spectrum_augite1 = Spectrum(fname='augite1', folder=FTIR_file_location, thickness_microns=876.4)
In [3]:
fig, ax = spectrum_augite1.plot_spectrum()
Next we'll plot two spectra on top of each other.
If you don't like the color scheme, you can either
The keyword offset moves your spectrum up and down by the specified amount.
In [4]:
# Making second spectrum called spectrum2 the same way as spectrum above
spectrum_augite2 = Spectrum(fname = 'augite2', thickness_microns = 876.4, folder=FTIR_folder)
# A larger plot with two spectra and a legend
fig, ax = spectrum_augite1.plot_spectrum(label='initial')
spectrum_augite2.plot_spectrum(axes=ax, label='3 hrs\n800$\degree$C', offset=0.5)
fig.set_size_inches(5, 5)
# the y-axis limits should be set to whatever is appropriate for your spectra
ax.set_ylim(0, 4.)
ax.set_title('Ferriss augite dehydration spectra', fontsize=14)
hleg = ax.legend(loc=1, fontsize=14)
# change font sizes
# this line is why we imported matplotlib.pyplot at the beginning
plt.tick_params(axis='both', which='major', labelsize=14)
ax.xaxis.label.set_size(14)
ax.yaxis.label.set_size(14)
This example uses data from Ferriss et al. 2016
In [5]:
from pynams import Sample
Kunlun_diopside_K4 = Sample(length_a_microns=7000, length_b_microns=[2185, 2190, 2188, 2185, 2188],
length_c_microns=[1546, 1551, 1536, 1548, 1548])
Kunlun_folder = os.path.join(pynams.__path__[0], 'example_FTIR_spectra\\KunlunDiopside\\')
spectrum_diopside = Spectrum(fname='K4q_cdb01', folder=Kunlun_folder, sample=Kunlun_diopside_K4, raypath='b')
print(spectrum_diopside.thickness_microns)
In [6]:
spectrum_olivine = Spectrum(folder=FTIR_folder, fname='olivine1')
spectrum_olivine.get_thickness_from_SiO()
Out[6]:
Lemaire et al. 2004 Figure 1 shows typical O-H stretches and 2nd order Si-O overtones for FTIR spectra for olivine with light polarized in different directions. spectrum.orientation(), produces a plot zoomed in on the relevant Si-O peaks with labeled lines for the major peaks shown in the Lemaire figure. This example is probably oriented with the electric vector || c.
In [7]:
fig = spectrum_olivine.orientation()
fig.set_size_inches(8, 8)
In [8]:
fig, ax = spectrum_augite1.plot_showbaseline()
The wavenumber range of the baseline is set by the keywords wn_high and wn_low, which default to 3700 and 3200 cm-1.
The keyword force_though_wn will force a fit though quadratic line through the curve at that wavenumber or list of wavenumbers.
Or tell it how much to shift the quadrtic curve away from being a line using the curvature keyword.
In [9]:
spectrum_augite1.make_baseline(force_through_wn=3500, wn_low=3400, show_fit_values=True)
fig, ax = spectrum_augite1.plot_showbaseline()
In [10]:
spectrum_augite1.make_baseline(curvature=0.3, wn_low=3000, show_fit_values=True, show_plot=False)
fig, ax = spectrum_augite1.plot_showbaseline(offset=1)
ax.set_ylim(0, 5)
Out[10]:
Uses the scipy 1D interpolations and a default spline kind='cubic'.
In [11]:
spectrum_augite1.make_baseline(linetype='spline')
fig, ax = spectrum_augite1.plot_showbaseline()
In [12]:
spectrum_augite1.make_baseline(force_through_wn=[3500, 3300], wn_low=3200, polynomial_order=3)
fig, ax = spectrum_augite1.plot_showbaseline()
In [13]:
my_1st_baseline = spectrum_augite1.make_baseline(force_through_wn=[3500, 3300], wn_low=3200, polynomial_order=3)
my_2nd_baseline = spectrum_augite1.make_baseline()
fig, ax = spectrum_augite1.plot_showbaseline()
spectrum_augite1.base_abs = my_1st_baseline
fig, ax = spectrum_augite1.plot_showbaseline()
In [14]:
spectrum_noisy_olivine = Spectrum(folder=FTIR_folder, fname='olivine4', thickness_microns=360.)
spectrum_noisy_olivine.make_baseline(wn_high=3655, abs_smear_high=10, wn_low=3350, abs_smear_low=10, curvature=0.03)
fig, ax = spectrum_noisy_olivine.plot_showbaseline()
fig.set_size_inches(6, 6)
ax.set_ylim(0, 0.3)
Out[14]:
In [15]:
spectrum_augite1.save_baseline();
In [16]:
spectrum_augite1.get_baseline();
In [17]:
fig, ax = spectrum_augite1.plot_subtractbaseline();
In [18]:
spectrum_augite1.get_baseline();
fig, ax, area = spectrum_augite1.make_area(show_plot=True, wn_low=3400, wn_high=3500)
Defaults are phase='cpx', calibration='Bell', and scale_water (how much to multiply the resulting water concentration by to account for areas in spectra in different direction that weren't measured) of 3.
These results are not really to be trusted, but will give you a very rough estimate. See, e.g., Rossman 2006 review.
The reported errors are based only on propagated error from the absorption coefficients and do not include the large (30%?) uncertainties associated with baseline choice and the potentially even larger errors associated with using unpolarized radiation (Withers 2013).
In [19]:
spectrum_augite1.water(printout=True, phase='cpx', calibration='Bell')
spectrum_noisy_olivine.water(printout=True, phase='ol', calibration='Withers', scale_water=2.1);
In [20]:
from pynams import Profile
profile = Profile(name='My first FTIR profile',
fnames=['olivine1', 'olivine2', 'olivine3', 'olivine4'],
folder=FTIR_folder,
thicknesses_microns=300.,
positions_microns=[10., 20., 30., 40.])
In [21]:
profile.get_thicknesses_from_SiO()
from pynams import styles
fig, ax, ax_water = styles.plot_area_profile_outline(show_water_ppm=False)
ax.set_ylabel('thickness $\mu$m')
ax.plot(profile.positions_microns, profile.thicknesses_microns, '^k')
ax.set_ylim(300, 400)
ax.set_xlim(0, 50)
ax.set_title('Estimated thickness across profile');
In [22]:
# make quadratic baselines. Keywords are the same as in spectrum.make_baseline() above
profile.make_baselines(wn_low=3100, wn_high=3600., force_through_wn=3300)
#profile.plot_showbaselines() # uncomment this out if you want to look at the baselines
# This happens automatically during plotting if profiles.areas hasn't been set yet
profile.make_areas()
# plot the profile area. Use top keyword to change the y-axis on both axes
fig, ax, ax_water = profile.plot_area_profile(show_water_ppm=True, centered=False, phase='ol',
calibration='Bell', scale_water=3., ytop=40)
# If you want to use your own area estimates, you can set them directly
profile.areas = [30, 30, 32, 30]
profile.plot_area_profile(axes=ax, centered=False, label='areas set directly',
style=styles.style_points1)
# Now let's add your lovely SIMS data
# Note that I'm plotting on ax_water. If you use ax, the SIMS data plots on the left axis.
SIMS_data = [20, 18, 16, 14]
ax_water.plot(profile.positions_microns, SIMS_data, 'sk', markersize=10, label='SIMS')
ax.legend(loc=3);
In [23]:
fig, ax, ax_water = profile.plot_area_profile(style=styles.style_points1, ignore_idx=[2])
The right and left hand axes are independent. To scale them up simultaneouly, use the ytop keyword when making the plot.
A dictionary of syle terms can be passed in through the keyword style. Several such dictionaries are stored in pynams.styles
In [24]:
fig, ax, ax2 = profile.plot_area_profile(normalize_areas=True, normalize_positions=True,
ytop=1.5, style=styles.style_points, scale_water=7,
calibration='Withers')
In [25]:
avespec = profile.average_spectra()
avespec.get_thickness_from_SiO()
fig, ax = avespec.plot_spectrum()
Perhaps you want to make a baseline and then not have to re-make it every time you want to plot the areas. Perhaps you are using different baselines to come up with rough peak-specific area profiles. Perhaps you want a sense for the range of the calculated area given different baselines. You need a way to handle multiple sets of baselines.
In [26]:
spectrum_augite1.make_baseline(wn_low=3400, force_through_wn=3500)
spectrum_augite1.plot_showbaseline()
spectrum_augite1.save_baseline(folder=FTIR_folder, baseline_ending='-baseline2.CSV')
In [27]:
# make a linear baseline and save it under a different name
spectrum_augite1.make_baseline()
spectrum_augite1.plot_showbaseline()
spectrum_augite1.save_baseline(folder=FTIR_folder, baseline_ending='-linearbaseline.CSV')
In [28]:
spectrum_augite1.get_baseline(folder=FTIR_folder)
fig, ax = spectrum_augite1.plot_showbaseline()
In [29]:
spectrum_augite1.get_baseline(folder=FTIR_folder, baseline_ending='-linearbaseline.CSV')
fig, ax = spectrum_augite1.plot_showbaseline()
In [30]:
profile.save_baselines(folder=FTIR_folder)
In [31]:
profile.get_baselines(folder=FTIR_folder)