<jupyter_start><jupyter_text>Statistical bias correction for satellite precipitation in British Columbia - an application of Gaussian process regressionYingkai ShaEarth, Ocean and Atmospheric Sciences, University of British Columbia **General steps:*** Process gridded GPM data in Nov 21st 2017 (`3B-HHR.MS.MRG.3IMERG`)* Process BC Hydro in-situ observations (daily precip)* Calculate the additive bias between in-situ obs and GPM (the closest GPM grid at station location).* Split station data (bias and obs) into train/test groups* Configure GPR (lat/lon --> bias) and learn hyper-parameters through training set (at station locations)* Apply GPR to produce bias for all GPM grids, and use gridded bias for GPM correction* Compare the corrected GPM with testing set, report metrics<jupyter_code># import modules
import GPy
import time
import numpy as np
from scipy import stats
from scipy.spatial import cKDTree
from sklearn.preprocessing import StandardScaler
# import plotting libraries
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
# KDTree function
def KDTree_wraper(dist_lon, dist_lat):
    '''
    A warper of scipy.spatial.cKDTree
    Tree = KDTree_wraper(dist_lon, dist_lat)
    '''
    return cKDTree(list(zip(dist_lon.ravel(), dist_lat.ravel())))

# Cartopy map settings
scale_param = '10m' # 10m for publication quality
# BC and Ca-US boundary
PROVINCE = cfeature.NaturalEarthFeature(
    category='cultural',
    name='admin_1_states_provinces_lines',
    scale=scale_param,
    facecolor='none')

% matplotlib inline
B = [0.5, 0.5, 1] # color<jupyter_output><empty_output><jupyter_text>Data<jupyter_code># Preprocessed GPM data
temp = np.load('/Users/yingkai/Repository/WRFT/_data/Satellite/GPM/GPM_BC_Nov_2017.npy')
lon = temp[()]['lon']
lat = temp[()]['lat']
GPM = temp[()]['GPM'][20, ...]
ERROR = temp[()]['ERROR'][20, ...]
TOPO = temp[()]['TOPO'] # TOPO for plotting
ocean_flag = TOPO<0 # flag out ocean
# BC Hydro station observation
temp = np.load('/Users/yingkai/Repository/WRFT/_data/Satellite/GPM/BCH_obs_Nov_2017.npy')
stn_lon = temp[()]['stn_lon']
stn_lat = temp[()]['stn_lat']
stn_obs = temp[()]['stn_obs']
stn_code = np.array(temp[()]['stn_code'])
# train/test split: 20 data poiunts for test
# ===== #
# For reproducing the result, here the "random" is hard coded. Cannon use seed because of the Monte Carlo later. 
## withhold = sorted(np.random.randint(len(stn_lon), size=20)) # <--- original version
stn_test_temp = ['ALU', 'BRI', 'COX', 'DAI', 'DOW', 'ERC', 'FDL', 'FST', 'HEB', 'HFF', \
                 'MSK', 'PNK', 'PYN', 'STV', 'WOL', 'YPW', 'YXJ', 'YZY']
withhold = []
for test_code in stn_test_temp:
    withhold.append(np.where(stn_code==test_code))
# ==== #
withhold_flag = np.in1d(np.arange(len(stn_lon)).astype(int), withhold)
# test
test_lon = stn_lon[withhold_flag]
test_lat = stn_lat[withhold_flag]
test_obs = stn_obs[withhold_flag]
test_code = np.array(stn_code)[withhold_flag]
# train
train_lon = stn_lon[~withhold_flag]
train_lat = stn_lat[~withhold_flag]
train_obs = stn_obs[~withhold_flag]
train_code = np.array(stn_code)[~withhold_flag]
# identify the closest GPM grid for each station
gridTree = KDTree_wraper(lon, lat)
grid_shape = lat.shape
_, inds = gridTree.query(list(zip(stn_lon, stn_lat)))
stn_indx, stn_indy = np.unravel_index(inds, grid_shape)
test_indx = stn_indx[withhold_flag]
test_indy = stn_indy[withhold_flag]
train_indx = stn_indx[~withhold_flag]
train_indy = stn_indy[~withhold_flag]
train_GPM = GPM[train_indx, train_indy]
#train_elev = TOPO[train_indx, train_indy] # elevation as an optional input (not used as default)
# ERROR = the "uncertainty" in GPM retrieval 
train_ERROR = ERROR[train_indx, train_indy]
test_ERROR = ERROR[test_indx, test_indy]
# an example of errors at training stations (see the strong heteroscedasticity)
train_ERROR
# Station plot is not showing here<jupyter_output><empty_output><jupyter_text>GPR model<jupyter_code># lat/lon as input features
GRID = np.vstack([train_lon, train_lat]).T # shape=(stn_num, 3)
# in-situ obs minus GPM (additive bias) as target
Y = (train_obs-train_GPM)[:, None]<jupyter_output><empty_output><jupyter_text>**Configuration**<jupyter_code># GPR configuration:
## zero-mean
## isotropic RBF + anisotropic Matern32 + heteroscedasticity defined by GPM error
## GPM error fixed

kernel = GPy.kern.RBF(input_dim=1, ARD=True)*GPy.kern.Matern32(input_dim=2, ARD=True) + \
         GPy.kern.White(input_dim=1)
# set priors
model = GPy.models.GPHeteroscedasticRegression(GRID, Y, kernel)
model['.*het_Gauss.variance'] = train_ERROR[:, None]**2; model.het_Gauss.variance.fix() # <--- fixed
model['sum.mul.rbf.lengthscale'] = 400 # prior can change, leads to different results
#model['sum.mul.Mat32.lengthscale'] = #
# # ===== #
# # MLE optimization for hyper-parameters. Loop over untill converge
# # Replaced by HMC
# status = 'test'
# while status != 'Converged':
#     out = model.optimize(messages=True)
#     status = out.status
# # ===== #<jupyter_output><empty_output><jupyter_text>**HMC hyper-param learning**<jupyter_code># Hybrid Monte Carlo, momentum of random walk set as 0.025 (adjusted by hand), Markov Chain sample num set as 50000 
start = time.time() # roughly 25 min
hmc = GPy.inference.mcmc.HMC(model, stepsize=2.5e-2) # stepsize sensitive 
s = hmc.sample(num_samples=20000)
end = time.time()
print('HMC wall time: {} sec'.format(end - start))
# 15000-20000 is confirmed to be the range that all hyper-params are stable
plt.plot(s)
plt.axvline(x=15000, ls='--', lw=2.5, c='k')
#plt.axvline(x=20000, ls='--', lw=2.5, c='k')
# update paramters into the GPR model

model.sum.mul.rbf.variance = s[-15000:, 0].mean()
model.sum.mul.rbf.lengthscale = s[-15000:, 1].mean()
model.sum.mul.Mat32.variance = s[-15000:, 2].mean()
model.sum.mul.Mat32.lengthscale = s[-15000:, 3:5].mean(axis=0)
model.sum.white.variance = s[-15000:, 5].mean()
#model.Gaussian_noise.variance = s[-10000:, 6].mean()
# export std
s[-15000:, :].std(axis=0)
model
model.sum.mul.Mat32.lengthscale
model.sum.mul.rbf.lengthscale<jupyter_output><empty_output><jupyter_text>Result<jupyter_code>mean = np.empty(lon.shape)
var = np.empty(lon.shape)
for i in range(lon.shape[0]):
    for j in range(lon.shape[1]):
        mean[i, j], var[i, j] = model._raw_predict(np.array([lon[i, j], lat[i, j]])[None, :])
mean_land = np.copy(mean)
mean_land[ocean_flag] = 0
clevel = np.array([0, 1, 2, 3.5, 5, 7.5, 10, 12.5, 15, 20, 25, 30, 40, 60])


fig = plt.figure(figsize=(16*0.75, 12*0.75))
ax1 = plt.subplot(2, 1, 1, projection=ccrs.PlateCarree())
ax2 = plt.subplot(2, 1, 2, projection=ccrs.PlateCarree())
#ax3 = plt.subplot(3, 1, 3, projection=ccrs.PlateCarree())
AX = [ax1, ax2]
for ax in AX:
    ax.set_extent([-140, -110, 48, 59], ccrs.Geodetic())
    ax.add_feature(cfeature.LAND.with_scale(scale_param), facecolor='w')
    ax.add_feature(cfeature.COASTLINE.with_scale(scale_param), zorder=5)
    ax.add_feature(cfeature.BORDERS.with_scale(scale_param), linestyle='--', linewidth=2.5, zorder=5)
    ax.add_feature(PROVINCE, edgecolor='k', linestyle='-', linewidth=1.5, zorder=5)
    ax.outline_patch.set_linewidth(2.5)
    
    GL = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True, 
                 linewidth=0.5, color='k', linestyle=':')
    GL.xlabels_top = False
    GL.ylabels_right = False
    GL.xlabels_bottom = False
    GL.xformatter = LONGITUDE_FORMATTER
    GL.yformatter = LATITUDE_FORMATTER
    GL.xlabel_style = {'size': 14}
    GL.ylabel_style = {'size': 14}
    GL.ylocator = mticker.FixedLocator([50, 52, 54, 56, 58])
GL.xlabels_bottom = True

CS1 = ax1.contourf(lon, lat, GPM, clevel, cmap=plt.cm.gist_ncar_r, extend='both')
CS2 = ax2.contourf(lon, lat, GPM+mean, clevel, cmap=plt.cm.gist_ncar_r, extend='both')
#CS3 = ax3.contourf(lon, lat, GPM+mean, clevel, cmap=plt.cm.gist_ncar_r, extend='both')
#fig.colorbar(CS, ax=ax1)
ax1.scatter(stn_lon, stn_lat, c=stn_obs, edgecolor='k', linewidth=1.5, s=75, vmin=clevel[0], vmax=clevel[-1], \
           cmap=plt.cm.gist_ncar_r)
ax2.scatter(stn_lon, stn_lat, c=stn_obs, edgecolor='k', linewidth=1.5, s=75, vmin=clevel[0], vmax=clevel[-1], \
           cmap=plt.cm.gist_ncar_r)
#ax3.scatter(stn_lon, stn_lat, c=stn_obs, edgecolor='k', linewidth=1.5, s=75, vmin=clevel[0], vmax=clevel[-1], \
#           cmap=plt.cm.gist_ncar_r)

ax1.set_title('(a) GPM retrieval (original data)', fontsize=14)
ax2.set_title('(b) GPR Prediction', fontsize=14)
#ax3.set_title('(c) GPR Prediction', fontsize=14)

cax = fig.add_axes([0.9, 0.205, 0.015, 0.6]) 
CBar = plt.colorbar(CS1, cax=cax)  
CBar.ax.tick_params(axis='y', labelsize=14, direction='in', length=14)
CBar.set_label('Daily precipitation [mm]', fontsize=14)
plt.tight_layout()

fig.savefig('/Users/yingkai/Desktop/GPR_high_quality/fig2.png', dpi=250, orientation='portrait', \
           papertype='a4', format='png', bbox_inches='tight', pad_inches=0.1)
break_id = [21, 42, 63, 85] 

def break_data(data, break_id):
    return data[:break_id[0]], data[break_id[0]:break_id[1]], \
    data[break_id[1]:break_id[2]], data[break_id[2]:break_id[3]]

fake_x = np.arange(len(train_lon))
fake_x_b = break_data(fake_x, break_id)
train_code_b = break_data(train_code, break_id)
train_GPM_b = break_data(train_GPM, break_id)
train_ERROR_b = break_data(train_ERROR, break_id)
train_obs_b = break_data(train_obs, break_id)

fig = plt.figure(figsize=(16*0.75, 11*0.75))
ax1 = plt.subplot(4, 1, 1)
ax2 = plt.subplot(4, 1, 2)
ax3 = plt.subplot(4, 1, 3)
ax4 = plt.subplot(4, 1, 4)
AX = [ax1, ax2, ax3, ax4]

for i, ax in enumerate(AX):
    ax.grid(linestyle=':'); ax.xaxis.grid(False)
    ax.spines["top"].set_visible(False)
    #ax.spines["bottom"].set_visible(False)
    ax.spines["right"].set_visible(False)
    ax.xaxis.set_tick_params(labelsize=10)
    ax.yaxis.set_tick_params(labelsize=14)
    #ax.set_xlabel('False positive rate', fontsize=14)
    #ax.set_ylabel('Rate of ', fontsize=14)
    [j.set_linewidth(2.5) for j in ax.spines.values()]
    ax.tick_params(axis="both", which="both", bottom=False, top=False, \
               labelbottom=True, left=False, right=False, labelleft=True)

    ax.set_xticks(fake_x_b[i])
    ax.set_xticklabels(train_code_b[i])
    ax.bar(fake_x_b[i], train_obs_b[i], color='w', edgecolor='k', linewidth=1.5, hatch='/', label='Station observation')
    ax.errorbar(fake_x_b[i], train_GPM_b[i], yerr=train_ERROR_b[i], fmt='o', color='r', \
                ms=12, mec='k', capsize=5, elinewidth=3, markeredgewidth=1.5, label='GPM retrieval')
    ax.set_ylim([-10, 90])
    ax.set_yticks([0, 40, 80])
LG = ax4.legend(bbox_to_anchor=(0.75, -0.2), ncol=2, prop={'size':14});
LG.get_frame().set_facecolor('white')
LG.get_frame().set_edgecolor('k')
LG.get_frame().set_linewidth(1.5)
#ax1.set_title('', fontsize=14, y=0.85)
ax3.text(-0.065, 1.05, 'Daily precipitation [mm]', va='center', ha='center', rotation='90', \
         fontsize=14, transform=ax3.transAxes)

fig.savefig('/Users/yingkai/Desktop/GPR_high_quality/fig1.png', dpi=250, orientation='portrait', \
           papertype='a4', format='png', bbox_inches='tight', pad_inches=0.1)
from sklearn.metrics import mean_squared_error
from sklearn.metrics import mean_absolute_error
def cal_verif_stats(stn_assimi, stn_obs):
    # verification stats
    RSME_assim = np.sqrt(mean_squared_error(stn_assimi, stn_obs))
    MAE_assim  = mean_absolute_error(stn_assimi, stn_obs)
    return MAE_assim, RSME_assim
test_mean = (GPM+mean)[test_indx, test_indy]
#test_mean_kg = (GPM+delta_kg)[test_indx, test_indy]
test_ori = (GPM)[test_indx, test_indy]
test_var = np.sqrt((var)[test_indx, test_indy])*1.96
fake_x2 = np.arange(len(test_mean))

fig = plt.figure(figsize=(16*0.75, 4))
ax = fig.gca()

ax.grid(linestyle=':'); ax.xaxis.grid(False)
ax.spines["top"].set_visible(False)
#ax.spines["bottom"].set_visible(False)
ax.spines["right"].set_visible(False)
ax.xaxis.set_tick_params(labelsize=12)
ax.yaxis.set_tick_params(labelsize=14)
ax.set_ylabel('Daily precipitation [mm]', fontsize=14)
[j.set_linewidth(2.5) for j in ax.spines.values()]
ax.tick_params(axis="both", which="both", bottom=False, top=False, \
           labelbottom=True, left=False, right=False, labelleft=True)
ax.set_ylim([-25, 100])
ax.bar(fake_x2, test_obs, color='w', edgecolor='k', linewidth=1.5, hatch='/', label='Station observation')
ax.errorbar(fake_x2, test_mean, yerr=test_var, fmt='o', color=B, \
           ms=12, mec='k', capsize=5, elinewidth=3, markeredgewidth=1.5, label='GPR Prediction')
ax.errorbar(fake_x2, test_ori, yerr=test_ERROR, fmt='o', color='r', \
                ms=12, mec='k', capsize=5, elinewidth=3, markeredgewidth=1.5, label='GPM retrieval')
#ax.set_title('Comparisons of GPR prediction and GPM retrievals')
LG = ax.legend(bbox_to_anchor=(0.375, 1.025), prop={'size':14});
LG.get_frame().set_facecolor('white')
LG.get_frame().set_edgecolor('k')
LG.get_frame().set_linewidth(1.5)

ax.set_xticks(fake_x2)
ax.set_xticklabels(test_code);

fig.savefig('/Users/yingkai/Desktop/GPR_high_quality/fig3.png', dpi=250, orientation='portrait', \
           papertype='a4', format='png', bbox_inches='tight', pad_inches=0.1)
print('GPR bias-corrected:')
cal_verif_stats(test_mean, test_obs)
print('Original:')
cal_verif_stats(test_ori, test_obs)<jupyter_output>Original:
{'fraction_non_alphanumeric': 0.16087833295989246, 'fraction_numerical': 0.0345059377100605, 'mean_word_length': 5.471725471242146, 'pattern_counts': {'":': 0, '<': 27, '<?xml version=': 0, '>': 25, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '161664', 'license': 'Unlicense', 'n_tokens_mistral': 5625, 'n_tokens_neox': 5276, 'n_words': 1025, 'path': 'Gaussian_Process_Project/Statistical_bias_correction_for_satellite_precipitation_in_British_Columbia__an_application_of_Gaussian_process_regression.ipynb', 'repo_name': 'yingkaisha/Homework', 'size': '13389'}