如何为NumPy二维散点图拟合直线并筛选离群点？

问题重现 #

用户尝试对二维网格上的两类观测数据（TIR1和TIR2亮温）拟合直线，使用numpy.polyfit时触发类型错误，核心报错为：

TypeError: expected 1D vector for x

用户代码如下：

import h5py
import sys
import numpy as np
import matplotlib.pyplot as plt

# First get data from HDF5 file with h5py:
fn = '/home/swadhin/project/insat/data/3RIMG_30MAR2018_0014_L1B_STD_V01R00.h5'
x = np.load('/home/swadhin/project/insat/land_mask.npy')

#Reading the L2 SST data
with np.load('sst.npz') as npz:
    arr_sst = np.ma.MaskedArray(**npz)

with h5py.File(fn) as f:
    print(list(f.keys()))
    # retrieve image data:
    image = 'IMG_TIR1'
    image2 = 'IMG_TIR2'
    lon = f['Longitude'][:]*0.01
    lat = f['Latitude'][:]*0.01
    img_arr = f[image][0, :, :]
    img_arr2 = f[image2][0, :, :]
    
    # get _FillValue for data masking
    img_arr_fill = f[image].attrs['_FillValue'][0]
    img_arr_fill2 = f[image2].attrs['_FillValue'][0]

    #LUT 
    radiance_lut_tir1 = np.array(f[image+str('_RADIANCE')])
    radiance_lut_tir2 = np.array(f[image2+str('_RADIANCE')])
    bt_lut_tir1 = np.array(f[image+str('_TEMP')])
    bt_lut_tir2 = np.array(f[image2+str('_TEMP')])
    
# retrieve the extent of the plot from file attributes:
    left_lon = f.attrs['left_longitude'][0]
    right_lon = f.attrs['right_longitude'][0]
    lower_lat = f.attrs['lower_latitude'][0]
    upper_lat = f.attrs['upper_latitude'][0]
    sat_long = f.attrs['Nominal_Central_Point_Coordinates(degrees)_Latitude_Longitude'][1]
    sat_hght = f.attrs['Nominal_Altitude(km)'][0] * 1000.0  # (for meters)
print('Done reading HDF5 file')

# Use np.ma.masked_equal with integer values to
# mask '_FillValue' data in corners:
img_arr_m = np.ma.masked_equal(img_arr, img_arr_fill)
img_arr_m2 = np.ma.masked_equal(img_arr2, img_arr_fill2)

lon_m = np.ma.masked_equal(lon, 327.67)
lat_m = np.ma.masked_equal(lat, 327.67)
#########BT from LUT

def count2bt_lut_tir1(count):
    return bt_lut_tir1[count]
def count2bt_lut_tir2(count):
    return bt_lut_tir2[count]


bt_tir1_lut = count2bt_lut_tir1(img_arr_m)
bt_tir2_lut = count2bt_lut_tir2(img_arr_m2)

bt_tir1_lut_m = np.ma.array(bt_tir1_lut, mask=x, dtype=np.float64)
bt_tir2_lut_m = np.ma.array(bt_tir2_lut, mask=x, dtype=np.float64)

################# Creating a Scatter Plot of Bts from TIR1 and TIR2

plt.scatter(bt_tir1_lut_m,bt_tir2_lut_m,s = 10)
np.shape(bt_tir1_lut_m)
m, b = np.polyfit(bt_tir1_lut_m, bt_tir2_lut_m, 1)
plt.plot(bt_tir1_lut_m, m*bt_lut_tir1 + b,color='k')


plt.xlabel("Brightness Temperature(K) in TIR1 (10.8 micron)")
plt.ylabel("Brightness Temperature(K) in TIR2 (12 micron)")

完整报错栈：

TypeError                                 Traceback (most recent call last)
/home/swadhin/project/insat/bt_comparison.py in <cell line: 156>()
    154 plt.scatter(bt_tir1_lut_m,bt_tir2_lut_m,s = 10)
    155 np.shape(bt_tir1_lut_m)
--> 156 m, b = np.polyfit(bt_tir1_lut_m, bt_tir2_lut_m, 1)
    157 plt.plot(bt_tir1_lut_m, m*bt_lut_tir1 + b,color='k')
    160 plt.xlabel("Brightness Temperature(K) in TIR1 (10.8 micron)")

File <__array_function__ internals>:180, in polyfit(*args, **kwargs)

File ~/anaconda3/envs/rttov/lib/python3.9/site-packages/numpy/lib/polynomial.py:636, in polyfit(x, y, deg, rcond, full, w, cov)
    634     raise ValueError("expected deg >= 0")
    635 if x.ndim != 1:
--> 636     raise TypeError("expected 1D vector for x")
    637 if x.size == 0:
    638     raise TypeError("expected non-empty vector for x")

TypeError: expected 1D vector for x

解决方案 #

1. 拟合直线：将二维MaskedArray转为一维有效数据 #

np.polyfit要求输入为一维数组，而bt_tir1_lut_m和bt_tir2_lut_m是二维MaskedArray。需要先提取其中非掩码的有效数据并展平为一维：

# 提取非掩码的一维有效数据
x_data = bt_tir1_lut_m.compressed()
y_data = bt_tir2_lut_m.compressed()

# 拟合直线（一阶多项式）
m, b = np.polyfit(x_data, y_data, 1)

修正拟合线绘制

原代码中plt.plot(bt_tir1_lut_m, m*bt_lut_tir1 + b)存在错误：bt_lut_tir1是LUT数组，不是当前的亮温数据。正确做法是基于x数据的范围生成拟合线的点：

# 生成拟合线的x轴范围（覆盖所有有效数据）
x_fit = np.linspace(x_data.min(), x_data.max(), 100)
y_fit = m * x_fit + b

# 绘制拟合线
plt.plot(x_fit, y_fit, color='k', label=f'Fit: y = {m:.2f}x + {b:.2f}')
plt.legend()

2. 获取远离拟合线的点（异常点） #

通过计算残差（实际值与拟合值的差），设定阈值筛选异常点：

# 计算所有有效点的拟合值与残差
y_pred = m * x_data + b
residuals = y_data - y_pred

# 设定异常点阈值（示例：残差绝对值超过2倍标准差）
threshold = 2 * np.std(residuals)
outlier_mask = np.abs(residuals) > threshold

# 提取异常点
x_outliers = x_data[outlier_mask]
y_outliers = y_data[outlier_mask]

# 在散点图中标记异常点
plt.scatter(x_outliers, y_outliers, s=15, color='r', label='Outliers')
plt.legend()

完整修正后的代码片段 #

################# Creating a Scatter Plot of Bts from TIR1 and TIR2

plt.scatter(bt_tir1_lut_m, bt_tir2_lut_m, s=10, alpha=0.5)

# 提取有效数据并拟合直线
x_data = bt_tir1_lut_m.compressed()
y_data = bt_tir2_lut_m.compressed()
m, b = np.polyfit(x_data, y_data, 1)

# 绘制拟合线
x_fit = np.linspace(x_data.min(), x_data.max(), 100)
y_fit = m * x_fit + b
plt.plot(x_fit, y_fit, color='k', label=f'Fit: y = {m:.2f}x + {b:.2f}')

# 检测并绘制异常点
y_pred = m * x_data + b
residuals = y_data - y_pred
threshold = 2 * np.std(residuals)
outlier_mask = np.abs(residuals) > threshold
plt.scatter(x_data[outlier_mask], y_data[outlier_mask], s=15, color='r', label='Outliers')

plt.xlabel("Brightness Temperature(K) in TIR1 (10.8 micron)")
plt.ylabel("Brightness Temperature(K) in TIR2 (12 micron)")
plt.legend()
plt.show()

内容的提问来源于stack exchange，提问作者The Emerging Star