Scipy linregress回归线与散点图拟合不准 相关系数异常问题

我尝试对数据集ocean_rr_1d和ocean_PBLH_1d（各含330万个数据点）绘制线性回归线，但生成的回归线完全不符合数据分布。两组数据分别是行星边界层高度（PBLH，单位米）和降雨率（单位mm/hr）的日度测量值，PBLH原始为小时级数据，我将逐小时数据拼接为完整日度数组，处理逻辑见下方代码。

我同时测试了两组数据的相关性，使用np.corrcoef(ocean_rr_1d, ocean_PBLH_1d)计算得到的相关系数为0.9998，接近完全相关，但绘图结果显然不符合完全相关的特征，理想的R=1相关数据集分布见下图：

我使用的是掩膜数据，数据处理和绘图代码如下。我已经尝试掩膜np.nan值但问题未解决，也通过print('Daily minimum PBLH:' + ' ' + str(np.min(ocean_PBLH_1d)))检查PBLH最小值为253，和图中最小PBLH点一致，不存在异常值拉偏回归线的情况。

import os
import subprocess, shlex
import conda
import datetime as dt
import numpy as np
import matplotlib.pyplot as plt
import glob
import pandas as pd
import netCDF4
import sys
from global_land_mask import globe
from netCDF4 import Dataset
from scipy import stats

# Fixes KeyError: "PROJ_LIB" issue
conda_file_dir         = conda.__file__
conda_dir              = conda_file_dir.split('lib')[0]
proj_lib               = os.path.join(os.path.join(conda_dir, 'share'), 'proj')
os.environ["PROJ_LIB"] = proj_lib
from mpl_toolkits.basemap import Basemap, cm

curDT = dt.datetime(2020, 7, 20, 0) # Start time 
endDT = dt.datetime(2020, 7, 20, 23) # End time

directory        = '/uufs/chpc.utah.edu/common/home/zpu-group10/pu03/mpletch/thesis/IMERG/July/gpm1.gesdisc.eosdis.nasa.gov/opendap/GPM_L3/GPM_3IMERGHH.06/2020/202/*.nc4' # Directory for IMERG data
regridded_merra2 = Dataset('/uufs/chpc.utah.edu/common/home/zpu-group10/pu03/mpletch/thesis/regridded_MERRA2_pblh/07' + curDT.strftime("%d") + '_new_regridded_MERRA2.nc4.nc4')
hourly           = ("00.V06B.HDF5.nc4", "20.V06B.HDF5.nc4", "40.V06B.HDF5.nc4", "60.V06B.HDF5.nc4", "80.V06B.HDF5.nc4") # Strings that represent tshe pattern of hourly rain-rate data since 30-minute data is included in the IMERG data.

rain_rate_all      = []
regridded_pblh_all = []
lat_regridded_all  = []
lon_regridded_all  = []

while curDT < endDT: # Loops through each hourly IMERG case by adding an hour to the current time (curDT). 
  for filepath in glob.iglob(directory): # Navigates to the IMERG data directory
    if filepath.endswith(hourly): # Ensures that only hourly and not half-hourly data is used
      
      imerg_data    = Dataset(filepath, 'r')
      imerg_data.variables.keys()
      rain_rate     = np.array(imerg_data.variables['precipitationCal'][0,100:650,0:250]*1).T # Reads in IMERG rain rate data
      lat_imerg     = np.array(imerg_data.variables['lat'][0:250]*1)
      lon_imerg     = np.array(imerg_data.variables['lon'][100:650]*1)
      levels_imerg  = [2 + 2*n for n in range(25)] # Defines contours for rain rate plots
 
      regridded_pblh = np.array(regridded_merra2.variables['PBLH'][curDT.strftime("%-H"),0:250,100:650])
      regridded_lat  = regridded_merra2.variables['lat'][0:250]
      regridded_lon  = regridded_merra2.variables['lon'][100:650]
      regridded_lat  = np.tile(regridded_lat,(len(regridded_lon),1)).T
      regridded_lon  = np.tile(regridded_lon,(len(regridded_lat),1))

      rain_rate_all.append(rain_rate)
      pblh_all.append(pblh)
      regridded_pblh_all.append(regridded_pblh)
      lat_regridded_all.append(regridded_lat)
      lon_regridded_all.append(regridded_lon)
   
    curDT = curDT + dt.timedelta(hours = 0.5) # Loops through hourly PBLH data instead of every two hours

# Converts all arrays to np.arrays.
regridded_pblh_all = np.array(regridded_pblh_all)
rain_rate_all      = np.array(rain_rate_all)
lat_regridded_all  = np.array(lat_regridded_all)
lon_regridded_all  = np.array(lon_regridded_all)

masked_rr        = np.where(rain_rate_all >=10, rain_rate_all, -1) # Masks any rain rate that is less than or equal to 10 mm/hr and replaces them with -1
masked_rr        = np.ma.masked_where(masked_rr == -1, masked_rr) # Masks -1 PBLH values.
onedim_masked_rr = masked_rr.flatten()
PBLH_col         = np.where(rain_rate_all > 10, regridded_pblh_all, -1) # Locates the PBLH where the rainfall rate is 10 mm/hr and greater
PBLH_col         = np.ma.masked_where(PBLH_col == -1, PBLH_col) # Masks the -1 fill values
PBLH_col         = np.ma.masked_where(PBLH_col == 1000000000000000.0, PBLH_col) # Masks the 1E15 values

on_land = globe.is_land(lat_regridded_all, lon_regridded_all) # Determines if a point is on land

ocean_PBLH    = np.where(on_land == False, PBLH_col, -9999) # Finds the PBLH over ocean instead of land.
ocean_PBLH    = np.ma.masked_where(ocean_PBLH < 0, ocean_PBLH) # Masks any PBLH below 0
ocean_PBLH    = np.ma.masked_where(ocean_PBLH == 1000000000000000.0, ocean_PBLH) # Masks 1E15 PBLH values that are for some reason in the data
ocean_PBLH_1d = ocean_PBLH.flatten() # Creates a 1-D array to use for correlation test

ocean_rr      = np.where(on_land == False, masked_rr, -9999) # Masks all data over land.
ocean_rr      = np.ma.masked_where(ocean_rr < 0, ocean_rr) # Masks out any negative values.
ocean_rr_1d   = ocean_rr.flatten() # Creates a 1-D array to use for the correlation test.

corr = np.corrcoef(ocean_rr_1d, ocean_PBLH_1d) # Runs correlation between the two data sets.

regression_mask = ~np.isnan(ocean_rr_1d) & ~np.isnan(ocean_PBLH_1d)
slope, intercept, r_value, p_value, std_err = stats.linregress(ocean_rr_1d[regression_mask], ocean_PBLH_1d[regression_mask]) # Masks any NaN values to be used in the correlation
plt.scatter(ocean_rr_1d, ocean_PBLH_1d, s = 2, c = 'blue')
plt.plot(ocean_rr_1d, intercept + slope*ocean_rr_1d, 'r')
plt.title('Scatter plot of rain rate (10+ mm/hr) and PBLH (m)\n 07-04-2020')
plt.xlabel('Rain rate (mm/hr)')
plt.ylabel('PBLH (m)')
plt.show()

以下是两个数据集的样本输出，由于掩膜值较多，我选取了一小段非掩膜数据展示格式：

PBLH数据 #

3235075          NaN
3235076          NaN
3235077   605.557373
3235078   545.571350
3235079   485.585358
3235080   425.599365
3235081          NaN
3235082          NaN

降雨率数据 #

3235075        NaN
3235076        NaN
3235077  26.847431
3235078  27.215666
3235079  14.650899
3235080  32.269173
3235081        NaN
3235082        NaN

内容的提问来源于stack exchange，提问作者mpletch1