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 data_mean,bootstrapping_CI,computing_confidence_intervals,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 os.chdir('/home/vsansi01/R_boxes/') def reading_boxes(files): xx = pd.DataFrame(['{:.2f}'.format(np.load(file)) for file in files]).to_excel('/home/vsansi01/output.xlsx') for file in files: print(file, '{:.4f}'.format(np.load(file))) #print('-----> Choose model, must be S2S, GraphFT or GraphOP') #model = str(input()) print('-----> Choose variable, must be u10, v10, t2m or tp') considered_variable = str(input()) print('-----> Choose how many weeks ahead, must be 2, 3 or 4') weeksahead = str(input()) #files_S2S = glob.glob(f'mean_value_{model}_{considered_variable}_{weeksahead}week' + '*',recursive=True) #xx = pd.DataFrame(['{:.2f}'.format(np.load(file)) for file in files]).to_excel('/home/vsansi01/output.xlsx') #for file in files: # print(file, '{:.4f}'.format(np.load(file))) files_S2S = sorted(glob.glob(f'mean_value_S2S_{considered_variable}_{weeksahead}week' + '*',recursive=True)) files_GraphFT = sorted(glob.glob(f'mean_value_GraphFT_{considered_variable}_{weeksahead}week' + '*',recursive=True)) files_S2S_ensemble = sorted(glob.glob(f'mean_value_S2S_ensemble_{considered_variable}_{weeksahead}week' + '*',recursive=True)) files_error_S2S = sorted(glob.glob(f'boxe_S2S_{considered_variable}_{weeksahead}week' + '*',recursive=True)) files_error_GraphFT = sorted(glob.glob(f'boxe_GraphFT_{considered_variable}_{weeksahead}week' + '*',recursive=True)) files_error_S2S_ensemble = sorted(glob.glob(f'boxe_S2S_ensemble_{considered_variable}_{weeksahead}week' + '*',recursive=True)) if weeksahead in [2,3]: files_GraphOP = sorted(glob.glob(f'mean_value_GraphOP_{considered_variable}_{weeksahead}week' + '*',recursive=True)) #reading_boxes(files_S2S_ensemble) for file1,file2,error1,error2 in zip(files_error_S2S_ensemble,files_error_GraphFT,bootstrapping_CI(files_error_S2S_ensemble,R_maps=True),bootstrapping_CI(files_error_GraphFT,R_maps=True)): print(f'{file1}, error = {round(error1,4)}')#, {file2}, error = {round(error2,4)}') #acc_values_week_ensemble_mean = data_mean(files_S2S_ensemble) #error_bar_week_ensemble_mean = bootstrapping_CI(files_S2S_ensemble) geographical_boxes = ('Northen Hemisphere','Tropical band','Southern Hemisphere','South Pole',\ 'North Pole','North Atlantic','Hawaii Archipelagos','French Polynesia',\ 'Europe','North America','Northeastern China - Japan','Indian Subcontinent',\ 'Tibetan Plateau') R_S2S = [0.18,0.32,0.09,0.0424,0.22,0.18,0.32,0.41,0.17,0.21,0.2,0.21,0.34] R_GraphFT = [0.35,0.33,0.24,0.31,0.41,0.33,0.5,0.45,0.36,0.32,0.27,0.36,0.48] R_S2S_ensemble = [0.30,0.42,0.19,0.13,0.34,0.24,0.53,0.50,0.21,0.31,0.22,0.32,0.52] R_error_S2S = [0.029,0.043,0.024,0.034,0.044,0.084,0.0204,0.0203,0.0076,0.0081,0.0071,0.0196,0.0335] R_error_GraphFT = [0.026,0.035,0.029,0.031,0.034,0.0096,0.0162,0.0096,0.0068,0.0066,0.0059,0.0171,0.0453] R_error_S2S_ensemble = [0.0034,0.0044,0.0027,0.0032,0.0049,0.008,0.0232,0.018,0.0074,0.0085,0.0076,0.0247,0.0438] assert len(R_S2S) == len(R_GraphFT) and len(R_S2S) == len(geographical_boxes), 'Not the same length' acc_values_week = { 'S2S': R_S2S, 'GraphFT': R_GraphFT } error_bar_week = { 'S2S': R_error_S2S, 'GraphFT': R_error_GraphFT } #x = np.arange(len(geographical_boxes)) # the label locations #width = 0.2 # the width of the bars #multiplier = 0 #fig, ax = plt.subplots(layout='constrained') #for attribute, measurement in acc_values_week.items(): # print(attribute) # print(measurement) # offset = width * multiplier #if attribute == 'S2S': #rect = ax.bar(x + offset, acc_values_week_ensemble_mean, width, label='S2S ensemble mean', alpha=0.2, color = 'b') #yerr_S2S = np.array([computing_confidence_intervals(r,0.95) for r in measurement]) #yerr_S2S = np.swapaxes(yerr_S2S,0,1) #plt.errorbar(x + offset,acc_values_week_ensemble_mean, yerr=error_bar_week_ensemble_mean, ecolor='k',fmt = 'None', capsize=3, capthick=0.5) #yerr = np.array([computing_confidence_intervals(r,0.95) for r in measurement]) #yerr = np.swapaxes(yerr,0,1) #plt.errorbar(x + offset,measurement, yerr=errors, ecolor='k',fmt = 'None', capsize=3, capthick=0.5) # rects = ax.bar(x + offset, measurement, width, label=attribute) # multiplier += 1 # Add some text for labels, title and custom x-axis tick labels, etc. #ymax = 0.7 if weeksahead == 2 else 0.4 #ax.set_ylabel('R') #ax.set_xlabel('Geographical boxes') #ax.set_xticks(x + width, geographical_boxes) #ax.legend(loc='upper right') #ax.set_ylim(0, 0.6) #plt.xticks(rotation=90) #plt.savefig(f'/home/vsansi01/Barplot_R_geographical_boxes_test.jpeg',dpi=300) x = np.arange(len(geographical_boxes)) # the label locations width = 0.2 # the width of the bars multiplier = 0 fig, ax = plt.subplots(figsize=(10, 6), layout='constrained') # Adjust figure size colors = plt.cm.get_cmap('tab10') # Using a color-blind-friendly colormap colors = [colors(0),colors(3)] # Plot bars with error bars (if available) and improve the aesthetics for i, ((attribute, measurement), (att, error)) in enumerate(zip(acc_values_week.items(),error_bar_week.items())): if attribute == 'GraphFT': continue else: offset = width * multiplier rects = ax.bar(x + offset, measurement, width, label=attribute, color=colors[i], alpha=1) if attribute == 'S2S': rect = ax.bar(x + offset, R_S2S_ensemble, width, label='S2S ensemble mean', alpha=0.2, color = 'b') #yerr_S2S = np.array([computing_confidence_intervals(r,0.95) for r in measurement]) #yerr_S2S = np.swapaxes(yerr_S2S,0,1) plt.errorbar(x + offset, R_S2S_ensemble, yerr=R_error_S2S_ensemble, ecolor='k',fmt = 'None', capsize=3, capthick=0.5) # Optionally, add error bars here (uncomment and customize if you want to add): #plt.errorbar(x + offset,measurement, yerr=errors, ecolor='k',fmt = 'None', capsize=3, capthick=0.5) ax.errorbar(x + offset, measurement, yerr=error, fmt='none', capsize=3, color='black') # Annotate the bars with their values #for rect in rects: # height = rect.get_height() + 0.02 # ax.annotate(f'{height:.2f}', # xy=(rect.get_x() + rect.get_width() / 2, height), # xytext=(0, 3), # 3 points vertical offset # textcoords="offset points", # ha='center', va='bottom', fontsize=8) multiplier += 1 # Add labels, title, and custom x-axis tick labels ax.set_ylabel('R', fontsize=12) ax.set_xlabel('Geographical Boxes', fontsize=12) #ax.set_title('R Values Across Geographical Boxes', fontsize=14, weight='bold') ax.set_xticks(x + width, geographical_boxes) ax.set_ylim(0, max([max(v) for v in acc_values_week.values()]) + 0.1) # Dynamic y-axis limit ax.legend(loc='upper left')#, bbox_to_anchor=(1, 1), fontsize=10) # Move legend outside the plot plt.xticks(rotation=90, fontsize=10) # Add gridlines to improve readability ax.grid(True, axis='y', linestyle='--', alpha=0.7) # Save the plot with a higher resolution plt.tight_layout() plt.savefig('/home/vsansi01/Barplot_R_geographical_boxes_readable.jpeg', dpi=300)