End of Analysis

notebooks/.virtual_documents/paleoPCA.ipynb

@@ -42,6 +42,8 @@ from eofs.xarray import Eof
 import cartopy.crs as ccrs
 import matplotlib.pyplot as plt
 import nc_time_axis
+from matplotlib import gridspec
+from matplotlib.colors import Normalize
 
 
 
@@ -244,7 +246,7 @@ plt.show(fig)
 
 
 
-mgs_common = mgs.common_time().standardize()
+mgs_common = mgs.sel(time=slice(850,2000)).common_time().standardize()
 
 
 pca = mgs_common.pca()
@@ -334,7 +336,7 @@ ds
 
 
 
-ds_geo = ds.sel(lat=slice(-27,27), lon=slice(250,330))
+ds_geo_all = ds.sel(lat=slice(-27,27), lon=slice(250,330))
 
 
 
@@ -346,6 +348,25 @@ ds_geo = ds.sel(lat=slice(-27,27), lon=slice(250,330))
 
 
 
+print("The minimum time is: "+str(np.min(mgs_common.series_list[0].time)))
+print("The maximum time is: "+str(np.max(mgs_common.series_list[0].time)))
+print("The resolution is: "+str(np.mean(np.diff(mgs_common.series_list[0].time))))
+
+
+ds_geo_time = ds_geo_all.sel(time=slice("0910-01-01 00:00:00" ,"1642-12-01 00:00:00"))
+
+
+
+
+
+ds_geo = ds_geo_time.resample(time='20A').mean()
+
+ds_geo
+
+
+
+
+
 %%time
 p16O = ds_geo['PRECRC_H216Or'] + ds_geo['PRECSC_H216Os'] + ds_geo['PRECRL_H216OR'] + ds_geo['PRECSL_H216OS']
 p18O = ds_geo['PRECRC_H218Or'] + ds_geo['PRECSC_H218Os'] + ds_geo['PRECRL_H218OR'] + ds_geo['PRECSL_H218OS']
@@ -397,7 +418,6 @@ pcs = solver.pcs(npcs=3, pcscaling=1)
 
 fig, ax = plt.subplots(figsize=[20,4])
 pcs[:, 0].plot(ax=ax, linewidth=1)
-
 ax = plt.gca()
 ax.axhline(0, color='k')
 ax.set_ylim(-3, 3)
@@ -409,6 +429,47 @@ ax.set_title('PC1 Time Series', fontsize=16)
 
 
 
+# Create a figure
+
+fig = plt.figure(figsize=[20,8])
+
+# Define the GridSpec
+gs = gridspec.GridSpec(1, 2, figure=fig)
+
+# Add a geographic map in the first subplot using Cartopy
+
+ax1 = fig.add_subplot(gs[0, 0], projection=ccrs.PlateCarree(central_longitude=290))
+ax1.coastlines()  # Add coastlines to the map
+
+# Plot the model results
+norm = Normalize(vmin=-1, vmax=1)
+eof1[0].plot.contourf(ax=ax1, levels = clevs, cmap=plt.cm.RdBu_r,
+                         transform=ccrs.PlateCarree(), add_colorbar=True, norm=norm)
+ax1.set_title('EOF1 expressed as covariance', fontsize=16)
+fig.axes[1].set_ylabel('')
+fig.axes[1].set_yticks(np.arange(-1,1.2,0.2))
+
+#Now let's scatter the proxy data
+EOF = pca.eigvecs[:, 0]
+ax1.scatter(filtered_df2['lon'],filtered_df2['lat'], c =EOF, cmap=plt.cm.RdBu_r, transform=ccrs.PlateCarree(), norm=norm, s=400, edgecolor='k', linewidth=3)
+
+## Let's plot the PCS!
+PC = pca.pcs[:, 0]
+
+
+ax2 = fig.add_subplot(gs[0, 1:],)
+time_model = np.arange(910,1660,20)
+ts1 = pyleo.Series(time = time_model, value =  pcs[:, 0], time_name = 'Years', time_unit = 'CE', value_name='PC1', label = 'CESM-LEM',verbose=False)
+ts2 = pyleo.Series(time = mgs_common.series_list[0].time, value =  PC, time_name = 'Years', time_unit = 'CE', value_name='PC1', label = 'Proxy Data', verbose = False)
+ts1.plot(ax=ax2, legend = True)
+ts2.plot(ax=ax2, legend = True)
+ax2.set_ylim([-1,1])
+ax2.legend()
+# Layout adjustments and display the figure
+plt.tight_layout()
+plt.show()
+
+