<jupyter_start><jupyter_text>This series of notebooks focuses on using multiple linear regression to represent 2007 HRDPS eigenvectors as a linear combination of 2007 CANRCM eigenvectors. I split this into multiple notebooks. This notebook simply completes the PCA analysis and exports the results for later notebooks.   PCA on HRDPS data<jupyter_code>files = glob.glob('/results/forcing/atmospheric/GEM2.5/gemlam/gemlam_y2007m??d??.nc')
##list of data files, febrary 2007 - each file covers one day, with hourly data
##the grid is 266 by 256 points, with resolution 2.5km
##relevant data in this file is u_wind and v_wind data
files.sort()
print(len(files))
##data must be converted into a 2D matrix for pca analysis
##transform takes a 3D data array (time, a, b) -> (a*b, time)
##(the data grid is flattened a column using numpy.flatten)

def transform(xarr):
    arr = np.array(xarr)
    data3 = np.empty((arr.shape[1]*arr.shape[2], arr.shape[0])) 
    for i in range(arr.shape[0]):
        flattened_arr = arr[i,:,:].flatten(order = 'C')
        data3[:, i] = flattened_arr[:]; 
    return data3

##reverse does the inverse of transform
def reverse(arr, arr_orig):
    a = arr_orig.shape[1]
    b = arr_orig.shape[2]
    data4 = np.empty((arr.shape[1], a, b))
    for i in range(arr.shape[1]):
        temp = arr[:,i].reshape(a,b)
        data4[i,:,:] = temp
    ##arr_new = arr_orig.copy(data = data4) ##copies the metadata from arr_orig, and the data from data4
    return data4
## 3-hour averaged matrix

u_data3 = np.empty( (8*len(files), 266, 256)) 
v_data3 = np.empty( (8*len(files), 266, 256)) 

for i in range(len(files)):
    dayX = xr.open_dataset(files[i])
    for j in range(8):
        u_3 = np.array( dayX['u_wind'][j*3:j*3 + 3, 0:266, 0:256] )
        v_3 = np.array( dayX['v_wind'][j*3:j*3 + 3, 0:266, 0:256] )
        avg_u = u_3.mean(axis = 0) ##average of 3-hour time interval
        avg_v = v_3.mean(axis = 0)
        u_data3[8*i + j, : , : ] = avg_u[0:266, 0:256] ##adding 3-hour average to new data array
        v_data3[8*i + j, : , : ] = avg_v[0:266, 0:256]
    
u_mat3 = transform(u_data3)
v_mat3 = transform(v_data3)
c_mat3 = np.concatenate((u_mat3,v_mat3), axis=0)

del u_mat3
del v_mat3
del u_3
del v_3
del avg_u
del avg_v
del dayX
u_data3.shape
##conducting PCA analysis on combined data matrix
pca = PCA(n_components = 50) ##adjust the number of principle conponents to be calculated
PCs_h = pca.fit_transform(c_mat3)
eigvecs_h = pca.components_
fracVar_h = pca.explained_variance_ratio_
##graphing percentage of original data represented by the first n principle conponents
def graph_variance(matrix, n):
    pcaG = PCA(n_components = n) ##Number of principle conponents to show
    PCsG = pcaG.fit_transform(matrix)
    plt.plot(np.cumsum(pcaG.explained_variance_ratio_))
    plt.xlabel('number of components')
    plt.ylabel('cumulative explained variance');
    plt.show()
graph_variance(c_mat3, 30)
def get_mode(PCs, n, arr_orig): ##converts PCs (column vectors) to 2d conpoents for u and v wind
    split = int(PCs.shape[0]/2)
    mode_u = PCs[:split, n].reshape(arr_orig.shape[1], arr_orig.shape[2])
    mode_v = PCs[split:, n].reshape(arr_orig.shape[1], arr_orig.shape[2])
    return (mode_u, mode_v)

def graph_nPCs(PCs, eigvecs, n, arr_orig):
    fig, ax = plt.subplots(n, 3, figsize=(10, 3*n))
    
    ax[0, 0].set_title("u-conponent")
    ax[0, 1].set_title("v-component")
    ax[0, 2].set_title("time-loadings")
    
    for i in range(n):
        mode_u, mode_v = get_mode(PCs, i, arr_orig)
        colors = ax[i, 0].pcolormesh(mode_u, cmap = 'bwr')
        fig.colorbar(colors, ax = ax[i,0])
        colors = ax[i, 1].pcolormesh(mode_v, cmap = 'bwr')
        fig.colorbar(colors, ax = ax[i,1])
        ax[i, 2].plot(eigvecs[i])

    plt.tight_layout()
    plt.show()
graph_nPCs(PCs_h, eigvecs_h, 5, u_data3)<jupyter_output><empty_output><jupyter_text>Pressure PCA on CANRCM data<jupyter_code>p3 = '/home/arandhawa/canrcm_sea_level_pressure_2007.nc'
d3 = xr.open_dataset(p3)
print(d3)
pressure = np.array(d3['psl'][16:,140:165,60:85])
pres_mat = transform(pressure)
#plots of original data at two timestamps
fig, axs = plt.subplots(1, 2, figsize=(10, 5))
axs[0].pcolormesh(pressure[0,:,:])
axs[1].pcolormesh(pressure[500,:,:])
fig.tight_layout()
##conducting PCA analysis on combined data matrix
pca = PCA(n_components = 50) ##adjust the number of principle conponents to be calculated
PCs = pca.fit_transform(pres_mat)
eigvecs = pca.components_
fracVar = pca.explained_variance_ratio_
def graph_pressure(PCs, eigvecs, n, arr_orig):
    fig, ax = plt.subplots(n, 2, figsize=(7, 3*n))
    
    ax[0, 0].set_title("Modes")
    ax[0, 1].set_title("Time-loadings")
    
    for i in range(n):
        mode = PCs[:,i].reshape(arr_orig.shape[1], arr_orig.shape[2])
        colors = ax[i, 0].pcolormesh(mode, cmap = 'bwr')
        fig.colorbar(colors, ax = ax[i,0])
    
        ax[i, 1].plot(eigvecs[i])

    plt.tight_layout()
    plt.show()
graph_variance(pres_mat, 30)
graph_pressure(PCs, eigvecs, 5, pressure)
d3.close()<jupyter_output><empty_output><jupyter_text>CANRCM Winds<jupyter_code>p1 = '/home/arandhawa/canrcm_2007_uas.nc'
p2 = '/home/arandhawa/canrcm_2007_vas.nc'

d1 = xr.open_dataset(p1)
d2 = xr.open_dataset(p2)

A2 = d1['uas'][16:,140:165,60:85]
B2 = d2['vas'][16:,140:165,60:85]
#plots of original data at two timestamps
fig, axs = plt.subplots(1, 2, figsize=(10, 5))
axs[0].pcolormesh(np.array(A2)[0,:,:])
axs[1].pcolormesh(np.array(A2)[0,:,:])
fig.tight_layout()
A_mat2 = transform(A2)
B_mat2 = transform(B2)
C_mat2 = np.concatenate((A_mat2,B_mat2), axis=0) 

##conducting PCA analysis on combined data matrix
pca2 = PCA(n_components = 30) ##adjust the number of principle conponents to be calculated
PCs2 = pca2.fit_transform(C_mat2)
eigvecs2 = pca2.components_
fracVar2 = pca2.explained_variance_ratio_

##reconstructing data using the chosen number of principle conponents
C_new2 = pca2.inverse_transform(pca2.transform(C_mat2))
split2 = int( C_new2.shape[0]/2 )
A_new2 = reverse(C_new2[:split2,:], A2) ##reconstructing u_winds from 10 PCs
B_new2 = reverse(C_new2[split2:,:], B2) ##reconstructing v_winds from 10 PCs
##graphing percentage of original data represented by the first n principle conponents
graph_variance(C_mat2, 30)
graph_nPCs(PCs2, eigvecs2, 5, A2)<jupyter_output><empty_output><jupyter_text>Exports<jupyter_code>##HRDPS exports
np.savetxt("eigvecs_hrdps_annual.csv", eigvecs_h, delimiter=",")
np.savetxt("PCs_hrdps_annual.csv", PCs_h, delimiter=",")
np.savetxt("raw_data_3h.csv", c_mat3, delimiter=",")
##CANRCM pressure exports
np.savetxt("eigvecs_psl.csv", eigvecs, delimiter=",") 
np.savetxt("PCs_psl.csv", PCs, delimiter=",")
np.savetxt("raw_data_psl.csv", pres_mat, delimiter=",")
##CANRCM wind exports
np.savetxt("eigvecs_winds.csv", eigvecs2, delimiter=",") 
np.savetxt("PCs_winds.csv", PCs2, delimiter=",")
np.savetxt("raw_data_winds.csv", C_mat2, delimiter=",")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14531138153185397, 'fraction_numerical': 0.04409448818897638, 'mean_word_length': 4.702857142857143, 'pattern_counts': {'":': 0, '<': 17, '<?xml version=': 0, '>': 18, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '163319', 'license': 'Apache-2.0', 'n_tokens_mistral': 3007, 'n_tokens_neox': 2792, 'n_words': 633, 'path': 'notebooks/reconstructingDataP1.ipynb', 'repo_name': 'UBC-MOAD/analysis-armaan', 'size': '6985'}