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 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 from data_utils import Reshape_arrays,geographical_boxes 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 print('-----> Remove seasonality ? Yes or No') remove_seasonality = '{}'.format(input()) remove_seasonality = True if remove_seasonality == 'Yes' else False assert remove_seasonality in [True,False], 'Must be boolean' 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!' #South pacific box including french polynesia and Fiji islands min_lon = 180 min_lat = 0 max_lon = 240 max_lat = -35 #Importing Pangu Weather's forecasts predicted_surface_2020,original_frames_indices_2020,dates_frames_2020 = opening_data_files_pangu(months=2,year=2020,considered_variable=considered_variable_pangu) predicted_surface_2021,original_frames_indices_2021,dates_frames_2021 = opening_data_files_pangu(months=2,year=2021,considered_variable=considered_variable_pangu) predicted_surface = xr.concat([predicted_surface_2020,predicted_surface_2021],dim='time') print('-----> Pangu Forecasts downloaded') lat = predicted_surface.latitude.values lon = predicted_surface.longitude.values PanguForecasts = predicted_surface #if considered_variable_pangu in ['msl','u10','v10','t2m','tp'] if considered_variable in ['u','v']: print('-----> Choose pressure level between [925, 200] hPa') considered_level = int(input()) assert considered_level in [925,200], 'Choose between 925 or 200.' PanguForecasts = PanguForecasts.sel(pressure_level=considered_level) print("-----> Processing Pangu's forecasts for mean fields in week 2 and 3 of the forecasts") weekly_PanguForecasts = weekly_mean_operation(PanguForecasts,considered_variable_pangu) weekly_PanguForecasts = np.swapaxes(weekly_PanguForecasts,0,1) print("-----> Importing ERA-5 reanalysis") #pressure_levels_observations = xr.open_dataset('/home/hopuare/pressure_fields_pangu_weather.grib',engine='cfgrib') surface_fields_observations = xr.open_dataset('/home/vsansi01/era5_2020_2021.nc') dates = np.array(surface_fields_observations.time.values) date_start = '2020-01-01' date_stop = '2020-12-31' ####!!!!! A changer a la fin data_processed = surface_fields_observations[VariablesDict[considered_variable]]#.sel(time=slice(date_start,date_stop)) print('-----> Processing data at 0.25° resolution') original_frames = [np.expand_dims(data_processed.isel(time=slice(i,i+29)).values,axis=0) for i in range(data_processed['time'].size-30)] Dates_Observations = [data_processed.isel(time=slice(i,i+29)).time.values[0] for i in range(data_processed['time'].size-30)] original_frames = np.concatenate(original_frames,axis=0) Dates_Observations = pd.DatetimeIndex(Dates_Observations).floor('D') weekly_original_frames = weekly_mean_operation(original_frames) print('-----> Importing S2S forecasts') #real-time forecasts #if considered_variable in ['u10','v10','tp']: # path_to_S2S_data = '/home/vsansi01/S2S_forecasts/real_time_forecasts/ecmwf_forecasts_winds_precipitation.nc' #elif considered_variable in ['t2m','t2d']: # path_to_S2S_data = '/home/vsansi01/S2S_forecasts/real_time_forecasts/ecmwf_forecasts_surface_temperature.nc' #real-time forecasts ensemble mean #if considered_variable in ['u10','v10','tp']: # path_to_S2S_data = '/home/vsansi01/S2S_forecasts/real_time_forecasts/ecmwf_forecasts_winds_precipitation_ensemble_mean.nc' #elif considered_variable in ['t2m','t2d']: # path_to_S2S_data = '/home/vsansi01/S2S_forecasts/real_time_forecasts/ecmwf_forecasts_surface_temperature_ensemble_mean.nc' #S2S_nc = xr.open_dataset(path_to_S2S_data)[considered_variable]#.sel(time=slice(date_start,date_stop))#.sel(latitude=slice(min_lat,max_lat), longitude=slice(min_lon,max_lon)) #reforecasts if considered_variable in ['u10','v10','tp']: path = '/home/vsansi01/S2S_forecasts/reforecasts/deterministic/' S2S_nc = processing_hindcasts(path,False,False)[considered_variable] elif considered_variable in ['t2m']: path = '/home/vsansi01/S2S_forecasts/reforecasts/deterministic/' S2S_nc = processing_hindcasts(path,True,False)[considered_variable] print('-----> S2S imported and processed') S2S_lat = S2S_nc.latitude S2S_lon = S2S_nc.longitude S2S_forecasts = S2S_nc.values weekly_S2S_forecasts = weekly_mean_operation(S2S_forecasts) print('-----> Importing graphcast forecasts at 0.25°') os.chdir('/home/vsansi01/graphcastOP_forecasts') files = glob.glob(considered_variable+'_*',recursive=True) predictions_graphcast_operational_nc = open_data_files_graphcast(files,considered_variable,False)#.sel(time=slice(date_start,date_stop)) DatesForecastGraphOperational = pd.DatetimeIndex(predictions_graphcast_operational_nc.time.dt.date).floor('D') OperationalGraphForecasts = np.squeeze(predictions_graphcast_operational_nc.to_array().fillna(0).values) Weekly_OperationalGraphForecasts = weekly_mean_operation(OperationalGraphForecasts,graphcast=True) print('-----> Imported graphcast forecasts at 0.25°') print('-----> Importing forecasts from Fine-tuned Graphcast') #changed for now the dates so that we can have a correlation predictions_FT_2019_14_days_unified FineTunedGraphForecasts = xr.open_dataset('/home/vsansi01/predictions_FT_2019_21_days_reforecasts_unified.nc') FineTunedGraphForecasts = FineTunedGraphForecasts[VariablesDict_Graphcast[considered_variable]] Dates_ForecastsFineTunedGraphcast = np.load('/home/vsansi01/predictions_dates_4WeeksAhead_graphcast_reforecasts_unified.npy') Dates_ForecastsFineTunedGraphcast = pd.DatetimeIndex(Dates_ForecastsFineTunedGraphcast).floor('D') #manually removing seasonlality for FinetunedGraphcast if remove_seasonality: days_ahead = 21 if week in [0,1,2] else 29 climatology = calling_climatology(considered_variable) climatology = RegriddingData(climatology,FineTunedGraphForecasts) #convert dates to dayofyear as climatology is expressed as dayofyear dates_DayOfYear = [pd.to_datetime(date).timetuple().tm_yday for date in Dates_ForecastsFineTunedGraphcast] WeeklyClimato = [np.expand_dims(climatology.sel(dayofyear=slice(day,day+days_ahead)).values,axis=0) for day in dates_DayOfYear] #21 WeeklyClimato = np.concatenate(WeeklyClimato,axis=0) WeeklyClimato = weekly_mean_operation(WeeklyClimato) WeeklyClimato = WeeklyClimato[:,week,...] FineTunedGraphForecasts = FineTunedGraphForecasts - WeeklyClimato print('-----> Fine-tuned Graphcast model imported') print('-----> Regridding 1.5° for comparison with S2S forecasts') #regridding raw observations to 1.5° #data processed being ERA5 reanalysis Observations150km = RegriddingData(data_processed,S2S_nc,weekly_original_frames) #regridding Pangu weather forecasts and observations to 1.5° weekly_PanguForecasts = RegriddingData(PanguForecasts,S2S_nc,weekly_PanguForecasts) #regridding graphcast forecasts at 0.25° to 1.5° Weekly_OperationalGraphForecasts = RegriddingData(predictions_graphcast_operational_nc,S2S_nc,Weekly_OperationalGraphForecasts) #regridding fine-tuned Graphcast forecasts from 1° to 1.5° weekly_FineTunedGraphForecasts = RegriddingData(FineTunedGraphForecasts,S2S_nc) weekly_FineTunedGraphForecasts = np.squeeze(weekly_FineTunedGraphForecasts.values) print('-----> Implementation of correlation between data') #1rst method #correlation between observation and graphcast fine tuned weekly forecasts on Mondays fit_dates = Dates_Observations.isin(Dates_ForecastsFineTunedGraphcast.date) second_fit_dates = Dates_ForecastsFineTunedGraphcast.isin(Dates_Observations[fit_dates].date) #updating correlated data and correlated dates Observations150km_GraphFT = Observations150km[fit_dates] Dates_Observations_GraphFT = Dates_Observations[fit_dates] weekly_FineTunedGraphForecasts = weekly_FineTunedGraphForecasts[second_fit_dates] Dates_ForecastsFineTunedGraphcast = Dates_ForecastsFineTunedGraphcast[second_fit_dates] #correlating S2S and Graphcast FT correlated_dates = pd.DatetimeIndex(S2S_nc.time.dt.date).floor('D').isin(Dates_ForecastsFineTunedGraphcast) assert True in fit_dates and True in second_fit_dates, 'No correlation between observation and graphcast fine tuned' assert True in correlated_dates, 'No correlation between S2S and Graphcast Fine Tuned' third_fit_dates = Dates_Observations.isin(pd.DatetimeIndex(S2S_nc.time.dt.date).floor('D')[correlated_dates]) fourth_fit_dates = pd.DatetimeIndex(S2S_nc.time.dt.date).floor('D')[correlated_dates].isin(Dates_Observations[third_fit_dates]) Observations150km_S2S = Observations150km[third_fit_dates] Dates_Observations_S2S = Dates_Observations[third_fit_dates] WeeklyS2S = weekly_S2S_forecasts[correlated_dates][fourth_fit_dates] Dates_S2S = S2S_nc.time.dt.date[correlated_dates][fourth_fit_dates] #correlating S2S and GraphOP correlated_dates = DatesForecastGraphOperational.isin(Dates_Observations_S2S) #assert True in fit_dates and True in second_fit_dates, 'No correlation between observation and S2S' assert True in third_fit_dates and True in fourth_fit_dates, 'No correlation between observation and S2S' assert True in correlated_dates, 'No correlation between S2S and Graph Operational' #correlation between observation and graphcast operational fifth_fit_dates = Dates_Observations.isin(DatesForecastGraphOperational[correlated_dates]) sixth_fit_dates = DatesForecastGraphOperational[correlated_dates].isin(Dates_Observations[fifth_fit_dates]) Observations150km_GraphOP = Observations150km[fifth_fit_dates] Dates_Observations_GraphOP = Dates_Observations[fifth_fit_dates] Weekly_OperationalGraphForecasts = Weekly_OperationalGraphForecasts[correlated_dates][sixth_fit_dates] DatesForecastGraphOperational = DatesForecastGraphOperational[correlated_dates][sixth_fit_dates] #correlating GraphOP and Pangu correlated_dates = pd.DatetimeIndex(PanguForecasts.time.dt.date).floor('D').isin(DatesForecastGraphOperational) #assert True in fit_dates and True in second_fit_dates, 'No correlation between observation and graphcast operational' assert True in fifth_fit_dates and True in sixth_fit_dates, 'No correlation between observations and graphcast operational' assert True in correlated_dates, 'No correlation between Graph Operational and Pangu' seventh_fit_dates = Dates_Observations.isin(pd.DatetimeIndex(PanguForecasts.time.dt.date).floor('D')[correlated_dates]) eigth_fit_dates = pd.DatetimeIndex(PanguForecasts.time.dt.date).floor('D')[correlated_dates].isin(Dates_Observations[seventh_fit_dates]) Observations150km_Pangu = Observations150km[seventh_fit_dates] Dates_Observations_Pangu = Dates_Observations[seventh_fit_dates] weekly_PanguForecasts = np.squeeze(weekly_PanguForecasts)[correlated_dates][eigth_fit_dates] #assert True in fit_dates and True in second_fit_dates, 'No correlation between observation and pangu' assert True in seventh_fit_dates and True in eigth_fit_dates, 'No correlation between observations and Pangu-Weather' print(Observations150km_GraphFT.shape) print(weekly_FineTunedGraphForecasts.shape) print(Observations150km_S2S.shape) print(WeeklyS2S.shape) print(Observations150km_GraphOP.shape) print(Weekly_OperationalGraphForecasts.shape) print(Observations150km_Pangu.shape) print(weekly_PanguForecasts.shape) print("-----> Done processing data") np.save('/home/vsansi01/observations_predictions_fields/Dates_FineTunedGraphcast_week4.npy',Dates_ForecastsFineTunedGraphcast) np.save('/home/vsansi01/observations_predictions_fields/Dates_S2S_week4.npy',Dates_Observations_S2S) np.save('/home/vsansi01/observations_predictions_fields/Dates_GraphOP_week4.npy',Dates_Observations_GraphOP) np.save('/home/vsansi01/observations_predictions_fields/Dates_Pangu_week4.npy',Dates_Observations_Pangu) stop coarse_shape = (-1,Observations150km_GraphOP.shape[2],Observations150km_GraphOP.shape[3]) #Fine-tuned graphcast observations first week chosen, no need to choose for GraphFT as the predictions is directly weekly mean forecasts. Observations150km_GraphFT = ChoosingWeek(Observations150km_GraphFT, week, coarse_shape) #S2S forecasts against ERA5 reanalysis at 1.5° Observations150km_S2S = ChoosingWeek(Observations150km_S2S, week, coarse_shape) WeeklyS2S = ChoosingWeek(WeeklyS2S, week, coarse_shape) ######## ### for 4th weekly forecasts which is not present in pangu and graphcast op ######## week_1= week-1 if week == 3 else week #Operational Pangu weather forecast against reanalysis at 1.5° Observations150km_Pangu = ChoosingWeek(Observations150km_Pangu, week_1, coarse_shape) weekly_PanguForecasts = ChoosingWeek(weekly_PanguForecasts, week_1, coarse_shape) #Graphcast OP against reanalysis at 1.5° Observations150km_GraphOP = ChoosingWeek(Observations150km_GraphOP, week_1, coarse_shape) Weekly_OperationalGraphForecasts = ChoosingWeek(Weekly_OperationalGraphForecasts, week_1, coarse_shape) if remove_seasonality: print('-----> Removing climatology') weekly_PanguForecasts = removing_climatology(considered_variable,weekly_PanguForecasts,\ dates = Dates_Observations_Pangu,week=week_1,template=S2S_nc) Observations150km_Pangu = removing_climatology(considered_variable,Observations150km_Pangu,\ dates = Dates_Observations_Pangu,week=week_1,template=S2S_nc) assert len(weekly_PanguForecasts) == len(Dates_Observations_Pangu) assert len(Observations150km_Pangu) == len(Dates_Observations_Pangu) Weekly_OperationalGraphForecasts = removing_climatology(considered_variable,Weekly_OperationalGraphForecasts,\ dates = DatesForecastGraphOperational,week=week_1,template=S2S_nc) Observations150km_GraphOP = removing_climatology(considered_variable,Observations150km_GraphOP,\ dates = DatesForecastGraphOperational,week=week_1,template=S2S_nc) assert len(Weekly_OperationalGraphForecasts) == len(DatesForecastGraphOperational) assert len(Observations150km_GraphOP) == len(DatesForecastGraphOperational) #weekly_FineTunedGraphForecasts = removing_climatology(considered_variable,weekly_FineTunedGraphForecasts, #dates = Dates_ForecastsFineTunedGraphcast,week=week,template=S2S_nc) Observations150km_GraphFT = removing_climatology(considered_variable,Observations150km_GraphFT,\ dates = Dates_ForecastsFineTunedGraphcast,week=week,template=S2S_nc) #assert len(weekly_FineTunedGraphForecasts) == len(Dates_ForecastsFineTunedGraphcast) assert len(Observations150km_GraphFT) == len(Dates_ForecastsFineTunedGraphcast) Observations150km_S2S = removing_climatology(considered_variable,Observations150km_S2S,dates = Dates_S2S,week=week,template=S2S_nc) WeeklyS2S = removing_climatology(considered_variable,WeeklyS2S,dates = Dates_S2S,week=week,template=S2S_nc) assert len(WeeklyS2S) == len(Dates_S2S) assert len(Observations150km_S2S) == len(Dates_S2S) #Change S2S_lat and S2S_lon CoordsPangu = { 'time' : Dates_Observations_Pangu, 'latitude' : S2S_lat, 'longitude' : S2S_lon } CoordsS2S = { 'time' : Dates_Observations_S2S, 'latitude' : S2S_lat, 'longitude' : S2S_lon } CoordsGraphOP = { 'time' : Dates_Observations_GraphOP, 'latitude' : S2S_lat, 'longitude' : S2S_lon } CoordsGraphFT = { 'time' : Dates_Observations_GraphFT, 'latitude' : S2S_lat, 'longitude' : S2S_lon } #Graphcast Fine tuned Observations150km_GraphFT = xr.DataArray(Observations150km_GraphFT,coords=CoordsGraphFT) weekly_FineTunedGraphForecasts = xr.DataArray(weekly_FineTunedGraphForecasts,coords=CoordsGraphFT) #Observations150km_GraphFT.to_netcdf(f'/home/vsansi01/observations_predictions_fields/Observations_GraphFT_{considered_variable}_{week+1}.nc') #weekly_FineTunedGraphForecasts.to_netcdf(f'/home/vsansi01/observations_predictions_fields/Predictions_GraphFT_{considered_variable}_{week+1}.nc') #S2S forecasts Observations150km_S2S = xr.DataArray(Observations150km_S2S,coords=CoordsS2S) WeeklyS2S = xr.DataArray(WeeklyS2S,coords=CoordsS2S) #Observations150km_S2S.to_netcdf(f'/home/vsansi01/observations_predictions_fields/Observations_S2S_ensemble_{considered_variable}_{week+1}.nc') #WeeklyS2S.to_netcdf(f'/home/vsansi01/observations_predictions_fields/Predictions_S2S_ensemble_{considered_variable}_{week+1}.nc') #Graphcast OP Observations150km_GraphOP = xr.DataArray(Observations150km_GraphOP,coords=CoordsGraphOP) Weekly_OperationalGraphForecasts = xr.DataArray(Weekly_OperationalGraphForecasts,coords=CoordsGraphOP) #Observations150km_GraphOP.to_netcdf(f'/home/vsansi01/observations_predictions_fields/Observations_GraphOP_{considered_variable}_{week+1}.nc') #Weekly_OperationalGraphForecasts.to_netcdf(f'/home/vsansi01/observations_predictions_fields/Predictions_GraphOP_{considered_variable}_{week+1}.nc') #Pangu weather forecasts Observations150km_Pangu = xr.DataArray(Observations150km_Pangu,coords=CoordsPangu) weekly_PanguForecasts = xr.DataArray(weekly_PanguForecasts,coords=CoordsPangu) #Observations150km_Pangu.to_netcdf(f'/home/vsansi01/observations_predictions_fields/Observations_Pangu_{considered_variable}_{week+1}.nc') #weekly_PanguForecasts.to_netcdf(f'/home/vsansi01/observations_predictions_fields/Predictions_Pangu_{considered_variable}_{week+1}.nc') #computing TCC spatial maps RCoeff_GraphFT = xr.corr(Observations150km_GraphFT,weekly_FineTunedGraphForecasts,dim='time') RCoeff_GraphOP = xr.corr(Observations150km_GraphOP,Weekly_OperationalGraphForecasts,dim='time') RCoeff_Pangu = xr.corr(Observations150km_Pangu,weekly_PanguForecasts,dim='time') RCoeff_S2S = xr.corr(Observations150km_S2S,WeeklyS2S,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}.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 = rmse(Observations150km_GraphFT,weekly_FineTunedGraphForecasts,dim='time') RMSE_GraphOP = rmse(Observations150km_GraphOP,Weekly_OperationalGraphForecasts,dim='time') RMSE_Pangu = rmse(Observations150km_Pangu,weekly_PanguForecasts,dim='time') RMSE_S2S = 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,RMSE_GraphFT.mean('longitude'),label='GraphFT') plt.plot(S2S_lat,RMSE_GraphOP.mean('longitude'),label='GraphOP') plt.plot(S2S_lat,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,RMSE_GraphFT_OP,label='GraphFT - GraphOP') plt.plot(S2S_lat,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_2021') ACC_GraphOP,std_GraphOP = anomaly_correlation_coefficient(Observations150km_GraphOP,Weekly_OperationalGraphForecasts, Dates_Observations_GraphOP,S2S_nc,considered_variable,week,'GraphOP_2021') #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_2021') 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) stop #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) stop if week in [0,1]: #Hawaii 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_GraphFT,min_lon=195,min_lat=14,max_lon=212,max_lat=25,name=f'S2S_{considered_variable}_{week+1}week_hawaii') #USA geographical_boxes(RCoeff_GraphFT,min_lon=232,min_lat=23,max_lon=281,max_lat=48,name=f'GraphFT_{considered_variable}_{week+1}week_usa') geographical_boxes(RCoeff_GraphOP,min_lon=232,min_lat=23,max_lon=281,max_lat=48,name=f'GraphOP_{considered_variable}_{week+1}week_usa') geographical_boxes(RCoeff_GraphFT,min_lon=232,min_lat=23,max_lon=281,max_lat=48,name=f'S2S_{considered_variable}_{week+1}week_usa') #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_GraphFT,min_lon=0,min_lat=60,max_lon=360,max_lat=90,name=f'S2S_{considered_variable}_{week+1}week_north_pole') #Northern India geographical_boxes(RCoeff_GraphFT,min_lon=66,min_lat=19,max_lon=87,max_lat=28,name=f'GraphFT_{considered_variable}_{week+1}week_north_india') geographical_boxes(RCoeff_GraphOP,min_lon=66,min_lat=19,max_lon=87,max_lat=28,name=f'GraphOP_{considered_variable}_{week+1}week_north_india') geographical_boxes(RCoeff_GraphFT,min_lon=66,min_lat=19,max_lon=87,max_lat=28,name=f'S2S_{considered_variable}_{week+1}week_north_india') #Northern India geographical_boxes(RCoeff_GraphFT,min_lon=66,min_lat=19,max_lon=87,max_lat=28,name=f'GraphFT_{considered_variable}_{week+1}week_north_india') geographical_boxes(RCoeff_GraphOP,min_lon=66,min_lat=19,max_lon=87,max_lat=28,name=f'GraphOP_{considered_variable}_{week+1}week_north_india') geographical_boxes(RCoeff_GraphFT,min_lon=66,min_lat=19,max_lon=87,max_lat=28,name=f'S2S_{considered_variable}_{week+1}week_north_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_GraphFT,83,28,91,34,f'S2S_{considered_variable}_{week+1}week_tibetan_plateau') #Tibetan plateau geographical_boxes(RCoeff_GraphFT,89,29,139,48,f'GraphFT_{considered_variable}_{week+1}week_northeastern_asia') geographical_boxes(RCoeff_GraphOP,89,29,139,48,f'GraphOP_{considered_variable}_{week+1}week_northeastern_asia') geographical_boxes(RCoeff_GraphFT,89,29,139,48,f'S2S_{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_northeastern_asia') geographical_boxes(RCoeff_GraphFT,299,14,351,50,f'S2S_{considered_variable}_{week+1}week_northeastern_asia') #Europe !!! #coords = [longitude if longitude <= 180 else longitude - 360 for longitude in RCoeff_GraphFT.longitude] #geographical_boxes(RCoeff_GraphFT.assign_coords({'longitude':coords}),-9,35,38,67,f'GraphFT_{considered_variable}_{week+1}week_europe') #geographical_boxes(RCoeff_GraphOP.assign_coords({'longitude':coords}),-9,35,38,67,f'GraphOP_{considered_variable}_{week+1}week_europe') #geographical_boxes(RCoeff_GraphFT.assign_coords({'longitude':coords}),-9,35,38,67,f'S2S_{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_GraphFT,279,-34,314,8,f'S2S_{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_GraphFT,117,-38,151,-18,f'S2S_{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_GraphFT,72,35,103,50,f'S2S_{considered_variable}_{week+1}week_central_asia') #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('/home/vsansi01/R_boxes/RCoeff_GraphFT_Box21.npy',RCoeff_GraphFT_Box21) np.save('/home/vsansi01/R_boxes/RCoeff_GraphOP_Box21.npy',RCoeff_GraphOP_Box21) np.save('/home/vsansi01/R_boxes/RCoeff_S2S_Box21.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('/home/vsansi01/R_boxes/RCoeff_GraphFT_Box22.npy',RCoeff_GraphFT_Box22) np.save('/home/vsansi01/R_boxes/RCoeff_GraphOP_Box22.npy',RCoeff_GraphOP_Box22) np.save('/home/vsansi01/R_boxes/RCoeff_S2S_Box22.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('/home/vsansi01/R_boxes/RCoeff_GraphFT_Box23.npy',RCoeff_GraphFT_Box23) np.save('/home/vsansi01/R_boxes/RCoeff_GraphOP_Box23.npy',RCoeff_GraphOP_Box23) np.save('/home/vsansi01/R_boxes/RCoeff_S2S_Box23.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('/home/vsansi01/R_boxes/RCoeff_GraphFT_Box24.npy',RCoeff_GraphFT_Box24) np.save('/home/vsansi01/R_boxes/RCoeff_GraphOP_Box24.npy',RCoeff_GraphOP_Box24) np.save('/home/vsansi01/R_boxes/RCoeff_S2S_Box24.npy',RCoeff_S2S_Box24) if week in [2,3]: #Hawaii 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_GraphFT,min_lon=195,min_lat=14,max_lon=212,max_lat=25,name=f'S2S_{considered_variable}_{week+1}week_hawaii') #USA geographical_boxes(RCoeff_GraphFT,min_lon=232,min_lat=23,max_lon=281,max_lat=48,name=f'GraphFT_{considered_variable}_{week+1}week_usa') geographical_boxes(RCoeff_GraphOP,min_lon=232,min_lat=23,max_lon=281,max_lat=48,name=f'GraphOP_{considered_variable}_{week+1}week_usa') geographical_boxes(RCoeff_GraphFT,min_lon=232,min_lat=23,max_lon=281,max_lat=48,name=f'S2S_{considered_variable}_{week+1}week_usa') #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_GraphFT,min_lon=0,min_lat=60,max_lon=360,max_lat=90,name=f'S2S_{considered_variable}_{week+1}week_north_pole') #Northern India geographical_boxes(RCoeff_GraphFT,min_lon=66,min_lat=19,max_lon=87,max_lat=28,name=f'GraphFT_{considered_variable}_{week+1}week_north_india') geographical_boxes(RCoeff_GraphOP,min_lon=66,min_lat=19,max_lon=87,max_lat=28,name=f'GraphOP_{considered_variable}_{week+1}week_north_india') geographical_boxes(RCoeff_GraphFT,min_lon=66,min_lat=19,max_lon=87,max_lat=28,name=f'S2S_{considered_variable}_{week+1}week_north_india') #Northern India geographical_boxes(RCoeff_GraphFT,min_lon=66,min_lat=19,max_lon=87,max_lat=28,name=f'GraphFT_{considered_variable}_{week+1}week_north_india') geographical_boxes(RCoeff_GraphOP,min_lon=66,min_lat=19,max_lon=87,max_lat=28,name=f'GraphOP_{considered_variable}_{week+1}week_north_india') geographical_boxes(RCoeff_GraphFT,min_lon=66,min_lat=19,max_lon=87,max_lat=28,name=f'S2S_{considered_variable}_{week+1}week_north_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_GraphFT,83,28,91,34,f'S2S_{considered_variable}_{week+1}week_tibetan_plateau') #Tibetan plateau geographical_boxes(RCoeff_GraphFT,89,29,139,48,f'GraphFT_{considered_variable}_{week+1}week_northeastern_asia') geographical_boxes(RCoeff_GraphOP,89,29,139,48,f'GraphOP_{considered_variable}_{week+1}week_northeastern_asia') geographical_boxes(RCoeff_GraphFT,89,29,139,48,f'S2S_{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_northeastern_asia') geographical_boxes(RCoeff_GraphFT,299,14,351,50,f'S2S_{considered_variable}_{week+1}week_northeastern_asia') #Europe !!! #coords = [longitude if longitude <= 180 else longitude - 360 for longitude in RCoeff_GraphFT.longitude] #print(RCoeff_GraphFT.assign_coords({'longitude':coords}).longitude.values) #geographical_boxes(RCoeff_GraphFT.assign_coords({'longitude':coords}),-9,35,38,67,f'GraphFT_{considered_variable}_{week+1}week_europe') #geographical_boxes(RCoeff_GraphOP.assign_coords({'longitude':coords}),-9,35,38,67,f'GraphOP_{considered_variable}_{week+1}week_europe') #geographical_boxes(RCoeff_GraphFT.assign_coords({'longitude':coords}),-9,35,38,67,f'S2S_{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_GraphFT,279,-34,314,8,f'S2S_{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_GraphFT,117,-38,151,-18,f'S2S_{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_GraphFT,72,35,103,50,f'S2S_{considered_variable}_{week+1}week_central_asia') #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('/home/vsansi01/R_boxes/RCoeff_GraphFT_Box31.npy',RCoeff_GraphFT_Box31) np.save('/home/vsansi01/R_boxes/RCoeff_GraphOP_Box31.npy',RCoeff_GraphOP_Box31) np.save('/home/vsansi01/R_boxes/RCoeff_S2S_Box31.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('/home/vsansi01/R_boxes/RCoeff_GraphFT_Box32.npy',RCoeff_GraphFT_Box32) np.save('/home/vsansi01/R_boxes/RCoeff_GraphOP_Box32.npy',RCoeff_GraphOP_Box32) np.save('/home/vsansi01/R_boxes/RCoeff_S2S_Box32.npy',RCoeff_S2S_Box32) print('-----> Plotting') 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') fig = plt.figure(figsize=(20,10)) ax1 = fig.add_subplot(2,9,(7,9), projection=ccrs.Robinson(central_longitude=180)) im1 = ax1.contourf(S2S_lon,S2S_lat,RCoeff_GraphFT,np.linspace(-1,1,40),transform=ccrs.PlateCarree(),vmin=-1,vmax=1,cmap=newcmp,extend='both') #if week in [0,1]: # ax1.add_patch(mpatches.Rectangle(xy=[-180, -35], width=20, height=35,facecolor='none', edgecolor='k',transform=ccrs.PlateCarree(),linewidth=3)) # ax1.add_patch(mpatches.Rectangle(xy=[-160, -12], width=40, height=12,facecolor='none', edgecolor='k',transform=ccrs.PlateCarree(),linewidth=3)) # ax1.add_patch(mpatches.Rectangle(xy=[-160, -20], width=40, height=8,facecolor='none', edgecolor='k',transform=ccrs.PlateCarree(),linewidth=3)) # ax1.add_patch(mpatches.Rectangle(xy=[-160, -35], width=40, height=15,facecolor='none', edgecolor='k',transform=ccrs.PlateCarree(),linewidth=3)) # ax1.text(-180,-30,r'$\overline{R1}$ = ' + '{:.2f}'.format(RCoeff_GraphFT_Box21),transform=ccrs.PlateCarree()) # ax1.text(-135,-5,r'$\overline{R2}$ = ' + '{:.2f}'.format(RCoeff_GraphFT_Box22),transform=ccrs.PlateCarree()) # ax1.text(-135,-17.5,r'$\overline{R3}$ = ' + '{:.2f}'.format(RCoeff_GraphFT_Box23),transform=ccrs.PlateCarree()) # ax1.text(-135,-30,r'$\overline{R4}$ = ' + '{:.2f}'.format(RCoeff_GraphFT_Box24),transform=ccrs.PlateCarree()) #elif week == 2: # ax1.add_patch(mpatches.Rectangle(xy=[-180, -35], width=20, height=35,facecolor='none', edgecolor='k',transform=ccrs.PlateCarree(),linewidth=3)) # ax1.add_patch(mpatches.Rectangle(xy=[-160, -35], width=40, height=35,facecolor='none', edgecolor='k',transform=ccrs.PlateCarree(),linewidth=3)) # ax1.text(-175,-10,r'$\overline{R1}$ = ' + '{:.2f}'.format(RCoeff_GraphFT_Box31),transform=ccrs.PlateCarree()) # ax1.text(-135,-30,r'$\overline{R2}$ = ' + '{:.2f}'.format(RCoeff_GraphFT_Box32),transform=ccrs.PlateCarree()) ax1.xformatter = LONGITUDE_FORMATTER ax1.yformatter = LATITUDE_FORMATTER ax1.gridlines(draw_labels=["x", "y", "bottom", "left", "geo"]) ax1.add_feature(cfeature.LAND,facecolor="none",edgecolor='black',linewidths=1.5,zorder=2) #ax1.set_extent([-180, -120, -35, 0],crs=ccrs.PlateCarree()) plt.title('TCC for W {} for GraphFT'.format(week+1)) ax2 = fig.add_subplot(2,9,(4,6), projection=ccrs.Robinson(central_longitude=180)) im2 = ax2.contourf(S2S_lon,S2S_lat,RCoeff_GraphOP,np.linspace(-1,1,40),transform=ccrs.PlateCarree(),vmin=-1,vmax=1,cmap=newcmp,extend='both') #if week in [0,1]: # ax2.add_patch(mpatches.Rectangle(xy=[-180, -35], width=20, height=35,facecolor='none', edgecolor='k',transform=ccrs.PlateCarree(),linewidth=3)) # ax2.add_patch(mpatches.Rectangle(xy=[-160, -12], width=40, height=12,facecolor='none', edgecolor='k',transform=ccrs.PlateCarree(),linewidth=3)) # ax2.add_patch(mpatches.Rectangle(xy=[-160, -20], width=40, height=8,facecolor='none', edgecolor='k',transform=ccrs.PlateCarree(),linewidth=3)) # ax2.add_patch(mpatches.Rectangle(xy=[-160, -35], width=40, height=15,facecolor='none', edgecolor='k',transform=ccrs.PlateCarree(),linewidth=3)) # ax2.text(-180,-30,r'$\overline{R1}$ = ' + '{:.2f}'.format(RCoeff_GraphOP_Box21),transform=ccrs.PlateCarree()) # ax2.text(-135,-5,r'$\overline{R2}$ = ' + '{:.2f}'.format(RCoeff_GraphOP_Box22),transform=ccrs.PlateCarree()) # ax2.text(-135,-17.5,r'$\overline{R3}$ = ' + '{:.2f}'.format(RCoeff_GraphOP_Box23),transform=ccrs.PlateCarree()) # ax2.text(-135,-30,r'$\overline{R4}$ = ' + '{:.2f}'.format(RCoeff_GraphOP_Box24),transform=ccrs.PlateCarree()) #elif week == 2: # ax2.add_patch(mpatches.Rectangle(xy=[-180, -35], width=20, height=35,facecolor='none', edgecolor='k',transform=ccrs.PlateCarree(),linewidth=3)) # ax2.add_patch(mpatches.Rectangle(xy=[-160, -35], width=40, height=35,facecolor='none', edgecolor='k',transform=ccrs.PlateCarree(),linewidth=3)) # ax2.text(-175,-10,r'$\overline{R1}$ = ' + '{:.2f}'.format(RCoeff_GraphOP_Box31),transform=ccrs.PlateCarree()) # ax2.text(-135,-30,r'$\overline{R2}$ = ' + '{:.2f}'.format(RCoeff_GraphOP_Box32),transform=ccrs.PlateCarree()) ax2.xformatter = LONGITUDE_FORMATTER ax2.yformatter = LATITUDE_FORMATTER ax2.gridlines(draw_labels=["x", "y", "bottom", "left", "geo"]) ax2.add_feature(cfeature.LAND,facecolor="none",edgecolor='black',linewidths=1.5,zorder=2) #ax2.set_extent([-180, -120, -35, 0],crs=ccrs.PlateCarree()) plt.title('Rt for W {} for GraphOP'.format(week+1)) #ax3 = fig.add_subplot(223, projection=ccrs.Robinson(central_longitude=180)) #im3 = ax3.contourf(S2S_lon,S2S_lat,RCoeff_Pangu,np.linspace(-1,1,40),transform=ccrs.PlateCarree(),vmin=-1,vmax=1,cmap=newcmp,extend='both') #if week in [0,1]: # ax3.add_patch(mpatches.Rectangle(xy=[-180, -35], width=20, height=35,facecolor='none', edgecolor='k',transform=ccrs.PlateCarree(),linewidth=3)) # ax3.add_patch(mpatches.Rectangle(xy=[-160, -12], width=40, height=12,facecolor='none', edgecolor='k',transform=ccrs.PlateCarree(),linewidth=3)) # ax3.add_patch(mpatches.Rectangle(xy=[-160, -20], width=40, height=8,facecolor='none', edgecolor='k',transform=ccrs.PlateCarree(),linewidth=3)) # ax3.add_patch(mpatches.Rectangle(xy=[-160, -35], width=40, height=15,facecolor='none', edgecolor='k',transform=ccrs.PlateCarree(),linewidth=3)) # ax3.text(-180,-30,r'$\overline{R1}$ = ' + '{:.2f}'.format(RCoeff_Pangu_Box21),transform=ccrs.PlateCarree()) # ax3.text(-135,-5,r'$\overline{R2}$ = ' + '{:.2f}'.format(RCoeff_Pangu_Box22),transform=ccrs.PlateCarree()) # ax3.text(-135,-17.5,r'$\overline{R3}$ = ' + '{:.2f}'.format(RCoeff_Pangu_Box23),transform=ccrs.PlateCarree()) # ax3.text(-135,-30,r'$\overline{R4}$ = ' + '{:.2f}'.format(RCoeff_Pangu_Box24),transform=ccrs.PlateCarree()) #elif week == 2: # ax3.add_patch(mpatches.Rectangle(xy=[-180, -35], width=20, height=35,facecolor='none', edgecolor='k',transform=ccrs.PlateCarree(),linewidth=3)) # ax3.add_patch(mpatches.Rectangle(xy=[-160, -35], width=40, height=35,facecolor='none', edgecolor='k',transform=ccrs.PlateCarree(),linewidth=3)) # ax3.text(-175,-10,r'$\overline{R1}$ = ' + '{:.2f}'.format(RCoeff_Pangu_Box31),transform=ccrs.PlateCarree()) # ax3.text(-135,-30,r'$\overline{R2}$ = ' + '{:.2f}'.format(RCoeff_Pangu_Box32),transform=ccrs.PlateCarree()) #ax3.xformatter = LONGITUDE_FORMATTER #ax3.yformatter = LATITUDE_FORMATTER #ax3.gridlines(draw_labels=["x", "y", "bottom", "left", "geo"]) #ax3.add_feature(cfeature.LAND,facecolor="none",edgecolor='black',linewidths=1.5,zorder=2) #ax3.set_extent([-180, -120, -35, 0],crs=ccrs.PlateCarree()) #plt.title('R correlation for W+{} Pangu'.format(week+1)) ax4 = fig.add_subplot(2,9,(1,3), projection=ccrs.Robinson(central_longitude=180)) im4 = ax4.contourf(S2S_lon,S2S_lat,RCoeff_S2S,np.linspace(-1,1,40),transform=ccrs.PlateCarree(),vmin=-1,vmax=1,cmap=newcmp,extend='both') #if week in [0,1]: # ax4.add_patch(mpatches.Rectangle(xy=[-180, -35], width=20, height=35,facecolor='none', edgecolor='k',transform=ccrs.PlateCarree(),linewidth=3)) # ax4.add_patch(mpatches.Rectangle(xy=[-160, -12], width=40, height=12,facecolor='none', edgecolor='k',transform=ccrs.PlateCarree(),linewidth=3)) # ax4.add_patch(mpatches.Rectangle(xy=[-160, -20], width=40, height=8,facecolor='none', edgecolor='k',transform=ccrs.PlateCarree(),linewidth=3)) # ax4.add_patch(mpatches.Rectangle(xy=[-160, -35], width=40, height=15,facecolor='none', edgecolor='k',transform=ccrs.PlateCarree(),linewidth=3)) # ax4.text(-180,-30,r'$\overline{R1}$ = ' + '{:.2f}'.format(RCoeff_S2S_Box21),transform=ccrs.PlateCarree()) # ax4.text(-135,-5,r'$\overline{R2}$ = ' + '{:.2f}'.format(RCoeff_S2S_Box22),transform=ccrs.PlateCarree()) # ax4.text(-135,-17.5,r'$\overline{R3}$ = ' + '{:.2f}'.format(RCoeff_S2S_Box23),transform=ccrs.PlateCarree()) # ax4.text(-135,-30,r'$\overline{R4}$ = ' + '{:.2f}'.format(RCoeff_S2S_Box24),transform=ccrs.PlateCarree()) #elif week == 2: # ax4.add_patch(mpatches.Rectangle(xy=[-180, -35], width=20, height=35,facecolor='none', edgecolor='k',transform=ccrs.PlateCarree(),linewidth=3)) # ax4.add_patch(mpatches.Rectangle(xy=[-160, -35], width=40, height=35,facecolor='none', edgecolor='k',transform=ccrs.PlateCarree(),linewidth=3)) # ax4.text(-175,-10,r'$\overline{R1}$ = ' + '{:.2f}'.format(RCoeff_S2S_Box31),transform=ccrs.PlateCarree()) # ax4.text(-135,-30,r'$\overline{R2}$ = ' + '{:.2f}'.format(RCoeff_S2S_Box32),transform=ccrs.PlateCarree()) ax4.xformatter = LONGITUDE_FORMATTER ax4.yformatter = LATITUDE_FORMATTER ax4.gridlines(draw_labels=["x", "y", "bottom", "left", "geo"]) ax4.add_feature(cfeature.LAND,facecolor="none",edgecolor='black',linewidths=1.5,zorder=2) #ax4.set_extent([-180, -120, -35, 0],crs=ccrs.PlateCarree()) plt.title('Rt for week {} for S2S'.format(week+1)) #plt.savefig('/home/vsansi01/results_{}_{}WeeksAhead_cf.jpeg'.format(considered_variable,week+1),dpi=300) #fig1 = plt.figure(figsize=(10,5)) R1 = RCoeff_GraphFT.values - RCoeff_GraphOP.values ax5 = fig.add_subplot(2,9,(15,17), projection=ccrs.Robinson(central_longitude=180)) im5 = ax5.contourf(S2S_lon,S2S_lat,R1,np.linspace(-0.5,0.5,40),transform=ccrs.PlateCarree(),vmin=-0.5,vmax=0.5,cmap=newcmp,extend='both') ax5.xformatter = LONGITUDE_FORMATTER ax5.yformatter = LATITUDE_FORMATTER ax5.gridlines(draw_labels=["x", "y", "bottom", "left", "geo"]) ax5.add_feature(cfeature.LAND,facecolor="none",edgecolor='black',linewidths=1.5,zorder=2) #ax1.set_extent([-180, -120, -35, 0],crs=ccrs.PlateCarree()) plt.title('Rt difference for W {} for GraphFT - GraphOP'.format(week+1),y=1) R2 = RCoeff_GraphFT.values - RCoeff_S2S.values ax6 = fig.add_subplot(2,9,(11,13), projection=ccrs.Robinson(central_longitude=180)) im6 = ax6.contourf(S2S_lon,S2S_lat,R2,np.linspace(-0.5,0.5,40),transform=ccrs.PlateCarree(),vmin=-0.5,vmax=0.5,cmap=newcmp,extend='both') ax6.xformatter = LONGITUDE_FORMATTER ax6.yformatter = LATITUDE_FORMATTER ax6.gridlines(draw_labels=["x", "y", "bottom", "left", "geo"]) ax6.add_feature(cfeature.LAND,facecolor="none",edgecolor='black',linewidths=1.5,zorder=2) #ax2.set_extent([-180, -120, -35, 0],crs=ccrs.PlateCarree()) plt.title('Rt difference of W {} for GraphFT - S2S'.format(week+1),y=1) c1 = fig.colorbar(im1,ax=ax1,shrink=0.6) c2 = fig.colorbar(im2,ax=ax2,shrink=0.6) #c3 = fig.colorbar(im3,ax=ax3,shrink=0.6) c4 = fig.colorbar(im4,ax=ax4,shrink=0.6) c5 = fig.colorbar(im5,ax=ax5,shrink=0.6) c6 = fig.colorbar(im6,ax=ax6,shrink=0.6) if remove_seasonality: plt.savefig('/home/vsansi01/results_Rt_{}_{}WeeksAhead_rf.jpeg'.format(considered_variable,week+1),dpi=300) #TCC 4 weeks ahead if week == 3: fig = plt.figure(figsize=(20,10)) ax1 = fig.add_subplot(1,9,(4,6), projection=ccrs.Robinson(central_longitude=180)) im1 = ax1.contourf(S2S_lon,S2S_lat,RCoeff_GraphFT,np.linspace(-1,1,40),transform=ccrs.PlateCarree(),vmin=-1,vmax=1,cmap=newcmp,extend='both') ax1.xformatter = LONGITUDE_FORMATTER ax1.yformatter = LATITUDE_FORMATTER ax1.gridlines(draw_labels=["x", "y", "bottom", "left", "geo"]) ax1.add_feature(cfeature.LAND,facecolor="none",edgecolor='black',linewidths=1.5,zorder=2) #ax1.set_extent([-180, -120, -35, 0],crs=ccrs.PlateCarree()) plt.title('Rt for W {} for GraphFT'.format(week+1)) ax2 = fig.add_subplot(1,9,(1,3), projection=ccrs.Robinson(central_longitude=180)) im2 = ax2.contourf(S2S_lon,S2S_lat,RCoeff_S2S,np.linspace(-1,1,40),transform=ccrs.PlateCarree(),vmin=-1,vmax=1,cmap=newcmp,extend='both') ax2.xformatter = LONGITUDE_FORMATTER ax2.yformatter = LATITUDE_FORMATTER ax2.gridlines(draw_labels=["x", "y", "bottom", "left", "geo"]) ax2.add_feature(cfeature.LAND,facecolor="none",edgecolor='black',linewidths=1.5,zorder=2) #ax2.set_extent([-180, -120, -35, 0],crs=ccrs.PlateCarree()) plt.title('Rt for W {} for S2S'.format(week+1)) ax5 = fig.add_subplot(1,9,(7,9), projection=ccrs.Robinson(central_longitude=180)) im5 = ax5.contourf(S2S_lon,S2S_lat,R2,np.linspace(-0.5,0.5,40),transform=ccrs.PlateCarree(),vmin=-0.5,vmax=0.5,cmap=newcmp,extend='both') ax5.xformatter = LONGITUDE_FORMATTER ax5.yformatter = LATITUDE_FORMATTER ax5.gridlines(draw_labels=["x", "y", "bottom", "left", "geo"]) ax5.add_feature(cfeature.LAND,facecolor="none",edgecolor='black',linewidths=1.5,zorder=2) plt.title('Rt difference for W {} for GraphFT - S2S'.format(week+1),y=1) c1 = fig.colorbar(im1,ax=ax1,shrink=0.6) c2 = fig.colorbar(im2,ax=ax2,shrink=0.6) #c3 = fig.colorbar(im3,ax=ax3,shrink=0.6) c5 = fig.colorbar(im5,ax=ax5,shrink=0.6) plt.tight_layout() if remove_seasonality: plt.savefig('/home/vsansi01/results_Rt_{}_{}WeeksAhead_rf.jpeg'.format(considered_variable,week+1),dpi=300) #RMSE newcolors = top(np.linspace(0, 1, 128)[::-1]) newcmp = ListedColormap(newcolors, name='OrangeBlue') fig = plt.figure(figsize=(20,10)) ax1 = fig.add_subplot(2,9,(7,9), projection=ccrs.Robinson(central_longitude=180)) im1 = ax1.contourf(S2S_lon,S2S_lat,RMSE_GraphFT,np.linspace(0,7,40),transform=ccrs.PlateCarree(),vmin=0,vmax=7,cmap=newcmp,extend='both') ax1.xformatter = LONGITUDE_FORMATTER ax1.yformatter = LATITUDE_FORMATTER ax1.gridlines(draw_labels=["x", "y", "bottom", "left", "geo"]) ax1.add_feature(cfeature.LAND,facecolor="none",edgecolor='black',linewidths=1.5,zorder=2) #ax1.set_extent([-180, -120, -35, 0],crs=ccrs.PlateCarree()) plt.title('RMSE for W {} for GraphFT'.format(week+1)) ax2 = fig.add_subplot(2,9,(4,6), projection=ccrs.Robinson(central_longitude=180)) im2 = ax2.contourf(S2S_lon,S2S_lat,RMSE_GraphOP,np.linspace(0,7,40),transform=ccrs.PlateCarree(),vmin=0,vmax=7,cmap=newcmp,extend='both') ax2.xformatter = LONGITUDE_FORMATTER ax2.yformatter = LATITUDE_FORMATTER ax2.gridlines(draw_labels=["x", "y", "bottom", "left", "geo"]) ax2.add_feature(cfeature.LAND,facecolor="none",edgecolor='black',linewidths=1.5,zorder=2) #ax2.set_extent([-180, -120, -35, 0],crs=ccrs.PlateCarree()) plt.title('RMSE for W {} for GraphOP'.format(week+1)) ax4 = fig.add_subplot(2,9,(1,3), projection=ccrs.Robinson(central_longitude=180)) im4 = ax4.contourf(S2S_lon,S2S_lat,RMSE_S2S,np.linspace(0,7,40),transform=ccrs.PlateCarree(),vmin=0,vmax=7,cmap=newcmp,extend='both') ax4.xformatter = LONGITUDE_FORMATTER ax4.yformatter = LATITUDE_FORMATTER ax4.gridlines(draw_labels=["x", "y", "bottom", "left", "geo"]) ax4.add_feature(cfeature.LAND,facecolor="none",edgecolor='black',linewidths=1.5,zorder=2) #ax4.set_extent([-180, -120, -35, 0],crs=ccrs.PlateCarree()) plt.title('RMSE for week {} for S2S'.format(week+1)) #newcolors = np.vstack((bottom(np.linspace(0, 1, 128)[::-1]),top(np.linspace(0, 1, 128)[::-1]))) newcmp = ListedColormap(newcolors, name='OrangeBlue') RMSE1 = RMSE_GraphFT - RMSE_GraphOP ax5 = fig.add_subplot(2,9,(15,17), projection=ccrs.Robinson(central_longitude=180)) im5 = ax5.contourf(S2S_lon,S2S_lat,RMSE1,np.linspace(-10,1,40),transform=ccrs.PlateCarree(),vmin=-10,vmax=1,cmap=newcmp,extend='both') ax5.xformatter = LONGITUDE_FORMATTER ax5.yformatter = LATITUDE_FORMATTER ax5.gridlines(draw_labels=["x", "y", "bottom", "left", "geo"]) ax5.add_feature(cfeature.LAND,facecolor="none",edgecolor='black',linewidths=1.5,zorder=2) plt.title('RMSE difference for W {} for GraphFT - GraphOP'.format(week+1),y=1) RMSE2 = RMSE_GraphFT - RMSE_S2S ax6 = fig.add_subplot(2,9,(11,13), projection=ccrs.Robinson(central_longitude=180)) im6 = ax6.contourf(S2S_lon,S2S_lat,RMSE2,np.linspace(-10,1,40),transform=ccrs.PlateCarree(),vmin=-10,vmax=1,cmap=newcmp,extend='both') ax6.xformatter = LONGITUDE_FORMATTER ax6.yformatter = LATITUDE_FORMATTER ax6.gridlines(draw_labels=["x", "y", "bottom", "left", "geo"]) ax6.add_feature(cfeature.LAND,facecolor="none",edgecolor='black',linewidths=1.5,zorder=2) plt.title('RMSE difference of W {} for GraphFT - S2S'.format(week+1),y=1) c1 = fig.colorbar(im1,ax=ax1,shrink=0.6) c2 = fig.colorbar(im2,ax=ax2,shrink=0.6) #c3 = fig.colorbar(im3,ax=ax3,shrink=0.6) c4 = fig.colorbar(im4,ax=ax4,shrink=0.6) c5 = fig.colorbar(im5,ax=ax5,shrink=0.6) c6 = fig.colorbar(im6,ax=ax6,shrink=0.6) #if remove_seasonality: # plt.savefig('/home/vsansi01/results_RMSE_{}_{}WeeksAhead_rf.jpeg'.format(considered_variable,week+1),dpi=300) # print('-----> Plots saved')