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
how can this be done?
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.
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()