如何用估计的背景值替换二维数组的异常值?

How to replace outliers of 2d array with estimated background values?

我有一个名为 no2 的二维数组,它与其他两个二维数组 szavza.

相关

测试数据(test.npz, 450 KB)可以从Google Drive.

下载

概述如下:

import numpy as np
import matplotlib.pyplot as plt

data = np.load('test.npz')
sza = data['sza']
vza = data['vza']
no2 = data['no2']

fig, axs = plt.subplots(2, 2, figsize=(8, 6))

ax1, ax2, ax3, ax4 = axs.flat

m = ax1.pcolormesh(no2)
plt.colorbar(m, ax=ax1)
ax1.set_title('no2')

m = ax2.pcolormesh(sza)
plt.colorbar(m, ax=ax2)
ax2.set_title('sza')

m = ax3.pcolormesh(vza)
plt.colorbar(m, ax=ax3)
ax3.set_title('vza')

s = ax4.scatter(sza, no2, c=vza, s=1)
plt.colorbar(s, ax=ax4, label='vza')
ax4.set_xlabel('sza')
ax4.set_ylabel('no2')

plt.tight_layout()

我想根据周围的背景或低 no2 值替换两个高 no2 区域以获得如下内容:

因为看起来 no2 线性依赖于 sza,如上一个子图所示,我想出了三个想法:

曲线拟合

使用 no2sza 之间的拟合与几个 vza 箱来计算背景 no2 以替换高 no2 值:

fig, axs = plt.subplots(3, 4, figsize=(12, 6))
ax = axs.flat

for index,bin in enumerate(range(5, 65, 5)):
    mask = (vza>bin)&(vza<bin+5)
    # print(index)
    s = ax[index].scatter(sza[mask], no2[mask], c=vza[mask], s=1)
    plt.colorbar(s, ax=ax[index], label='vza')
    ax[index].set_title(str(bin)+'<vza<'+str(bin+5))

for ax in axs.flat:
    ax.set_xlabel('sza')
    ax.set_ylabel('no2')

plt.tight_layout()

我尝试拟合一个 bin 的曲线 (45

from scipy.optimize import curve_fit

def func(x, a, b, c):
    return a * np.exp(-b * x) + c

xdata = sza[(vza>45)&(vza<50)]
ydata = no2[(vza>45)&(vza<50)]
popt, pcov = curve_fit(func, xdata, ydata, p0=(1, 1e-5, 1))

plt.plot(xdata, ydata, 'b-', label='data')

plt.plot(xdata, func(xdata, *popt), 'r-',
         label='fit: a=%5.3f, b=%5.3f, c=%5.3f' % tuple(popt))
plt.legend()

然而,它没有得到我想要的:

是否可以满足以下两个条件?

- 拟合曲线并获取高值的背景值

- 将随机噪声添加到拟合的背景值(这可以 运行 多次以获得更真实的值,如周围的背景值)

或者其他更好的方法?

渐变

我检查了梯度,希望它能使高值更显着:

# 
grad = np.gradient(no2)
fulgrad = np.sqrt(grad[0]**2 + grad[1]**2)

fig, axs = plt.subplots(1, 2, figsize=(6, 3))

ax1, ax2  = axs.flat

m = ax1.pcolormesh(no2)
plt.colorbar(m, ax=ax1)
ax1.set_title('no2')

m = ax2.pcolormesh(fulgrad)
plt.colorbar(m, ax=ax2)
ax2.set_title('no2 gradient')

plt.tight_layout()

但是,它只能显示一些轮廓:

图像处理

我不知道如何使用 scikit-learn 只替换高值并保持背景不变。

最后,我想出了如何用估计的背景值替换高值。

只需使用 scikit-ued 中的 dual-tree complex wavelet transform

import numpy as np
import matplotlib.pyplot as plt
from skued import baseline_dt

data = np.load('../data/test.npz')

baseline = baseline_dt(data['no2'], wavelet = 'qshift3', level = 6, max_iter = 150)

fig, axs = plt.subplots(1, 3, figsize=(12, 4))

ax1, ax2, ax3 = axs.flat

m = ax1.imshow(data['no2'], vmin=0, vmax=7e-4)
plt.colorbar(m, ax=ax1)
ax1.set_title('no2')

m = ax2.imshow(baseline, vmin=0, vmax=7e-4)
plt.colorbar(m, ax=ax2)
ax2.set_title('baseline')

m = ax3.imshow(data['no2']-baseline, vmin=0, vmax=7e-4)
plt.colorbar(m, ax=ax3)
ax3.set_title('no2 - baseline')