Peak-specific profiles

New to pynams? Start with the intro and basic examples

Other peak-fitting programs

My personal experience with peakfit.m was that it often gave weird and clearly incorrect results and I ended up hand-fitting all the spectra individually for my cpx paper anyway. With pynams I've tried to make that process slightly less painful.

Some setup


In [1]:
%matplotlib inline
%config InlineBackend.figure_formats=["svg"]
from __future__ import print_function, division
import matplotlib.pyplot as plt
import pynams
from pynams import Spectrum, Profile
import os 

FTIR_file_location = os.path.join(pynams.__path__[0], 'example_FTIR_spectra\\')


spectrum_augite1 = Spectrum(fname='augite1', folder=FTIR_file_location, thickness_microns=876.4)
spectrum_augite2 = Spectrum(fname = 'augite2', folder=FTIR_file_location, thickness_microns=876.4)

profile = Profile(fnames=['olivine1', 'olivine2', 'olivine3', 'olivine4'],  
                  folder=FTIR_file_location, 
                  thicknesses_microns=300., 
                  positions_microns=[10., 20., 30., 40.])

profile.get_baselines()


Got baseline  C:\Users\Elizabeth Ferriss\Documents\Code\Pynams\pynams\example_FTIR_spectra\olivine1-baseline.CSV
Got baseline  C:\Users\Elizabeth Ferriss\Documents\Code\Pynams\pynams\example_FTIR_spectra\olivine2-baseline.CSV
Got baseline  C:\Users\Elizabeth Ferriss\Documents\Code\Pynams\pynams\example_FTIR_spectra\olivine3-baseline.CSV
Got baseline  C:\Users\Elizabeth Ferriss\Documents\Code\Pynams\pynams\example_FTIR_spectra\olivine4-baseline.CSV

Locate peaks


In [2]:
my_peak_locations = spectrum_augite1.find_peaks()
print('peak wavenumber locations')
print(my_peak_locations)


Making default linear baseline
Use Spectrum.make_baseline for other baseline options
peak wavenumber locations
[ 3245.663  3459.726  3546.509  3629.434]

Option 1: Use multiple baselines for well-separated peaks

You can make and save different baselines, one under each peak or set of peaks that you are interested in, and do profile.get_baselines() before each plotting step, as described in the basic pynams examples.

Option 2: Determine peak heights only


In [3]:
spectrum_augite1.make_peakheights(my_peak_locations)

for loc, height in zip(my_peak_locations, spectrum_augite1.peak_heights):
    print('peak at', loc, 'cm-1 is', height, 'cm-1 from baseline')


peak at 3245.663 cm-1 is 0.125212582944 cm-1 from baseline
peak at 3459.726 cm-1 is 2.78405925037 cm-1 from baseline
peak at 3546.509 cm-1 is 2.31683741179 cm-1 from baseline
peak at 3629.434 cm-1 is 2.74750112798 cm-1 from baseline

Option 3: Use Gaussians

The default uses the positions found by find_peaks and widths = 50/cm wavenumbers


In [4]:
fig, ax = spectrum_augite1.make_peakfit()


peak positions: [ 3245.663  3459.726  3546.509  3629.434]
heights: [ 0.12521258  2.78405925  2.31683741  2.74750113]

Make whatever Gaussians you want


In [5]:
fig, ax = spectrum_augite1.make_peakfit(peak_positions=[3616, 3560, 3460],
                                        peak_heights=[2.6, 1.3, 2.8],
                                        peak_widths=[60, 60, 150])


Save your peaks

Keywords are similar to save_baseline, except baseline_ending is peak_ending.


In [6]:
spectrum_augite1.save_peakfit()


Saved C:\Users\Elizabeth Ferriss\Documents\Code\Pynams\pynams\example_FTIR_spectra\augite1-peakfit.CSV

Get saved peaks


In [7]:
spectrum_augite1.get_peakfit()


Got peak info from C:\Users\Elizabeth Ferriss\Documents\Code\Pynams\pynams\example_FTIR_spectra\augite1-peakfit.CSV

Plot your peaks without changing anything


In [8]:
fig, ax = spectrum_augite1.plot_peakfit()


Copy peaks from one spectrum to another


In [9]:
spectrum_augite2.make_peakfit_like(spectrum=spectrum_augite1)
fig, ax = spectrum_augite2.plot_peakfit()
ax.set_ylim(0, 3)


Making default linear baseline
Use Spectrum.make_baseline for other baseline options
Out[9]:
(0, 3)

Peak fitting in profiles

For proper plotting later, all peak positions should be consistent across the profile. For spectra with no peaks, include a placeholder, e.g., by setting the height and/or width to zero while still including the peak position.

Suggestion: Start by fitting peaks to the average spectra across the profile


In [10]:
average_spectrum = profile.average_spectra()
average_spectrum.make_baseline(wn_low=3100, wn_high=3600, force_through_wn=3400)
#average_spectrum.plot_showbaseline()
fig, ax = average_spectrum.make_peakfit(peak_positions=[3230, 3329, 3358, 3436, 3452, 3488, 3525, 3542, 3563, 3573],
                                        peak_heights=[0.045, 0.015, 0.011, 0.01, 0.014, 0.055, 0.21, 0.05, 0.15, 0.33],
                                        peak_widths=[40, 15, 10, 10,  10, 30, 25, 12, 15, 10])
fig.set_size_inches(10, 4)


Fiddle with each spectrum in profile.spectra

Remember to keep the peak positions consistent across the profile.


In [11]:
profile.make_baselines(wn_low=3100, wn_high=3600, force_through_wn=3400)
profile.spectra[0].make_peakfit(peak_positions=[3230, 3329, 3358, 3436, 3452, 3488, 3525, 3542, 3563, 3573],
                                peak_heights=[0.06, 0.02, 0.02, 0.01, 0.014, 0.055, 0.21, 0.05, 0.15, 0.33],
                                peak_widths=[40, 15, 10, 10,  10, 30, 25, 12, 15, 10])
profile.spectra[1].make_peakfit(peak_positions=[3230, 3329, 3358, 3436, 3452, 3488, 3525, 3542, 3563, 3573],
                                peak_heights=[0.06, 0.03, 0.03, 0.02, 0.03, 0.055, 0.21, 0.05, 0.15, 0.33],
                                peak_widths=[40, 15, 10, 10,  10, 30, 25, 12, 15, 10])
profile.spectra[2].make_peakfit(peak_positions=[3230, 3329, 3358, 3436, 3452, 3488, 3525, 3542, 3563, 3573],
                                peak_heights=[0.06, 0.02, 0.02, 0.01, 0.02, 0.055, 0.21, 0.05, 0.15, 0.33],
                                peak_widths=[40, 15, 10, 10,  10, 30, 25, 12, 15, 10])
profile.spectra[3].make_peakfit(peak_positions=[3230, 3329, 3358, 3436, 3452, 3488, 3525, 3542, 3563, 3573],
                                peak_heights=[0.01, 0.0, 0.0, 0.01, 0.01, 0.05, 0.16, 0.03, 0.12, 0.22],
                                peak_widths=[40, 15, 10, 10,  10, 30, 25, 12, 15, 10]);


Save your peakfits

You may need to save multiple peakfit endings to go with multiple baselines.


In [12]:
profile.save_peakfits()


Saved C:\Users\Elizabeth Ferriss\Documents\Code\Pynams\pynams\example_FTIR_spectra\olivine1-peakfit.CSV
Saved C:\Users\Elizabeth Ferriss\Documents\Code\Pynams\pynams\example_FTIR_spectra\olivine2-peakfit.CSV
Saved C:\Users\Elizabeth Ferriss\Documents\Code\Pynams\pynams\example_FTIR_spectra\olivine3-peakfit.CSV
Saved C:\Users\Elizabeth Ferriss\Documents\Code\Pynams\pynams\example_FTIR_spectra\olivine4-peakfit.CSV

Get your saved peakfits


In [13]:
profile.get_peakfits()


Got peak info from C:\Users\Elizabeth Ferriss\Documents\Code\Pynams\pynams\example_FTIR_spectra\olivine1-peakfit.CSV
Got peak info from C:\Users\Elizabeth Ferriss\Documents\Code\Pynams\pynams\example_FTIR_spectra\olivine2-peakfit.CSV
Got peak info from C:\Users\Elizabeth Ferriss\Documents\Code\Pynams\pynams\example_FTIR_spectra\olivine3-peakfit.CSV
Got peak info from C:\Users\Elizabeth Ferriss\Documents\Code\Pynams\pynams\example_FTIR_spectra\olivine4-peakfit.CSV

Plot peaks across a profile by specifying the peak index number


In [14]:
fig, ax_left, ax_right = profile.plot_area_profile(peak_idx=0)


Plot the peak heights instead of the peak areas


In [15]:
fig, ax_right, ax_left = profile.plot_area_profile(peak_idx=0, heights_instead=True)


Group multiple peaks together

In this example, the two [Ti] peaks at wavenumbers 3573 and 2525/cm and are at peak indices 9 and 6. You can group them together to create a composite peak of index 10. The composite peak's 'position' will be the sum of the wavenumbers of the component peaks


In [16]:
profile.make_composite_peak(peak_idx_list=[6, 9])
fig, ax_left, ax_right = profile.plot_area_profile(peak_idx=10)


Peak-specific diffusion modeling

It's just like regular diffusion modeling in pynams, but you specify the peak index.


In [17]:
profile.fitD(log10Dm2s=-14, peak_idx=0, varyD=False, vary_time=True, vary_initial=False);


initial: 2.55+/-0
final: 0.00+/-0
log10D m2/s: -14.00+/-0
minutes: 45.53+/-25.72