scipy 如何在谱图中应用连续谱消除

b5lpy0ml  于 2022-11-10  发布在  其他
关注(0)|答案(1)|浏览(119)

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?
qaxu7uf2

qaxu7uf21#

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()

相关问题