scipyconvex-hull

how to apply continuum removal in spectral graph


I have to apply continuum removal on a graph and I have used scipy convexhull function to find convex hull, now i have to apply continuum removal.

here is the code-

import pandas as pd
import numpy as np
from scipy.spatial import ConvexHull
import matplotlib.pyplot as plt

data=open('15C80D4_00002.txt')
d=pd.read_table(data, sep=r'\t',header=None, names=['Wvl', 'Reflectance'],skiprows=1, 
engine='python')

x=d.iloc[:,:1]
a1=np.array(x)

y=d.iloc[:,1:]
b1=np.array(y)

points=np.concatenate((a1,b1), axis=1)


fig = plt.figure()
ax = fig.subplots()

hull = ConvexHull(points)
for simplex in hull.simplices:
    ax.plot(points[simplex,0], points[simplex,1], 'k-')

on plotting the graph i get convex hull graph

  1. i dont want the bottom line, only upper part
  2. i want the graph to be something like this picture, the graph should come in same axis after continuum removal

how can this be done?


Solution

  • Continuum removal equates to division of the spectrum by its convex hull. The basic idea in the implementation below is that we add two points, one at each end of the line, which have a lower y-value than any other point, and are thus guaranteed to form the "lower" part of the convex hull.
    After the convex hull computation, we simply remove them again and are left with the "upper" part of the hull. From there you just have to interpolate the hull at the given x-coordinates to get a corresponding y' value, which you subtract from the original y-values.

    enter image description here

    import numpy as np
    import matplotlib.pyplot as plt
    
    from scipy.spatial import ConvexHull
    from scipy.interpolate import interp1d
    
    def continuum_removal(points, show=False):
        x, y = points.T
        augmented = np.concatenate([points, [(x[0], np.min(y)-1), (x[-1], np.min(y)-1)]], axis=0)
        hull = ConvexHull(augmented)
        continuum_points = points[np.sort([v for v in hull.vertices if v < len(points)])]
        continuum_function = interp1d(*continuum_points.T)
        yprime = y / continuum_function(x)
    
        if show:
            fig, axes = plt.subplots(2, 1, sharex=True)
            axes[0].plot(x, y, label='Data')
            axes[0].plot(*continuum_points.T, label='Continuum')
            axes[0].legend()
            axes[1].plot(x, yprime, label='Data / Continuum')
            axes[1].legend()
    
        return np.c_[x, yprime]
    
    x = np.linspace(0, 1, 100)
    y = np.random.randn(len(x))
    points = np.c_[x, y]
    new_points = continuum_removal(points, show=True)
    
    plt.show()