import cartopy from cartopy import config import cartopy.feature as cfeature import cartopy.crs as ccrs from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER import numpy as np import datetime as dt import pandas as pd import matplotlib import matplotlib.pyplot as plt import matplotlib.lines as mlines import matplotlib.patches as mpatches from matplotlib import cm from matplotlib.colors import ListedColormap, LinearSegmentedColormap import xarray as xr import os import cfgrib import glob from sklearn.preprocessing import MinMaxScaler from sklearn.model_selection import train_test_split from sklearn.metrics import r2_score import xesmf as xe from xskillscore import rmse from scipy.stats import pearsonr,bootstrap from scipy.signal import correlate2d from data_utils import Reshape_arrays,geographical_boxes from processing_outputs import open_data_files_graphcast,opening_data_files_pangu,weekly_mean_operation,RegriddingData,\ ChoosingWeek,anomaly_correlation_coefficient,calling_climatology,removing_climatology,processing_hindcasts print('-----> Choose variable : [msl,u10,v10,t2m,tp] variables for surface fields') considered_variable = '{}'.format(input()) print('-----> Choose how many weeks ahead, i.e 1, 2, 3 or 4 weeks ahead') week = int(input())-1 remove_seasonality = True VariablesDict = {'u10':'10m_u_component_of_wind','v10':'10m_v_component_of_wind','t2m':'2m_temperature','tp':'total_precipitation'} VariablesDict_Graphcast = {'u10':'10m_u_component_of_wind','v10':'10m_v_component_of_wind','t2m':'2m_temperature','tp':'total_precipitation_6hr'} #settings u10 if chosen variable is 'tp' just as dummy data. Pangu Weather doesnt predict 'tp' #considered_variable_graphcast = considered_variable if considered_variable != 'tp' else 'tp' considered_variable_pangu = 'u10' if considered_variable == 'tp' else considered_variable assert considered_variable_pangu in ['msl','u10','v10','t2m','tp'], 'Wrong ! Enter correct variable!' S2S_lat,S2S_lon = np.load('/home/vsansi01/S2S_lat.npy'),np.load('/home/vsansi01/S2S_lon.npy') #GraphFT Observations150km_GraphFT = xr.open_dataset(f'/home/vsansi01/observations_predictions_fields/Observations_GraphFT_{considered_variable}_{week+1}.nc').to_array() weekly_FineTunedGraphForecasts = xr.open_dataset(f'/home/vsansi01/observations_predictions_fields/Predictions_GraphFT_{considered_variable}_{week+1}.nc').to_array() #S2S forecasts deterministic Observations150km_S2S = xr.open_dataset(f'/home/vsansi01/observations_predictions_fields/Observations_S2S_{considered_variable}_{week+1}.nc').to_array() WeeklyS2S = xr.open_dataset(f'/home/vsansi01/observations_predictions_fields/Predictions_S2S_{considered_variable}_{week+1}.nc').to_array() #S2S forecasts probabilistic #Observations150km_S2S_ensemble = xr.open_dataset(f'/home/vsansi01/observations_predictions_fields/Observations_S2S_ensemble_{considered_variable}_{week+1}.nc').to_array() #WeeklyS2S_ensemble = xr.open_dataset(f'/home/vsansi01/observations_predictions_fields/Predictions_S2S_ensemble_{considered_variable}_{week+1}.nc').to_array() #Graphcast OP Observations150km_GraphOP = xr.open_dataset(f'/home/vsansi01/observations_predictions_fields/Observations_GraphOP_{considered_variable}_{week+1}.nc').to_array() Weekly_OperationalGraphForecasts = xr.open_dataset(f'/home/vsansi01/observations_predictions_fields/Predictions_GraphOP_{considered_variable}_{week+1}.nc').to_array() #Pangu weather forecasts Observations150km_Pangu = xr.open_dataset(f'/home/vsansi01/observations_predictions_fields/Observations_Pangu_{considered_variable}_{week+1}.nc').to_array() weekly_PanguForecasts= xr.open_dataset(f'/home/vsansi01/observations_predictions_fields/Predictions_Pangu_{considered_variable}_{week+1}.nc').to_array() #lattitude weight when computing R correlation coefficient CosPhi = np.array([[np.cos(lat*np.pi/180) for lon in S2S_lon] for lat in S2S_lat]) #computing R spatial maps RCoeff_GraphFT = np.squeeze(xr.corr(Observations150km_GraphFT,weekly_FineTunedGraphForecasts,dim='time')) RCoeff_GraphOP = np.squeeze(xr.corr(Observations150km_GraphOP,Weekly_OperationalGraphForecasts,dim='time')) RCoeff_Pangu = np.squeeze(xr.corr(Observations150km_Pangu,weekly_PanguForecasts,dim='time')) RCoeff_S2S = np.squeeze(xr.corr(Observations150km_S2S,WeeklyS2S,dim='time')) #RCoeff_S2S_ensemble = np.squeeze(xr.corr(np.sqrt(CosPhi)*Observations150km_S2S_ensemble,np.sqrt(CosPhi)*WeeklyS2S_ensemble,dim='time')) #compute the RCoeff over RCoeff_xt_GraphFT = [pearsonr(observed_anomaly.reshape(-1),forecasted_anomaly.reshape(-1)).statistic\ for observed_anomaly,forecasted_anomaly in zip(Reshape_arrays(Observations150km_GraphFT),Reshape_arrays(weekly_FineTunedGraphForecasts))] RCoeff_xt_GraphOP = [pearsonr(observed_anomaly.reshape(-1),forecasted_anomaly.reshape(-1)).statistic\ for observed_anomaly,forecasted_anomaly in zip(Reshape_arrays(Observations150km_GraphOP),Reshape_arrays(Weekly_OperationalGraphForecasts))] RCoeff_xt_Pangu = [pearsonr(observed_anomaly.reshape(-1),forecasted_anomaly.reshape(-1)).statistic\ for observed_anomaly,forecasted_anomaly in zip(Reshape_arrays(Observations150km_Pangu),Reshape_arrays(weekly_PanguForecasts))] RCoeff_xt_S2S = [pearsonr(observed_anomaly.reshape(-1),forecasted_anomaly.reshape(-1)).statistic\ for observed_anomaly,forecasted_anomaly in zip(Reshape_arrays(Observations150km_S2S),Reshape_arrays(WeeklyS2S))] np.save(f'/home/vsansi01/ACC_R_values/RCoeff_xt_GraphFT_{considered_variable}_{week+1}week.npy',np.array(RCoeff_xt_GraphFT)) np.save(f'/home/vsansi01/ACC_R_values/RCoeff_xt_GraphOP_{considered_variable}_{week+1}week.npy',RCoeff_xt_GraphOP) np.save(f'/home/vsansi01/ACC_R_values/RCoeff_xt_S2S_{considered_variable}_{week+1}week.npy',np.array(RCoeff_xt_S2S)) #if remove_seasonality: # RCoeff_GraphFT.to_netcdf('/home/vsansi01/ACC_R_values/R_{}_{}_week{}.nc'.format('GraphFT',considered_variable,week+1)) # RCoeff_S2S.to_netcdf('/home/vsansi01/ACC_R_values/R_{}_{}_week{}.nc'.format('S2S',considered_variable,week+1)) # if week != 3: # RCoeff_GraphOP.to_netcdf('/home/vsansi01/ACC_R_values/R_{}_{}_week{}.nc'.format('GraphOP',considered_variable,week+1)) # RCoeff_Pangu.to_netcdf('/home/vsansi01/ACC_R_values/R_{}_{}_week{}.nc'.format('Pangu',considered_variable,week+1)) RCoeff_GraphFT_global = np.mean(RCoeff_GraphFT) RCoeff_S2S_global = np.mean(RCoeff_S2S) RCoeff_GraphOP_global = np.mean(RCoeff_GraphOP) RCoeff_Pangu_global = np.mean(RCoeff_Pangu) print('RCoeff GraphFT' , RCoeff_GraphFT_global.values) print('RCoeff GraphOP' , RCoeff_GraphOP_global.values) print('RCoeff Pangu' , RCoeff_Pangu_global.values) print('RCoeff S2S' , RCoeff_S2S_global.values) RCoeff_GraphFT_OP = RCoeff_GraphFT.mean('longitude') - RCoeff_GraphOP.mean('longitude') RCoeff_GraphFT_S2S = RCoeff_GraphFT.mean('longitude') - RCoeff_S2S.mean('longitude') RCoeff_GraphFT_OP_xt = np.array(RCoeff_xt_GraphFT) - np.array(RCoeff_xt_GraphOP) RCoeff_GraphFT_S2S_xt = np.array(RCoeff_xt_GraphFT) - np.array(RCoeff_xt_S2S) plt.figure()#figsize=(15,15)) plt.subplot(2,1,1) plt.plot(S2S_lat,RCoeff_xt_GraphFT,label='GraphFT') if week != 3: plt.plot(S2S_lat,RCoeff_xt_GraphOP,label='GraphOP') plt.plot(S2S_lat,RCoeff_xt_S2S,label='S2S') plt.legend(loc='upper right') plt.ylabel('Zonal mean TCC') plt.xlabel('Latitudes') plt.ylim(-0.35,1.2) plt.subplot(2,1,2) xx = np.linspace(-90,90,100) yy = [0 for x in xx] if week != 3: plt.plot(S2S_lat,RCoeff_GraphFT_OP_xt,label='GraphFT - GraphOP') plt.plot(S2S_lat,RCoeff_GraphFT_S2S_xt,label='GraphFT - S2S') plt.plot(xx,yy,'r') plt.ylabel('Difference of zonal mean TCC') plt.xlabel('Latitudes') plt.ylim(-0.7,1.2) plt.legend(loc='upper right')#bbox_to_anchor=(0.85, 0.5, 0.5, 0.5)) plt.tight_layout() #if remove_seasonality: # plt.savefig('/home/vsansi01/tcc_latitudes_{}_{}WeeksAhead_rf.jpeg'.format(considered_variable,week+1),dpi=300) #computing RMSE RMSE_GraphFT = np.squeeze(rmse(Observations150km_GraphFT,weekly_FineTunedGraphForecasts,dim='time')) RMSE_GraphOP = np.squeeze(rmse(Observations150km_GraphOP,Weekly_OperationalGraphForecasts,dim='time')) RMSE_Pangu = np.squeeze(rmse(Observations150km_Pangu,weekly_PanguForecasts,dim='time')) RMSE_S2S = np.squeeze(rmse(Observations150km_S2S,WeeklyS2S,dim='time')) RMSE_GraphFT_global = np.mean(RMSE_GraphFT) RMSE_GraphOP_global = np.mean(RMSE_GraphOP) RMSE_Pangu_global = np.mean(RMSE_Pangu) RMSE_S2S_global = np.mean(RMSE_S2S) print('RMSE GraphFT' , RMSE_GraphFT_global.values) print('RMSE GraphOP' , RMSE_GraphOP_global.values) print('RMSE Pangu' , RMSE_Pangu_global.values) print('RMSE S2S' , RMSE_S2S_global.values) RMSE_GraphFT_OP = RMSE_GraphFT.mean('longitude') - RMSE_GraphOP.mean('longitude') RMSE_GraphFT_S2S = RMSE_GraphFT.mean('longitude') - RMSE_S2S.mean('longitude') plt.figure() plt.subplot(2,1,1) plt.plot(S2S_lat,np.squeeze(RMSE_GraphFT.mean('longitude')),label='GraphFT') plt.plot(S2S_lat,np.squeeze(RMSE_GraphOP.mean('longitude')),label='GraphOP') plt.plot(S2S_lat,np.squeeze(RMSE_S2S.mean('longitude')),label='S2S') plt.legend(loc='upper right') plt.ylabel('Zonal mean RMSE') plt.xlabel('Latitudes') plt.ylim(0,10) plt.subplot(2,1,2) xx = np.linspace(-90,90,100) yy = [0 for x in xx] plt.plot(S2S_lat,np.squeeze(RMSE_GraphFT_OP),label='GraphFT - GraphOP') plt.plot(S2S_lat,np.squeeze(RMSE_GraphFT_S2S),label='GraphFT - S2S') plt.plot(xx,yy,'r') plt.ylabel('Difference of zonal mean RMSE') plt.xlabel('Latitudes') plt.ylim(-5,5) plt.legend(loc='upper right')#bbox_to_anchor=(0.85, 0.5, 0.5, 0.5)) plt.tight_layout() #if remove_seasonality: # plt.savefig('/home/vsansi01/rmse_latitudes_{}_{}WeeksAhead_rf.jpeg'.format(considered_variable,week+1),dpi=300) #if not remove_seasonality: #computing ACC # ACC_GraphFT,std_GraphFT = anomaly_correlation_coefficient(Observations150km_GraphFT,weekly_FineTunedGraphForecasts, Dates_Observations_GraphFT,S2S_nc,considered_variable,week,'GraphFT') #ACC_GraphOP,std_GraphOP = anomaly_correlation_coefficient(Observations150km_GraphOP,Weekly_OperationalGraphForecasts, Dates_Observations_GraphOP,S2S_nc,considered_variable,week,False,'GraphOP') #ACC_Pangu,std_Pangu = anomaly_correlation_coefficient(Observations150km_Pangu,weekly_PanguForecasts, Dates_Observations_Pangu,S2S_nc,considered_variable,week,False,'Pangu') #ACC_S2S,std_S2S = anomaly_correlation_coefficient(Observations150km_S2S,WeeklyS2S,Dates_Observations_S2S,S2S_nc,considered_variable,week,'S2S_ensemble_mean') #print('ACC_GraphFT',ACC_GraphFT, 'std_GraphFT',std_GraphFT) #print('ACC_GraphOP',ACC_GraphOP,'std_GraphOP',std_GraphOP) #print('ACC_Pangu',ACC_Pangu,'std_Pangu',std_Pangu) #print('ACC_S2S',ACC_S2S,'std_S2S',std_S2S) #ACC_GraphFT_OP = np.array(ACC_GraphFT) - np.array(ACC_GraphOP) #ACC_GraphFT_S2S = np.array(ACC_GraphFT) - np.array(ACC_S2S) #plt.figure() #plt.subplot(2,1,1) #plt.plot(S2S_lat,ACC_GraphFT,label='GraphFT') #plt.plot(S2S_lat,ACC_GraphOP,label='GraphOP') #plt.plot(S2S_lat,ACC_S2S,label='S2S') #plt.legend(loc='upper right') #plt.ylabel('Zonal mean TCC') #plt.xlabel('Latitudes') #plt.subplot(2,1,2) #xx = np.linspace(-90,90,100) #yy = [0 for x in xx] #plt.ylabel('Difference of zonal mean TCC') #plt.xlabel('Latitudes') #plt.ylim(-0.7,1.2) #plt.plot(xx,yy,'r') #plt.plot(S2S_lat,ACC_GraphFT_OP,label='GraphFT - GraphOP') #plt.plot(S2S_lat,ACC_GraphFT_S2S,label='GraphFT - S2S') #plt.legend() #plt.tight_layout() #plt.savefig('/home/vsansi01/acc_latitudes_{}_{}WeeksAhead_cf.jpeg'.format(considered_variable,week+1),dpi=300) if week in [0,1]: #hawaii archipelagos geographical_boxes(RCoeff_GraphFT,min_lon=195,min_lat=14,max_lon=212,max_lat=25,name=f'GraphFT_{considered_variable}_{week+1}week_hawaii') geographical_boxes(RCoeff_GraphOP,min_lon=195,min_lat=14,max_lon=212,max_lat=25,name=f'GraphOP_{considered_variable}_{week+1}week_hawaii') geographical_boxes(RCoeff_S2S,min_lon=195,min_lat=14,max_lon=212,max_lat=25,name=f'S2S_{considered_variable}_{week+1}week_hawaii') #geographical_boxes(RCoeff_S2S_ensemble,min_lon=195,min_lat=14,max_lon=212,max_lat=25,name=f'S2S_ensemble_{considered_variable}_{week+1}week_hawaii') #North America - Canada geographical_boxes(RCoeff_GraphFT,min_lon=186,min_lat=25,max_lon=281,max_lat=63,name=f'GraphFT_{considered_variable}_{week+1}week_north_america') geographical_boxes(RCoeff_GraphOP,min_lon=186,min_lat=25,max_lon=281,max_lat=63,name=f'GraphOP_{considered_variable}_{week+1}week_north_america') geographical_boxes(RCoeff_S2S,min_lon=186,min_lat=25,max_lon=281,max_lat=63,name=f'S2S_{considered_variable}_{week+1}week_north_america') #geographical_boxes(RCoeff_S2S_ensemble,min_lon=186,min_lat=25,max_lon=281,max_lat=63,name=f'S2S_ensemble_{considered_variable}_{week+1}week_north_america') #Nothern pole geographical_boxes(RCoeff_GraphFT,min_lon=0,min_lat=60,max_lon=360,max_lat=90,name=f'GraphFT_{considered_variable}_{week+1}week_north_pole') geographical_boxes(RCoeff_GraphOP,min_lon=0,min_lat=60,max_lon=360,max_lat=90,name=f'GraphOP_{considered_variable}_{week+1}week_north_pole') geographical_boxes(RCoeff_S2S,min_lon=0,min_lat=60,max_lon=360,max_lat=90,name=f'S2S_{considered_variable}_{week+1}week_north_pole') #geographical_boxes(RCoeff_S2S_ensemble,min_lon=0,min_lat=60,max_lon=360,max_lat=90,name=f'S2S_ensemble_{considered_variable}_{week+1}week_north_pole') #South pole geographical_boxes(RCoeff_GraphFT,min_lon=0,min_lat=-90,max_lon=360,max_lat=-60,name=f'GraphFT_{considered_variable}_{week+1}week_south_pole') geographical_boxes(RCoeff_GraphOP,min_lon=0,min_lat=-90,max_lon=360,max_lat=-60,name=f'GraphOP_{considered_variable}_{week+1}week_south_pole') geographical_boxes(RCoeff_S2S,min_lon=0,min_lat=-90,max_lon=360,max_lat=-60,name=f'S2S_{considered_variable}_{week+1}week_south_pole') #geographical_boxes(RCoeff_S2S_ensemble,min_lon=0,min_lat=-90,max_lon=360,max_lat=-60,name=f'S2S_ensemble_{considered_variable}_{week+1}week_south_pole') #Indian subcontinent geographical_boxes(RCoeff_GraphFT,min_lon=66,min_lat=5,max_lon=87,max_lat=28,name=f'GraphFT_{considered_variable}_{week+1}week_india') geographical_boxes(RCoeff_GraphOP,min_lon=66,min_lat=5,max_lon=87,max_lat=28,name=f'GraphOP_{considered_variable}_{week+1}week_india') geographical_boxes(RCoeff_S2S,min_lon=66,min_lat=5,max_lon=87,max_lat=28,name=f'S2S_{considered_variable}_{week+1}week_india') #geographical_boxes(RCoeff_S2S_ensemble,min_lon=66,min_lat=5,max_lon=87,max_lat=28,name=f'S2S_ensemble_{considered_variable}_{week+1}week_india') #Tibetan plateau geographical_boxes(RCoeff_GraphFT,83,28,91,34,f'GraphFT_{considered_variable}_{week+1}week_tibetan_plateau') geographical_boxes(RCoeff_GraphOP,83,28,91,34,f'GraphOP_{considered_variable}_{week+1}week_tibetan_plateau') geographical_boxes(RCoeff_S2S,83,28,91,34,f'S2S_{considered_variable}_{week+1}week_tibetan_plateau') #geographical_boxes(RCoeff_S2S_ensemble,83,28,91,34,f'S2S_ensemble_{considered_variable}_{week+1}week_tibetan_plateau') #Northeastern Asia - Japan geographical_boxes(RCoeff_GraphFT,90,29,150,70,f'GraphFT_{considered_variable}_{week+1}week_northeastern_asia') geographical_boxes(RCoeff_GraphOP,90,29,150,70,f'GraphOP_{considered_variable}_{week+1}week_northeastern_asia') geographical_boxes(RCoeff_S2S,90,29,150,70,f'S2S_{considered_variable}_{week+1}week_northeastern_asia') #geographical_boxes(RCoeff_S2S_ensemble,90,29,150,70,f'S2S_ensemble_{considered_variable}_{week+1}week_northeastern_asia') #North - Atlantic geographical_boxes(RCoeff_GraphFT,299,14,351,50,f'GraphFT_{considered_variable}_{week+1}week_north_atlantic') geographical_boxes(RCoeff_GraphOP,299,14,351,50,f'GraphOP_{considered_variable}_{week+1}week_north_atlantic') geographical_boxes(RCoeff_S2S,299,14,351,50,f'S2S_{considered_variable}_{week+1}week_north_atlantic') #geographical_boxes(RCoeff_S2S_ensemble,299,14,351,50,f'S2S_ensemble{considered_variable}_{week+1}week_north_atlantic') #Europe !!! coords = [longitude - 9 for longitude in RCoeff_GraphFT.longitude] geographical_boxes(RCoeff_GraphFT.assign_coords({'longitude':coords}),0,35,47,67,f'GraphFT_{considered_variable}_{week+1}week_europe') geographical_boxes(RCoeff_GraphOP.assign_coords({'longitude':coords}),0,35,47,67,f'GraphOP_{considered_variable}_{week+1}week_europe') geographical_boxes(RCoeff_S2S.assign_coords({'longitude':coords}),0,35,47,67,f'S2S_{considered_variable}_{week+1}week_europe') #geographical_boxes(RCoeff_S2S_ensemble.assign_coords({'longitude':coords}),0,35,47,67,f'S2S_ensemble_{considered_variable}_{week+1}week_europe') #South America geographical_boxes(RCoeff_GraphFT,279,-34,314,8,f'GraphFT_{considered_variable}_{week+1}week_south_america') geographical_boxes(RCoeff_GraphOP,279,-34,314,8,f'GraphOP_{considered_variable}_{week+1}week_south_america') geographical_boxes(RCoeff_S2S,279,-34,314,8,f'S2S_{considered_variable}_{week+1}week_south_america') #geographical_boxes(RCoeff_S2S_ensemble,279,-34,314,8,f'S2S_ensemble_{considered_variable}_{week+1}week_south_america') #Australia geographical_boxes(RCoeff_GraphFT,117,-38,151,-18,f'GraphFT_{considered_variable}_{week+1}week_australia') geographical_boxes(RCoeff_GraphOP,117,-38,151,-18,f'GraphOP_{considered_variable}_{week+1}week_australia') geographical_boxes(RCoeff_S2S,117,-38,151,-18,f'S2S_{considered_variable}_{week+1}week_australia') #geographical_boxes(RCoeff_S2S_ensemble,117,-38,151,-18,f'S2S_ensemble_{considered_variable}_{week+1}week_australia') #Central Asia geographical_boxes(RCoeff_GraphFT,72,35,103,50,f'GraphFT_{considered_variable}_{week+1}week_central_asia') geographical_boxes(RCoeff_GraphOP,72,35,103,50,f'GraphOP_{considered_variable}_{week+1}week_central_asia') geographical_boxes(RCoeff_S2S,72,35,103,50,f'S2S_{considered_variable}_{week+1}week_central_asia') #geographical_boxes(RCoeff_S2S_ensemble,72,35,103,50,f'S2S_ensemble_{considered_variable}_{week+1}week_central_asia') #Northern Hemisphere geographical_boxes(RCoeff_GraphFT,0,30,360,90,f'GraphFT_{considered_variable}_{week+1}week_north_hemisphere') geographical_boxes(RCoeff_GraphOP,0,30,360,90,f'GraphOP_{considered_variable}_{week+1}week_north_hemisphere') geographical_boxes(RCoeff_S2S,0,30,360,90,f'S2S_{considered_variable}_{week+1}week_north_hemisphere') #geographical_boxes(RCoeff_S2S_ensemble,0,30,360,90,f'S2S_ensemble_{considered_variable}_{week+1}week_north_hemisphere') #Equatorial band geographical_boxes(RCoeff_GraphFT,0,-30,360,30,f'GraphFT_{considered_variable}_{week+1}week_equatorial_band') geographical_boxes(RCoeff_GraphOP,0,-30,360,30,f'GraphOP_{considered_variable}_{week+1}week_equatorial_band') geographical_boxes(RCoeff_S2S,0,-30,360,30,f'S2S_{considered_variable}_{week+1}week_equatorial_band') #geographical_boxes(RCoeff_S2S_ensemble,0,-30,360,30,f'S2S_ensemble_{considered_variable}_{week+1}week_equatorial_band') #Southern Hemisphere geographical_boxes(RCoeff_GraphFT,0,-90,360,-30,f'GraphFT_{considered_variable}_{week+1}week_south_hemisphere') geographical_boxes(RCoeff_GraphOP,0,-90,360,-30,f'GraphOP_{considered_variable}_{week+1}week_south_hemisphere') geographical_boxes(RCoeff_S2S,0,-90,360,-30,f'S2S_{considered_variable}_{week+1}week_south_hemisphere') #geographical_boxes(RCoeff_S2S_ensemble,0,-90,360,-30,f'S2S_ensemble_{considered_variable}_{week+1}week_south_hemisphere') #French Polynesia geographical_boxes(RCoeff_GraphFT,180,-35,240,0,f'GraphFT_{considered_variable}_{week+1}week_french_polynesia') geographical_boxes(RCoeff_GraphOP,180,-35,240,0,f'GraphOP_{considered_variable}_{week+1}week_french_polynesia') geographical_boxes(RCoeff_S2S,180,-35,240,0,f'S2S_{considered_variable}_{week+1}week_french_polynesia') #geographical_boxes(RCoeff_S2S_ensemble,180,-35,240,0,f'S2S_ensemble_{considered_variable}_{week+1}week_french_polynesia') #box2_1 coordinates min_lon = 180 min_lat = 0 max_lon = 200 max_lat = -35 RCoeff_GraphFT_Box21 = RCoeff_GraphFT.sel(latitude=slice(min_lat,max_lat), longitude=slice(min_lon,max_lon)).mean().values RCoeff_GraphOP_Box21 = RCoeff_GraphOP.sel(latitude=slice(min_lat,max_lat), longitude=slice(min_lon,max_lon)).mean().values RCoeff_Pangu_Box21 = RCoeff_Pangu.sel(latitude=slice(min_lat,max_lat), longitude=slice(min_lon,max_lon)).mean().values RCoeff_S2S_Box21 = RCoeff_S2S.sel(latitude=slice(min_lat,max_lat), longitude=slice(min_lon,max_lon)).mean().values np.save(f'/home/vsansi01/R_boxes/RCoeff_GraphFT_Box21_{considered_variable}_{week+1}.npy',RCoeff_GraphFT_Box21) np.save(f'/home/vsansi01/R_boxes/RCoeff_GraphOP_Box21_{considered_variable}_{week+1}.npy',RCoeff_GraphOP_Box21) np.save(f'/home/vsansi01/R_boxes/RCoeff_S2S_Box21_{considered_variable}_{week+1}.npy',RCoeff_S2S_Box21) #box2_2 coordinates min_lon = 200 min_lat = 0 max_lon = 240 max_lat = -12 RCoeff_GraphFT_Box22 = RCoeff_GraphFT.sel(latitude=slice(min_lat,max_lat), longitude=slice(min_lon,max_lon)).mean().values RCoeff_GraphOP_Box22 = RCoeff_GraphOP.sel(latitude=slice(min_lat,max_lat), longitude=slice(min_lon,max_lon)).mean().values RCoeff_Pangu_Box22 = RCoeff_Pangu.sel(latitude=slice(min_lat,max_lat), longitude=slice(min_lon,max_lon)).mean().values RCoeff_S2S_Box22 = RCoeff_S2S.sel(latitude=slice(min_lat,max_lat), longitude=slice(min_lon,max_lon)).mean().values np.save(f'/home/vsansi01/R_boxes/RCoeff_GraphFT_Box22_{considered_variable}_{week+1}.npy',RCoeff_GraphFT_Box22) np.save(f'/home/vsansi01/R_boxes/RCoeff_GraphOP_Box22_{considered_variable}_{week+1}.npy',RCoeff_GraphOP_Box22) np.save(f'/home/vsansi01/R_boxes/RCoeff_S2S_Box22_{considered_variable}_{week+1}.npy',RCoeff_S2S_Box22) #box2_3 coordinates min_lon = 200 min_lat = -12 max_lon = 240 max_lat = -20 RCoeff_GraphFT_Box23 = RCoeff_GraphFT.sel(latitude=slice(min_lat,max_lat), longitude=slice(min_lon,max_lon)).mean().values RCoeff_GraphOP_Box23 = RCoeff_GraphOP.sel(latitude=slice(min_lat,max_lat), longitude=slice(min_lon,max_lon)).mean().values RCoeff_Pangu_Box23 = RCoeff_Pangu.sel(latitude=slice(min_lat,max_lat), longitude=slice(min_lon,max_lon)).mean().values RCoeff_S2S_Box23 = RCoeff_S2S.sel(latitude=slice(min_lat,max_lat), longitude=slice(min_lon,max_lon)).mean().values np.save(f'/home/vsansi01/R_boxes/RCoeff_GraphFT_Box23_{considered_variable}_{week+1}.npy',RCoeff_GraphFT_Box23) np.save(f'/home/vsansi01/R_boxes/RCoeff_GraphOP_Box23_{considered_variable}_{week+1}.npy',RCoeff_GraphOP_Box23) np.save(f'/home/vsansi01/R_boxes/RCoeff_S2S_Box23_{considered_variable}_{week+1}.npy',RCoeff_S2S_Box23) #box2_4 coordinates min_lon = 180 min_lat = -20 max_lon = 240 max_lat = -35 RCoeff_GraphFT_Box24 = RCoeff_GraphFT.sel(latitude=slice(min_lat,max_lat), longitude=slice(min_lon,max_lon)).mean().values RCoeff_GraphOP_Box24 = RCoeff_GraphOP.sel(latitude=slice(min_lat,max_lat), longitude=slice(min_lon,max_lon)).mean().values RCoeff_Pangu_Box24 = RCoeff_Pangu.sel(latitude=slice(min_lat,max_lat), longitude=slice(min_lon,max_lon)).mean().values RCoeff_S2S_Box24 = RCoeff_S2S.sel(latitude=slice(min_lat,max_lat), longitude=slice(min_lon,max_lon)).mean().values np.save(f'/home/vsansi01/R_boxes/RCoeff_GraphFT_Box24_{considered_variable}_{week+1}.npy',RCoeff_GraphFT_Box24) np.save(f'/home/vsansi01/R_boxes/RCoeff_GraphOP_Box24_{considered_variable}_{week+1}.npy',RCoeff_GraphOP_Box24) np.save(f'/home/vsansi01/R_boxes/RCoeff_S2S_Box24_{considered_variable}_{week+1}.npy',RCoeff_S2S_Box24) if week in [2,3]: #hawaii archipelagos geographical_boxes(RCoeff_GraphFT,min_lon=195,min_lat=14,max_lon=212,max_lat=25,name=f'GraphFT_{considered_variable}_{week+1}week_hawaii') geographical_boxes(RCoeff_GraphOP,min_lon=195,min_lat=14,max_lon=212,max_lat=25,name=f'GraphOP_{considered_variable}_{week+1}week_hawaii') geographical_boxes(RCoeff_S2S,min_lon=195,min_lat=14,max_lon=212,max_lat=25,name=f'S2S_{considered_variable}_{week+1}week_hawaii') #geographical_boxes(RCoeff_S2S_ensemble,min_lon=195,min_lat=14,max_lon=212,max_lat=25,name=f'S2S_ensemble_{considered_variable}_{week+1}week_hawaii') #North America - Canada geographical_boxes(RCoeff_GraphFT,min_lon=186,min_lat=25,max_lon=281,max_lat=63,name=f'GraphFT_{considered_variable}_{week+1}week_north_america') geographical_boxes(RCoeff_GraphOP,min_lon=186,min_lat=25,max_lon=281,max_lat=63,name=f'GraphOP_{considered_variable}_{week+1}week_north_america') geographical_boxes(RCoeff_S2S,min_lon=186,min_lat=25,max_lon=281,max_lat=63,name=f'S2S_{considered_variable}_{week+1}week_north_america') #geographical_boxes(RCoeff_S2S_ensemble,min_lon=186,min_lat=25,max_lon=281,max_lat=63,name=f'S2S_ensemble_{considered_variable}_{week+1}week_north_america') #Nothern pole geographical_boxes(RCoeff_GraphFT,min_lon=0,min_lat=60,max_lon=360,max_lat=90,name=f'GraphFT_{considered_variable}_{week+1}week_north_pole') geographical_boxes(RCoeff_GraphOP,min_lon=0,min_lat=60,max_lon=360,max_lat=90,name=f'GraphOP_{considered_variable}_{week+1}week_north_pole') geographical_boxes(RCoeff_S2S,min_lon=0,min_lat=60,max_lon=360,max_lat=90,name=f'S2S_{considered_variable}_{week+1}week_north_pole') #geographical_boxes(RCoeff_S2S_ensemble,min_lon=0,min_lat=60,max_lon=360,max_lat=90,name=f'S2S_ensemble_{considered_variable}_{week+1}week_north_pole') #South pole geographical_boxes(RCoeff_GraphFT,min_lon=0,min_lat=-90,max_lon=360,max_lat=-60,name=f'GraphFT_{considered_variable}_{week+1}week_south_pole') geographical_boxes(RCoeff_GraphOP,min_lon=0,min_lat=-90,max_lon=360,max_lat=-60,name=f'GraphOP_{considered_variable}_{week+1}week_south_pole') geographical_boxes(RCoeff_S2S,min_lon=0,min_lat=-90,max_lon=360,max_lat=-60,name=f'S2S_{considered_variable}_{week+1}week_south_pole') #geographical_boxes(RCoeff_S2S_ensemble,min_lon=0,min_lat=-90,max_lon=360,max_lat=-60,name=f'S2S_ensemble_{considered_variable}_{week+1}week_south_pole') #Indian subcontinent geographical_boxes(RCoeff_GraphFT,min_lon=66,min_lat=5,max_lon=87,max_lat=28,name=f'GraphFT_{considered_variable}_{week+1}week_india') geographical_boxes(RCoeff_GraphOP,min_lon=66,min_lat=5,max_lon=87,max_lat=28,name=f'GraphOP_{considered_variable}_{week+1}week_india') geographical_boxes(RCoeff_S2S,min_lon=66,min_lat=5,max_lon=87,max_lat=28,name=f'S2S_{considered_variable}_{week+1}week_india') #geographical_boxes(RCoeff_S2S_ensemble,min_lon=66,min_lat=5,max_lon=87,max_lat=28,name=f'S2S_ensemble_{considered_variable}_{week+1}week_india') #Tibetan plateau geographical_boxes(RCoeff_GraphFT,83,28,91,34,f'GraphFT_{considered_variable}_{week+1}week_tibetan_plateau') geographical_boxes(RCoeff_GraphOP,83,28,91,34,f'GraphOP_{considered_variable}_{week+1}week_tibetan_plateau') geographical_boxes(RCoeff_S2S,83,28,91,34,f'S2S_{considered_variable}_{week+1}week_tibetan_plateau') #geographical_boxes(RCoeff_S2S_ensemble,83,28,91,34,f'S2S_ensemble_{considered_variable}_{week+1}week_tibetan_plateau') #Northeastern Asia - Japan geographical_boxes(RCoeff_GraphFT,90,29,150,70,f'GraphFT_{considered_variable}_{week+1}week_northeastern_asia') geographical_boxes(RCoeff_GraphOP,90,29,150,70,f'GraphOP_{considered_variable}_{week+1}week_northeastern_asia') geographical_boxes(RCoeff_S2S,90,29,150,70,f'S2S_{considered_variable}_{week+1}week_northeastern_asia') #geographical_boxes(RCoeff_S2S_ensemble,90,29,150,70,f'S2S_ensemble_{considered_variable}_{week+1}week_northeastern_asia') #North - Atlantic geographical_boxes(RCoeff_GraphFT,299,14,351,50,f'GraphFT_{considered_variable}_{week+1}week_north_atlantic') geographical_boxes(RCoeff_GraphOP,299,14,351,50,f'GraphOP_{considered_variable}_{week+1}week_north_atlantic') geographical_boxes(RCoeff_S2S,299,14,351,50,f'S2S_{considered_variable}_{week+1}week_north_atlantic') #geographical_boxes(RCoeff_S2S_ensemble,299,14,351,50,f'S2S_ensemble_{considered_variable}_{week+1}week_north_atlantic') #Europe !!! coords = [longitude - 9 for longitude in RCoeff_GraphFT.longitude] geographical_boxes(RCoeff_GraphFT.assign_coords({'longitude':coords}),0,35,47,67,f'GraphFT_{considered_variable}_{week+1}week_europe') geographical_boxes(RCoeff_GraphOP.assign_coords({'longitude':coords}),0,35,47,67,f'GraphOP_{considered_variable}_{week+1}week_europe') geographical_boxes(RCoeff_S2S.assign_coords({'longitude':coords}),0,35,47,67,f'S2S_{considered_variable}_{week+1}week_europe') #geographical_boxes(RCoeff_S2S_ensemble.assign_coords({'longitude':coords}),0,35,47,67,f'S2S_ensemble_{considered_variable}_{week+1}week_europe') #South America geographical_boxes(RCoeff_GraphFT,279,-34,314,8,f'GraphFT_{considered_variable}_{week+1}week_south_america') geographical_boxes(RCoeff_GraphOP,279,-34,314,8,f'GraphOP_{considered_variable}_{week+1}week_south_america') geographical_boxes(RCoeff_S2S,279,-34,314,8,f'S2S_{considered_variable}_{week+1}week_south_america') #geographical_boxes(RCoeff_S2S_ensemble,279,-34,314,8,f'S2S_ensemble_{considered_variable}_{week+1}week_south_america') #Australia geographical_boxes(RCoeff_GraphFT,117,-38,151,-18,f'GraphFT_{considered_variable}_{week+1}week_australia') geographical_boxes(RCoeff_GraphOP,117,-38,151,-18,f'GraphOP_{considered_variable}_{week+1}week_australia') geographical_boxes(RCoeff_S2S,117,-38,151,-18,f'S2S_{considered_variable}_{week+1}week_australia') #geographical_boxes(RCoeff_S2S_ensemble,117,-38,151,-18,f'S2S_ensemble_{considered_variable}_{week+1}week_australia') #Central Asia geographical_boxes(RCoeff_GraphFT,72,35,103,50,f'GraphFT_{considered_variable}_{week+1}week_central_asia') geographical_boxes(RCoeff_GraphOP,72,35,103,50,f'GraphOP_{considered_variable}_{week+1}week_central_asia') geographical_boxes(RCoeff_S2S,72,35,103,50,f'S2S_{considered_variable}_{week+1}week_central_asia') #geographical_boxes(RCoeff_S2S_ensemble,72,35,103,50,f'S2S_ensemble_{considered_variable}_{week+1}week_central_asia') #Northern Hemisphere geographical_boxes(RCoeff_GraphFT,0,30,360,90,f'GraphFT_{considered_variable}_{week+1}week_north_hemisphere') geographical_boxes(RCoeff_GraphOP,0,30,360,90,f'GraphOP_{considered_variable}_{week+1}week_north_hemisphere') geographical_boxes(RCoeff_S2S,0,30,360,90,f'S2S_{considered_variable}_{week+1}week_north_hemisphere') #geographical_boxes(RCoeff_S2S_ensemble,0,30,360,90,f'S2S_ensemble_{considered_variable}_{week+1}week_north_hemisphere') #Equatorial band geographical_boxes(RCoeff_GraphFT,0,-30,360,30,f'GraphFT_{considered_variable}_{week+1}week_equatorial_band') geographical_boxes(RCoeff_GraphOP,0,-30,360,30,f'GraphOP_{considered_variable}_{week+1}week_equatorial_band') geographical_boxes(RCoeff_S2S,0,-30,360,30,f'S2S_{considered_variable}_{week+1}week_equatorial_band') #geographical_boxes(RCoeff_S2S_ensemble,0,-30,360,30,f'S2S_ensemble_{considered_variable}_{week+1}week_equatorial_band') #Southern Hemisphere geographical_boxes(RCoeff_GraphFT,0,-90,360,-30,f'GraphFT_{considered_variable}_{week+1}week_south_hemisphere') geographical_boxes(RCoeff_GraphOP,0,-90,360,-30,f'GraphOP_{considered_variable}_{week+1}week_south_hemisphere') geographical_boxes(RCoeff_S2S,0,-90,360,-30,f'S2S_{considered_variable}_{week+1}week_south_hemisphere') #geographical_boxes(RCoeff_S2S_ensemble,0,-90,360,-30,f'S2S_ensemble_{considered_variable}_{week+1}week_south_hemisphere') #French Polynesia geographical_boxes(RCoeff_GraphFT,180,-35,240,0,f'GraphFT_{considered_variable}_{week+1}week_french_polynesia') geographical_boxes(RCoeff_GraphOP,180,-35,240,0,f'GraphOP_{considered_variable}_{week+1}week_french_polynesia') geographical_boxes(RCoeff_S2S,180,-35,240,0,f'S2S_{considered_variable}_{week+1}week_french_polynesia') #geographical_boxes(RCoeff_S2S_ensemble,180,-35,240,0,f'S2S_ensemble_{considered_variable}_{week+1}week_french_polynesia') #box3_1 coordinates min_lon = 180 min_lat = 0 max_lon = 200 max_lat = -35 RCoeff_GraphFT_Box31 = RCoeff_GraphFT.sel(latitude=slice(min_lat,max_lat), longitude=slice(min_lon,max_lon)).mean().values RCoeff_GraphOP_Box31 = RCoeff_GraphOP.sel(latitude=slice(min_lat,max_lat), longitude=slice(min_lon,max_lon)).mean().values RCoeff_S2S_Box31 = RCoeff_S2S.sel(latitude=slice(min_lat,max_lat), longitude=slice(min_lon,max_lon)).mean().values np.save(f'/home/vsansi01/R_boxes/RCoeff_GraphFT_Box31_{considered_variable}_{week+1}.npy',RCoeff_GraphFT_Box31) np.save(f'/home/vsansi01/R_boxes/RCoeff_GraphOP_Box31_{considered_variable}_{week+1}.npy',RCoeff_GraphOP_Box31) np.save(f'/home/vsansi01/R_boxes/RCoeff_S2S_Box31.npy_{considered_variable}_{week+1}.npy',RCoeff_S2S_Box31) #box3_2 coordinates min_lon = 200 min_lat = 0 max_lon = 240 max_lat = -35 RCoeff_GraphFT_Box32 = RCoeff_GraphFT.sel(latitude=slice(min_lat,max_lat), longitude=slice(min_lon,max_lon)).mean().values RCoeff_GraphOP_Box32 = RCoeff_GraphOP.sel(latitude=slice(min_lat,max_lat), longitude=slice(min_lon,max_lon)).mean().values RCoeff_S2S_Box32 = RCoeff_S2S.sel(latitude=slice(min_lat,max_lat), longitude=slice(min_lon,max_lon)).mean().values np.save(f'/home/vsansi01/R_boxes/RCoeff_GraphFT_Box32_{considered_variable}_{week+1}.npy',RCoeff_GraphFT_Box32) np.save(f'/home/vsansi01/R_boxes/RCoeff_GraphOP_Box32_{considered_variable}_{week+1}.npy',RCoeff_GraphOP_Box32) np.save(f'/home/vsansi01/R_boxes/RCoeff_S2S_Box32_{considered_variable}_{week+1}.npy',RCoeff_S2S_Box32) print('-----> Plotting') import numpy as np import matplotlib.pyplot as plt import cartopy.crs as ccrs import cartopy.feature as cfeature from matplotlib.colors import ListedColormap, Normalize import matplotlib.ticker as mticker # Custom colormap blending blues and oranges top = matplotlib.colormaps.get_cmap('Oranges_r') bottom = matplotlib.colormaps.get_cmap('Blues') newcolors = np.vstack((bottom(np.linspace(0, 1, 128)[::-1]), top(np.linspace(0, 1, 128)[::-1]))) newcmp = ListedColormap(newcolors, name='OrangeBlue') contour_levels = np.linspace(-1, 1, 40) diff_levels = np.linspace(-0.5, 0.5, 40) norm = Normalize(vmin=-1, vmax=1) if week in [1, 2]: # Initialize figure and axes for subplots fig, axs = plt.subplots(2, 3, figsize=(20, 10), subplot_kw={'projection': ccrs.Robinson(central_longitude=180)}) # Flattening the array of axes for easy indexing axs = axs.flatten() # Plot 1: TCC for GraphFT im1 = axs[2].contourf(S2S_lon, S2S_lat, RCoeff_GraphFT, levels=contour_levels, cmap=newcmp, transform=ccrs.PlateCarree(), norm=norm, extend='both') axs[2].set_title(f'Fig. c) Rt for Week {week + 1} - GraphFT', fontsize=12, pad=10) axs[2].add_feature(cfeature.LAND, facecolor="none", edgecolor='black', linewidths=1.5, zorder=2) # Plot 2: TCC for GraphOP im2 = axs[1].contourf(S2S_lon, S2S_lat, RCoeff_GraphOP, levels=contour_levels, cmap=newcmp, transform=ccrs.PlateCarree(), norm=norm, extend='both') axs[1].set_title(f'Fig. b) Rt for Week {week + 1} - GraphOP', fontsize=12, pad=10) axs[1].add_feature(cfeature.LAND, facecolor="none", edgecolor='black', linewidths=1.5, zorder=2) # Plot 3: TCC for S2S im3 = axs[0].contourf(S2S_lon, S2S_lat, RCoeff_S2S, levels=contour_levels, cmap=newcmp, transform=ccrs.PlateCarree(), norm=norm, extend='both') axs[0].set_title(f'Fig. a) Rt for Week {week + 1} - S2S', fontsize=12, pad=10) axs[0].add_feature(cfeature.LAND, facecolor="none", edgecolor='black', linewidths=1.5, zorder=2) # Difference plots R1 = RCoeff_GraphFT.values - RCoeff_GraphOP.values R2 = RCoeff_GraphFT.values - RCoeff_S2S.values # Plot 4: Difference GraphFT - GraphOP, positioned to center in row im4 = axs[5].contourf(S2S_lon, S2S_lat, R1, levels=diff_levels, cmap=newcmp, transform=ccrs.PlateCarree(), extend='both') axs[5].set_title(f'Fig. e) Difference (GraphFT - GraphOP) for Week {week + 1}', fontsize=12, pad=10) axs[5].add_feature(cfeature.LAND, facecolor="none", edgecolor='black', linewidths=1.5, zorder=2) # Plot 5: Difference GraphFT - S2S, positioned to center in row im5 = axs[3].contourf(S2S_lon, S2S_lat, R2, levels=diff_levels, cmap=newcmp, transform=ccrs.PlateCarree(), extend='both') axs[3].set_title(f'Fig. d) Difference (GraphFT - S2S) for Week {week + 1}', fontsize=12, pad=10) axs[3].add_feature(cfeature.LAND, facecolor="none", edgecolor='black', linewidths=1.5, zorder=2) # Remove the unused middle subplot in the second row fig.delaxes(axs[4]) # Shared colorbar cbar_ax = fig.add_axes([0.92, 0.15, 0.02, 0.7]) fig.colorbar(im1, cax=cbar_ax, label='Correlation Coefficient (Rt)', orientation='vertical') # Add land features and gridlines to each axis for ax in [axs[0], axs[1], axs[2], axs[3], axs[5]]: ax.add_feature(cfeature.LAND, facecolor="none", edgecolor='black', linewidth=0.8) gl = ax.gridlines(draw_labels=True, color="gray", linestyle="--", linewidth=0.5) gl.xlocator = mticker.FixedLocator(np.arange(-180, 181, 60)) gl.ylocator = mticker.FixedLocator(np.arange(-90, 91, 30)) gl.xlabel_style = {'size': 10, 'color': 'gray'} gl.ylabel_style = {'size': 10, 'color': 'gray'} gl.top_labels = gl.right_labels = False plt.tight_layout(rect=[0, 0, 0.9, 1]) if remove_seasonality: plt.savefig(f'/home/vsansi01/results_Rt_{considered_variable}_{week + 1}WeeksAhead_rf_test.jpeg', dpi=300) # Special case for TCC 4 weeks ahead elif week == 3: fig, axs = plt.subplots(1, 3, figsize=(20, 10), subplot_kw={'projection': ccrs.Robinson(central_longitude=180)}) # Plot GraphFT for 4 weeks ahead im1 = axs[1].contourf(S2S_lon, S2S_lat, RCoeff_GraphFT, levels=contour_levels, cmap=newcmp, transform=ccrs.PlateCarree(), norm=norm, extend='both') axs[1].set_title(f'Fig. b) Rt for Week {week + 1} - GraphFT', fontsize=12) axs[1].add_feature(cfeature.LAND,facecolor="none",edgecolor='black',linewidths=1.5,zorder=2) # Plot S2S for 4 weeks ahead im2 = axs[0].contourf(S2S_lon, S2S_lat, RCoeff_S2S, levels=contour_levels, cmap=newcmp, transform=ccrs.PlateCarree(), norm=norm, extend='both') axs[0].set_title(f'Fig. a) Rt for Week {week + 1} - S2S', fontsize=12) axs[0].add_feature(cfeature.LAND,facecolor="none",edgecolor='black',linewidths=1.5,zorder=2) # Plot Difference GraphFT - S2S R2 = RCoeff_GraphFT.values - RCoeff_S2S.values im3 = axs[2].contourf(S2S_lon, S2S_lat, R2, levels=diff_levels, cmap=newcmp, transform=ccrs.PlateCarree(), extend='both') axs[2].set_title(f'Fig. c) Difference (GraphFT - S2S) for Week {week + 1}', fontsize=12) axs[2].add_feature(cfeature.LAND,facecolor="none",edgecolor='black',linewidths=1.5,zorder=2) # Shared colorbar cbar_ax = fig.add_axes([0.92, 0.15, 0.02, 0.7]) fig.colorbar(im1, cax=cbar_ax, label='Correlation Coefficient (R)', orientation='vertical') # Add colorbars for each plot for idx, ax in enumerate(axs): ax.add_feature(cfeature.LAND, facecolor="none", edgecolor='black', linewidth=0.8) gl = ax.gridlines(draw_labels=True, color="gray", linestyle="--", linewidth=0.5) gl.xlocator = mticker.FixedLocator(np.arange(-180, 181, 60)) gl.ylocator = mticker.FixedLocator(np.arange(-90, 91, 30)) gl.xlabel_style = {'size': 10, 'color': 'gray'} gl.ylabel_style = {'size': 10, 'color': 'gray'} gl.top_labels = gl.right_labels = False plt.tight_layout(rect=(0,0,0.92,1)) if remove_seasonality: plt.savefig(f'/home/vsansi01/results_Rt_{considered_variable}_{week + 1}WeeksAhead_test.jpeg', dpi=300) print('-----> Plotted')