I am attempting to plot a linear regression line in my data (ocean_rr_1d and ocean_PBLH_1d, 3,300,000 points in each data set), but the line does not match the pattern of the data at all. Below is a plot of the data. The data are daily measurements of planetary boundary layer height (PBLH, units of meters) and rain rate (mm/hr). Since the PBLH data are hourly data, I had to append the hourly data into an array to create a full day's worth of data. See code below for how I did this.
I am also testing the correlation between these two datasets. When I use np.corrcoef(ocean_rr_1d, ocean_PBLH_1d), it outputs 0.9998, which is a nearly perfect correlation between the two variables. However, the data are clearly not plotting in a fashion that represents a perfect correlation.
One would expect a perfectly correlated data set to look similar to the plot below:
The data that I am using is masked data, and below is the code for it as well as the code I used to create the plot. I have attempted to mask the np.nan values to see if that fixes the issue but it does not do anything. I have also checked to see what the minimum value of PBLH is by using print('Daily minimum PBLH:' + ' ' + str(np.min(ocean_PBLH_1d))) to check to see if the line is incorporating any outliers. However, it outputs 253 which is same as the smallest PBLH point on the plot.
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()
Here is a sample output of each data set. Since there are so many masked values, I took a small sample from each dataset to show the non-masked data format.
PBLH data:
3235075 NaN
3235076 NaN
3235077 605.557373
3235078 545.571350
3235079 485.585358
3235080 425.599365
3235081 NaN
3235082 NaN
Rain rate data:
3235075 NaN
3235076 NaN
3235077 26.847431
3235078 27.215666
3235079 14.650899
3235080 32.269173
3235081 NaN
3235082 NaN
Please let me know if you need any other information from me. I greatly appreciate any help.