为 numpy 2D 散点图绘制最佳拟合线

Question

我试图在散点图中拟合一条线，其中 X 和 Y 坐标是相同维度的 2D NumPy 数组。这里的 X 和 Y 是在同一组网格点上记录的两种不同类型的观测值。我尝试使用 numpy.polyfit 但出现以下错误：

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

如何在这些数据中拟合一条线？那么如何获得远离拟合线的点呢？

Answer 1

如错误中所示，

polyfit

需要一维数组，但

bt_tir1_lut_m

和

bt_tir2_lut_m

包含二维数组。您可以绘制它们，因为在内部，

plt.scatter

可能将它们展平为一维数组，但

polyfit

无法单独做到这一点。

因此，您需要使用

np.reshape(bt_tir1_lut_m, -1)

将数组展平为 1d。

这是一个看起来像这样的示例

import matplotlib.pyplot as plt
import numpy as np

# Generating random 2d arrays
a1 = np.random.rand(10,10)
a2 = np.random.rand(10,10)

# Generating a random mask
mask = np.random.randint(0, 2, size=a1.shape)

# Generating masked arrays
a1_m = np.ma.array(a1, mask=mask, dtype=np.float64)
a2_m = np.ma.array(a2, mask=mask, dtype=np.float64)

# Plotting the 2d arrays 
plt.scatter(a1_m, a2_m)
# Getting the coefficient of the interpolation polynomial, by reshaping
# the masked arrays
m, b = np.polyfit(np.reshape(a1_m, -1), np.reshape(a2_m, -1), 1)

# Plotting the resulting affine function
plt.plot(np.reshape(a1_m, -1), m*np.reshape(a1_m, -1) + b,color='k')
plt.show()