# -*- coding: utf-8 -*- """ Created on Tue Feb 7 09:43:41 2017 @author: sonja """ #============================================================================== # Start forecast model #============================================================================== # -*- coding: utf-8 -*- # needed imports import numpy as np import operator #for detrending the data import matplotlib matplotlib.rcParams['backend'] = "Qt4Agg" #from datetime import datetime import scipy import os #needed for loop over files from sklearn.cross_decomposition import CCA import functions #============================================================================== # Idea #============================================================================== # x(t) state vector (e.g. DJF precip) # x_i clusteroids # y(t) predictor (e.g. sst) #y_i predictor corresponding to x_i, obtained by averaging or regression such that |y_i|=1 #Given y(t_OCT), how do we predict x(t_DJF)? #define projection coeffs: #\alpha_i=Y_iT*y(t) #predict: x(t) = (Sum(\alpha_i * X_i) / (sum(\alpha_j) ) ) #Idea: using all data except for last year to predict last year: # 5 precursors: #snow cover (Eurasia) #ice (Eurasia) #SST (tropical pacific, north atl, mediterranean) #GPH #SLP #============================================================================== # Procedure #============================================================================== #1.)calculate clusters without year, which we want to predict #2.) calculate composites from clusters #3.) normalize composites #4.) add all composites together #5.) define projection coeffs #6.)predict X(t) for that year #7.) calculate pattern regression #8.) after all years were predicted, calculate time regression print ("\n") print ("Start Forecast model") print ("\n") #Load temperature temperature= functions.forecast_variable("temperature") # Load precursor data precursors_temperature=functions.precursors("Precursors for temperature forecast") # reshape precursors to 1D precursors_temperature_1D = functions.precursors_1D("1D precursors for temperature forecast",precursors_temperature) # parameters for forecast paramsForecast_temp=functions.paramsForecast("temperature forecast",0,50,1967,temperature.tmpDataAllYearsObs_1D,temperature.tmpData2) #possible precursors for which forecast should be calculated saveArray=["SIC","SCE","SST_tropics","SST_atl","SST_medi","SLP_nh","GPH500_EU"] #select from saveArray which precursors should be used to calculate forecast #0="SIC", 1="SCE", 2="SST_tropics",3="SST_atl",4="SST_medi",5="SLP_nh",6="GPH500_EU" listOfprecursors=[2,6] # begYear = 0 endYear = 50 yearStart=1967 clusterarray = np.zeros((2,temperature.tmpDataAllYearsObs_1D.shape[0]-1))#historic reasons for choice of dimensions patternCorr = np.zeros(endYear-begYear) #possible precursors for which forecast should be calculated saveArray=["SIC","SCE","SST_tropics","SST_atl","SST_medi","SLP_nh","GPH500_EU"] #select from saveArray which precursors should be used to calculate forecast #0="SIC", 1="SCE", 2="SST_tropics",3="SST_atl",4="SST_medi",5="SLP_nh",6="GPH500_EU" listOfprecursors=[2,6] # parameters for forecast paramsForecast_temp=functions.paramsForecast("temperature forecast",0,50,1967,temperature.tmpDataAllYearsObs_1D,temperature.tmpData2) for year in range(begYear,endYear): prcp_weights_training_Y = np.delete(temperature.tmpDataAllYearsObs_1D,year,axis=0) #remove 1 year & standardize data seaICData_1DSub = np.delete(precursors_temperature_1D.seaICData2_1D,year,axis=0) seaceData_1D_longTerm = precursors_temperature_1D.seaICData2_1D[year] seaICData_1DSub,seaceData_1D_longTerm=functions.calculate_standardized_precursors(seaICData_1DSub,seaceData_1D_longTerm) snowCEData2_1DSub = np.delete(precursors_temperature_1D.snowCEData2_1D,year,axis=0) snowCEData_1D_longTerm = precursors_temperature_1D.snowCEData2_1D[year] snowCEData2_1DSub,snowCEData_1D_longTerm=functions.calculate_standardized_precursors(snowCEData2_1DSub,snowCEData_1D_longTerm) mediData2_1DSub = np.delete(precursors_temperature_1D.seastData_Mediterranean2_1D,year,axis=0) mediData_1D_longTerm = precursors_temperature_1D.seastData_Mediterranean2_1D[year] mediData2_1DSub,mediData_1D_longTerm=functions.calculate_standardized_precursors(mediData2_1DSub,mediData_1D_longTerm) atlData2_1DSub = np.delete(precursors_temperature_1D.seastData_North_Atl2_1D,year,axis=0) atlData_1D_longTerm = precursors_temperature_1D.seastData_North_Atl2_1D[year] atlData2_1DSub,atlData_1D_longTerm=functions.calculate_standardized_precursors(atlData2_1DSub,atlData_1D_longTerm) tropData2_1DSub = np.delete(precursors_temperature_1D.seastData_South_Pac2_1D,year,axis=0) tropData_1D_longTerm = precursors_temperature_1D.seastData_South_Pac2_1D[year] tropData2_1DSub,tropData_1D_longTerm=functions.calculate_standardized_precursors(tropData2_1DSub,tropData_1D_longTerm) slpData2_1DSub = np.delete(precursors_temperature_1D.sealpData2_1D,year,axis=0) slpData_1D_longTerm = precursors_temperature_1D.sealpData2_1D[year] slpData2_1DSub,slpData_1D_longTerm=functions.calculate_standardized_precursors(slpData2_1DSub,slpData_1D_longTerm) gph500Data2_1DSub = np.delete(precursors_temperature_1D.geophp500Data2_1D,year,axis=0) gph500Data_1D_longTerm = precursors_temperature_1D.geophp500Data2_1D[year] gph500Data2_1DSub,gph500Data_1D_longTerm=functions.calculate_standardized_precursors(gph500Data2_1DSub,gph500Data_1D_longTerm) #concatenate data dataX=(\ seaICData_1DSub,\ snowCEData2_1DSub,\ tropData2_1DSub,\ atlData2_1DSub,\ mediData2_1DSub,\ slpData2_1DSub,\ gph500Data2_1DSub,\ ) dataForecast = (\ seaceData_1D_longTerm,\ snowCEData_1D_longTerm,\ tropData_1D_longTerm,\ atlData_1D_longTerm,\ mediData_1D_longTerm,\ slpData_1D_longTerm,\ gph500Data_1D_longTerm,\ ) #merge only those precursors which are selected if len(listOfprecursors)==1: All_Data_1D_training_X=np.array(dataX[listOfprecursors[0]]) All_Data_1DForecast=np.array(dataForecast[listOfprecursors[0]]) else: All_Data_1D_training_X=np.concatenate(operator.itemgetter(*listOfprecursors)(dataX), axis=1) All_Data_1DForecast = np.concatenate(operator.itemgetter(*listOfprecursors)(dataForecast), axis=0) #training cca = CCA(n_components=4) cca.fit(All_Data_1D_training_X, prcp_weights_training_Y) All_Data_1D_training_X_transform, prcp_weights_training_Y_transform = cca.fit_transform(All_Data_1D_training_X, prcp_weights_training_Y) #prediction forecast_precip2=cca.predict(All_Data_1DForecast.reshape(1, -1)) forecast_precip = forecast_precip2.reshape(-1,1)[:,0] paramsForecast_temp.tmpDataAllYearsPredict[year]=forecast_precip #pattern correlation paramsForecast_temp.patternCorr[year-begYear] = scipy.stats.pearsonr(forecast_precip, temperature.tmpDataAllYearsObs_1D[year])[0] # ''' print ("\n") print ("Finished year "+str(1967+year)+ "") print ("\n") #correlation paramsForecast_temp.tmpDataAllYearsObs_1D_Corr=functions.calculate_time_correlation(temperature.tmpDataAllYearsObs_1D,paramsForecast_temp) tmpDataAllYearsObs_1D_Corr_reshape = np.reshape(paramsForecast_temp.tmpDataAllYearsObs_1D_Corr,(temperature.tmpData2.shape[1],temperature.tmpData2.shape[2])) print(np.mean(paramsForecast_temp.patternCorr[0:paramsForecast_temp.end_year-paramsForecast_temp.beg_year])) print('cross-correlation:') tmpDataAllYearsObs_1D_Corr_reshape = np.reshape(paramsForecast_temp.tmpDataAllYearsObs_1D_Corr,(temperature.tmpData2.shape[1],temperature.tmpData2.shape[2])) tmpCorr = tmpDataAllYearsObs_1D_Corr_reshape.copy() print(np.nanmean(tmpDataAllYearsObs_1D_Corr_reshape)) print('pattern correlation:') # calculate time correlation for each point in forecast model paramsForecast_temp.tmpDataAllYearsObs_1D_Corr=functions.calculate_time_correlation(temperature.tmpDataAllYearsObs_1D,paramsForecast_temp) tmpDataAllYearsObs_1D_Corr_reshape = np.reshape(paramsForecast_temp.tmpDataAllYearsObs_1D_Corr,(temperature.tmpData2.shape[1],temperature.tmpData2.shape[2])) print(np.mean(paramsForecast_temp.patternCorr[0:paramsForecast_temp.end_year-paramsForecast_temp.beg_year])) #define outputname if len(listOfprecursors)==1: saveString=saveArray[listOfprecursors[0]] else: saveString=join(operator.itemgetter(*listOfprecursors)(saveArray),"-") a_perce_reshape = np.reshape(paramsForecast_temp.a_perce,(temperature.tmpData2.shape[1],temperature.tmpData2.shape[2])) file_path = str(len(listOfprecursors))+'-precursor/' directory = os.path.dirname(file_path) try: os.stat(directory) except: os.mkdir(directory) #save plot #functions.save_plot(file_path+saveString+'-1967-2016_CCA.png',tmpDataAllYearsObs_1D_Corr_reshape,a_perce_reshape,temperature.tempParams,temperature.lons,temperature.lats) #save time correlation #file = open(saveString+'-1967-2016_CCA.txt','w') #file.write(str(np.nanmean(tmpDataAllYearsObs_1D_Corr_reshape).round(5))) #save spatial correlation #functions.save_variable(file_path+'/pattern_time_correlation_GHCN_CAMS_'+saveString+'1967_2016_CCA.nc',tmpDataAllYearsObs_1D_Corr_reshape,temperature.lat_grid,temperature.lon_grid,'time_correlation') # save correlation significance #functions.save_variable(file_path+'/pattern_time_correlation_GHCN_CAMS_'+saveString+'1967_2016_CCA_significance.nc',tmpDataAllYearsObs_1D_Corr_reshape,temperature.lat_grid,temperature.lon_grid,'a_perce_reshape')