Smoothing data for determining peak values in Python

Question

I have a transect with peaks and trough, and want to determine the peak values of both. The dataset has quite some noise so currently, the peaks do not return as a single value. I tried to smooth the data with a rolling mean, and even though the outcome is better than without smoothing, there are still multiple 'peaks'. [CSV file here]

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from scipy.signal import argrelextrema
from pandas import read_csv
from numpy import mean
from matplotlib import pyplot
import csv

df = pd.read_csv('transect2.csv', delimiter=',', header=None, names=['x', 'y'])

plt.plot(df['x'], df['y'], label='Original Height')

rolling = df.rolling(window=100)
rolling_mean = rolling.mean()

plt.xlabel('Distance')
plt.ylabel('Height')

plt.plot(rolling_mean['x'], rolling_mean['y'], label='Mean Height 100')
plt.legend(loc='upper left')
plt.show()

n=1000

ilocs_min = argrelextrema(rolling_mean.y.values, np.less_equal, order=n)[0]
ilocs_max = argrelextrema(rolling_mean.y.values, np.greater_equal, order=n)[0]

df.y.plot (color='gray')
df.iloc[ilocs_max].y.plot(style='.', lw=10, color='red', marker="v");
df.iloc[ilocs_min].y.plot(style='.', lw=10, color='green', marker="^");

Smoothing the data further will not represent the reality so either I can improve this smoothing or use a different smoothing function.

Answer 1

My first instinct is to use Savitzky-Golay filter for smoothing. The second is to forget the argrelextrema when you have a noisy dataset. I have never had any good results using it this way. Better alternative is find_peaks or find_peaks_cwt.

I worked out:

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from scipy.signal import argrelextrema
from scipy.signal import savgol_filter, find_peaks, find_peaks_cwt
from pandas import read_csv
import csv

df = pd.read_csv('transect2.csv', delimiter=',', header=None, names=['x', 'y'])

plt.plot(df['x'], df['y'], label='Original Height')

#apply a Savitzky-Golay filter
smooth = savgol_filter(df.y.values, window_length = 351, polyorder = 5)

#find the maximums
peaks_idx_max, _ = find_peaks(smooth, prominence = 0.01)

#reciprocal, so mins will become max
smooth_rec = 1/smooth

#find the mins now
peaks_idx_mins, _ = find_peaks(smooth_rec, prominence = 0.01)

plt.xlabel('Distance')
plt.ylabel('Height')


plt.plot(df['x'], smooth, label='smoothed')

#plot them
plt.scatter(df.x.values[peaks_idx_max], smooth[peaks_idx_max], s = 55,
            c = 'green', label = 'max')
plt.scatter(df.x.values[peaks_idx_mins], smooth[peaks_idx_mins], s = 55,
            c = 'black', label = 'min')
plt.legend(loc='upper left')
plt.show()

which outputs to this