#!/usr/bin/env python3 import matplotlib.pyplot as plt import data_norcpm import data_hadisst import numpy as np import os from scipy.stats.stats import pearsonr import sys #plt.plot(range(5)) #plt.ylabel('some numbers') #plt.show() def mask2d_by_latlon(lats,lons,latbe,lonbe): ## create mask2d(nlat,nlon) array from axis(lats,lons) ## and set False in latbe,lonbe ## need to deal with lon flip nlat,nlon = len(lats),len(lons) mask2d = np.full((nlat,nlon),False) ## True: make value = filled_value llats = [ not min(latbe)<=i<=max(latbe) for i in lats] ## need more coding for -180~180 or 0~360 lons = [ i if i>0 else i+360. for i in lons] lonbe = [ i if i>0 else i+360. for i in lonbe] llons = [ not min(lonbe)<=i<=max(lonbe) for i in lons] for i in range(nlat): mask2d[i,:] = np.logical_or(mask2d[i,:],llons) for i in range(nlon): mask2d[:,i] = np.logical_or(mask2d[:,i],llats) return mask2d def weight2d_by_latlon(lats,lons): nlat,nlon = len(lats),len(lons) w2d = np.full((nlat,nlon),0.) w1d = np.cos(np.deg2rad(lats)) for i in range(nlon): w2d[:,i] = w1d return w2d def varTLL_ts(varTLL,mask2d='',latbe='',lonbe=''): ## cal area mean time series ## from 3d(TLL) var ## should be dims in varTLL.dims[] #if (not mask2d) and (latbe and lonbe): # mask2d = mask2d_by_latlon(varTLL.dims[1]['value'],varTLL.dims[2]['value'],latbe,lonbe) ## need more check on arguments ## weighting weight2d = weight2d_by_latlon(varTLL.dims[1]['value'],varTLL.dims[2]['value']) ts = list() maskTLL = np.ma.getmaskarray(varTLL) for i in range(varTLL.shape[0]): tslice = varTLL[i,] #print(np.count_nonzero(maskTLL[i,])) tslice.mask = np.logical_or(maskTLL[i,],mask2d) ts.append(np.ma.average(tslice,weights=weight2d)) return np.array(ts) def remove_seasonal_cycle(var1d,months,clm12ts=''): ## for monthly data, ## var1d: data time series ## months: data time, must be same dim as var1d ## if months all same, it becomes remove mean ## clm12ts: 1D climatology, should be 12 elements if type(clm12ts) != type('') and len(clm12ts) != 12: print('remove_seasonal_cycle(): len(clm12ts) != 12, exit... ') return False scrVar1d = np.array(var1d) ## incase var1d is list testclmts = list() for m in range(1,12+1): need = [ j for j in range(len(months)) if months[j] == m ] if not need: continue if type(clm12ts) != type(''): testclmts.append(clm12ts[m-1]) scrVar1d[need] -= clm12ts[m-1] else: testclmts.append(np.mean(scrVar1d[need])) scrVar1d[need] -= np.mean(scrVar1d[need]) if testclmts: testclmtsavg = np.mean(testclmts) #print([f'{i-testclmtsavg:.2f}' for i in testclmts]) return scrVar1d path = './test_set/ana_19800115_me_hindcasts' path = '/projects/NS9039K/shared/norcpm/cases/NorCPM/NorCPM_V1/ana_19800115_me_hindcasts' prefix = 'ana_19800115_me' path2 = '/cluster/shared/NS9039K/norcpm_ana_hindcast/hindcasts' prefix2 = 'norcpm_ana_hindcast' histpath = '/cluster/shared/NS9039K/archive/noresm_ctl_f09_tn14_19700101' leadmonth = 8 initmonths = [10] var, component = 'SST', 'atm' ## region setting tag, latbe, lonbe = 'iod_west', [-10.,10.], [50,70] tag, latbe, lonbe = 'nino3.4', [-5.,5.], [-170,-120] tag, latbe, lonbe = 'atl3', [-5.,5.], [-20,0] clmyb = 1985 clmye = 2010 initMon = [ 1, 4, 7, 10] initmonName = ['Jan','Apr','Jul','Oct'] def norcpm_hindcasts_ensmean_region_ts_cor(path,var='SST',component='atm',tag='atl3',latbe=[-5.,5.],lonbe=[-20,0],leadmonth=0,inityears=[],initmonths=[],prefix='',ploting=False): ## compare norcpm hindcasts SST regional time series to obs (HadISST) ## return correlation and persistance hindcasts correlation ## get norcpm hindcasts in ts cachefn = f'data/{var}_{component}_{tag}_leadMonth{leadmonth:02}.npy' if os.path.isfile(cachefn): with open(cachefn,'rb') as f: ts = np.ma.array(np.load(f)) vtime = np.ma.array(np.load(f)) ts_years = np.load(f) ts_months = np.load(f) else: b = data_norcpm.norcpm_forecasts(path,prefix) v = b.get_monthly_variable_ensmean_leadmonth(var,component,leadmonth=leadmonth,months=initmonths,years=inityears) vtime = v.time ts_years = [ int(i[0:4]) for i in vtime ] ts_months = [ int(i[4:6]) for i in vtime ] m2d = mask2d_by_latlon(v.dims[1]['value'],v.dims[2]['value'],latbe,lonbe) ts = varTLL_ts(v,m2d) with open(cachefn,'wb') as f: np.save(f,list(ts)) np.save(f,list(vtime)) np.save(f,ts_years) np.save(f,ts_months) ## obs: HadISST obs = data_hadisst.hadisst() ## obs ts with time of norcpm obsv = obs.read_monthly_data(yyyymm=vtime) m2d = mask2d_by_latlon(obsv.dims[1]['value'],obsv.dims[2]['value'],latbe,lonbe) obsts = varTLL_ts(obsv,m2d) ## obs persistant forecast perv = obs.read_monthly_data_persistant_leadmonth(yyyymm=vtime,leadmonth=leadmonth) pervts = varTLL_ts(perv,m2d) ## remove annual cycle scrts = remove_seasonal_cycle(ts,ts_months) scrobsts = remove_seasonal_cycle(obsts,ts_months) scrpervts = remove_seasonal_cycle(pervts,ts_months) ## correlation #r = np.corrcoef(scrts,scrobsts) ## Pearson product-moment correlation coefficients #r_per = np.corrcoef(scrpervts,scrobsts) ## Pearson product-moment correlation coefficients #r = np.correlate(scrts,scrobsts) ## cross-correlation #r_per = np.correlate(scrpervts,scrobsts) ## cross-correlation r = pearsonr(scrts,scrobsts) ## Pearson correlation coefficient r_per = pearsonr(scrpervts,scrobsts) ## Pearson correlation coefficient rmse = np.sqrt(((scrts-scrobsts)**2).mean()) rmse_per = np.sqrt(((scrpervts-scrobsts)**2).mean()) ## ploting if ploting: plt.plot(scrts,color='red',label=f'hindcast,r={r[0,1]:.2f}') plt.plot(scrpervts,color='green',label=f'persis,r={r_per[0,1]:.2f}') plt.plot(scrobsts,color='black',label='obs') plt.legend() plt.title(tag) plt.show() return r,r_per,rmse,rmse_per def norcpm_hindcast_ensmean_region_ts_bias(path,var='SST',component='atm',tag='atl3',latbe=[-5.,5.],lonbe=[-20,0],inityear=0,initmonth=0,prefix='',ploting=False,remove_hist_clm12=False,difobs=False,withobs=False): ## compare norcpm 1 hindcast SST regional time series to obs (HadISST) ## return difference of time series ## get norcpm hindcasts in ts ## not done yet fnmonth = 12 #cachefn = f'data/{var}_{component}_{tag}_{inityear}-{initmonth:02}-15.npy' cachefn = f'data/{var}_{component}_{tag}_{inityear}-{initmonth:02}-15.npy' if os.path.isfile(cachefn): with open(cachefn,'rb') as f: ts = np.ma.array(np.load(f)) vtime = np.ma.array(np.load(f)) ts_years = np.load(f) ts_months = np.load(f) dims = np.load(f,allow_pickle=True) else: b = data_norcpm.norcpm_forecasts(path,prefix) v = b.get_monthly_variable_ensmean_leadmonth(var,component,leadmonth=0,months=[initmonth],years=[inityear]) dims = v.dims vtime = v.time for i in range(1,fnmonth+1): ## from 0 to nmonth v1 = b.get_monthly_variable_ensmean_leadmonth(var,component,leadmonth=i,months=[initmonth],years=[inityear]) vtime = np.ma.append(vtime,v1.time) v = np.ma.append(v,v1,axis=0) v.dims = dims ts_years = [ int(i[0:4]) for i in vtime ] ts_months = [ int(i[4:6]) for i in vtime ] m2d = mask2d_by_latlon(dims[1]['value'],dims[2]['value'],latbe,lonbe) ts = varTLL_ts(v,m2d) with open(cachefn,'wb') as f: np.save(f,list(ts)) np.save(f,list(vtime)) np.save(f,ts_years) np.save(f,ts_months) np.save(f,dims) ## remove annual cycle if remove_hist_clm12: ## the m2d is not saved, need rewrite whole data model hist = data_norcpm.norcpm_data(histpath) v12 = hist.get_ensmean_clm12(var,component,years=list(range(clmyb,clmye+1)),tag=f'hist_clm12_{clmyb}-{clmye}') v12.dims = dims m2d = mask2d_by_latlon(dims[1]['value'],dims[2]['value'],latbe,lonbe) v12ts = varTLL_ts(v12,m2d) #print('ts:') #print(ts) #print('clm12:') #print(v12ts) ts = remove_seasonal_cycle(ts,ts_months,clm12ts=v12ts) #print('scr:') #print(ts) #return False else: ## K to degC ts = ts-273.15 if difobs or withobs: ## obs: HadISST obs = data_hadisst.hadisst() ## obs ts with time of norcpm obsv = obs.read_monthly_data(yyyymm=vtime) m2d = mask2d_by_latlon(obsv.dims[1]['value'],obsv.dims[2]['value'],latbe,lonbe) obsts = varTLL_ts(obsv,m2d) ## remove obs annual cycle if remove_hist_clm12: if difobs or withobs: obsclm12 = obs.read_clm12('sst',years=list(range(clmyb,clmye+1))) obsclm12.dims = obsv.dims obsclm12ts = varTLL_ts(obsclm12,m2d) #print('obsts:') #print(obsts) #print('obsclm12ts:') #print(obsclm12ts) obsts = remove_seasonal_cycle(obsts,ts_months,clm12ts=obsclm12ts) #print('scr obsts:') #print(obsts) if difobs: ts = ts-obsts ## ploting if ploting: plt.plot(ts,color='black',label=f'hindcast') #plt.legend() plt.title(f'hindcast bias(scr) {inityear}-{initmonth:02}') plt.show() if withobs: return ts,obsts return ts def norcpm_get_r_p(): cachefn = f'norcpm_1vs2_{tag}.npy' if os.path.isfile(cachefn): with open(cachefn,'rb') as f: r = np.load(f) p = np.load(f) rmse = np.load(f) else: r = np.empty([13,4,3]) ## leadmonth(0-12), initmonth, norcpm 1&2 and norcpm1 sst (atm) p = np.empty([13,4,3]) ## need to use len() rmse = np.empty([13,4,3]) for j in initMon: for i in range(13): r_norcpm, r_persis1, rmse_norcpm, rmse_persis = norcpm_hindcasts_ensmean_region_ts_cor(var='TS',path=path,prefix=prefix,tag=f'{tag}_im{j:02}',latbe=latbe,lonbe=lonbe,initmonths=[j],leadmonth=i) #r_norcpmSST, r_persis1 = norcpm_hindcasts_ensmean_region_ts_cor(var='SST',path=path,prefix=prefix,tag=f'{tag}_im{j:02}',latbe=latbe,lonbe=lonbe,initmonths=[j],leadmonth=i) r_norcpm2, r_persis2, rmse_norcpm2, rmse_persis = norcpm_hindcasts_ensmean_region_ts_cor(var='TS',path=path2,prefix=prefix2,tag=f'norcpm2_{tag}_im{j:02}',latbe=latbe,lonbe=lonbe,initmonths=[j],leadmonth=i,inityears=list(range(1985,2010+1))) #print(f'lm={i:02}/im={j} R of NorCPM_V1: {r_norcpm[0]:.2f}, NorCPM_V1 SST: {r_norcpmSST[0]:.2f}, NorCPM2: {r_norcpm2[0]:.2f}, persis: {r_persis1[0]:.2f}') print(f'lm={i:02}/im={j} R of NorCPM_V1: {r_norcpm[0]:.2f}, NorCPM2: {r_norcpm2[0]:.2f}, persis: {r_persis1[0]:.2f}') r[i,initMon.index(j),:] = r_persis1[0],r_norcpm[0],r_norcpm2[0] #,r_norcpmSST[0] p[i,initMon.index(j),:] = r_persis1[1],r_norcpm[1],r_norcpm2[1] #,r_norcpmSST[1] rmse[i,initMon.index(j),:] = rmse_persis,rmse_norcpm,rmse_norcpm2 with open(cachefn,'wb') as f: np.save(f,r) np.save(f,p) np.save(f,rmse) return r,p,rmse def plotting_4_cor(title,fn,r,ymin,ymax,p=False): fig, ax = plt.subplots(2,2,figsize=( 6.4 *1.2, 4.8 *1.2)) pthres = 0.05 monName = list('JFMAMJJASOND') for i in ax.ravel(): i.set_ylim(ymin,ymax) i.set_xticks(list(range(12+1))) colors = ['black' ,'green', 'red'] labels = ['Persis','NorCPM_V1','NorCPM2'] for i,lm in enumerate(initMon): ilm = initMon.index(lm) ii,jj = int(i/2),i%2 ax[ii,jj].set_title(f' init {initmonName[ilm]}') for j in range(3): ax[ii,jj].plot(r[:,ilm,j],color=colors[j],label=labels[j]) ax[ii,jj].set_xticklabels([monName[(i+lm-1)%12] for i in range(13)]) if type(p) != type(False): ## add dot where p <= pthres for jp in range(13): if p[jp,ilm,j] > pthres: continue ax[ii,jj].plot(jp,r[jp,ilm,j],marker='o',color=colors[j]) if False: ## to del lm = 1 ilm = initMon.index(lm) ax[0,0].set_title(f' init {initmonName[ilm]}') ax[0,0].plot(r[:,ilm,0],color='black',label='Persis') ax[0,0].plot(r[:,ilm,1],color='green',label='NorCPM_V1') ax[0,0].plot(r[:,ilm,2],color='red',label='NorCPM2') ax[0,0].set_xticklabels([monName[(i+lm-1)%12] for i in range(13)]) lm = 4 ilm = initMon.index(lm) ax[0,1].set_title(f' init {initmonName[ilm]}') ax[0,1].plot(r[:,ilm,0],color='black',label='Persis') ax[0,1].plot(r[:,ilm,1],color='green',label='NorCPM_V1') ax[0,1].plot(r[:,ilm,2],color='red',label='NorCPM2') ax[0,1].set_xticklabels([monName[(i+lm-1)%12] for i in range(13)]) lm = 7 ilm = initMon.index(lm) ax[1,0].set_title(f' init {initmonName[ilm]}') ax[1,0].plot(r[:,ilm,0],color='black',label='Persis') ax[1,0].plot(r[:,ilm,1],color='green',label='NorCPM_V1') ax[1,0].plot(r[:,ilm,2],color='red',label='NorCPM2') ax[1,0].set_xticklabels([monName[(i+lm-1)%12] for i in range(13)]) lm = 10 ilm = initMon.index(lm) ax[1,1].set_title(f' init {initmonName[ilm]}') ax[1,1].plot(r[:,ilm,0],color='black',label='Persis') ax[1,1].plot(r[:,ilm,1],color='green',label='NorCPM_V1') ax[1,1].plot(r[:,ilm,2],color='red',label='NorCPM2') ax[1,1].set_xticklabels([monName[(i+lm-1)%12] for i in range(13)]) fig.suptitle(title) ax[0,0].legend() ax[0,1].set_yticklabels([]) ax[1,1].set_yticklabels([]) ax[1,0].set_xlabel('lead month') ax[1,1].set_xlabel('lead month') plt.savefig(fn,bbox_inches='tight',pad_inches=0) plt.show() def plotting_4_lines(title,fn,lines,ymin,ymax,colors,show=False): ## lines: [nfig,nline,y] dims = lines.shape fig, ax = plt.subplots(2,2,figsize=( 6.4 *1.2, 4.8 *1.2)) colorindex = ['black','red','blue'] monName = list('JFMAMJJASOND') for i in ax.ravel(): i.set_ylim(ymin,ymax) i.set_xticks(list(range(12+1))) for i,lm in enumerate(initMon): ilm = initMon.index(lm) ii = int(i/2) jj = i%2 ax[ii,jj].set_title(f' init {initmonName[ilm]}') for l in range(dims[1]): color = colorindex[int(colors[i,l])] ax[ii,jj].plot(lines[i,l,:],color=color) ax[ii,jj].set_xticklabels([monName[(i+lm-1)%12] for i in range(13)]) fig.suptitle(title) #ax[0,0].legend() ax[0,1].set_yticklabels([]) ax[1,1].set_yticklabels([]) ax[1,0].set_xlabel('lead month') ax[1,1].set_xlabel('lead month') plt.savefig(fn,bbox_inches='tight',pad_inches=0) if show: plt.show() def plot_bias(): yb = 1985 ye = 2010 ny = ye-yb+1 hindts = np.empty([4,ny,13]) obsts = np.empty([4,ny,13]) colorthreshold = .3 colors = np.zeros_like(obsts[:,:,0]) for im,Moninit in enumerate(initMon): Moninitname = initmonName[initMon.index(Moninit)] #plt.ylim(-2.5,2.5) for i,y in enumerate(range(yb,ye+1)): hindts[im,i,:],obsts[im,i,:] = norcpm_hindcast_ensmean_region_ts_bias(path2,var='TS',component='atm',tag=tag,latbe=latbe,lonbe=lonbe,inityear=y,initmonth=Moninit,prefix='norcpm_ana_hindcast',ploting=False,remove_hist_clm12=True,difobs=True,withobs=True) #print(','.join([f'{i:.2f}' for i in a])) m = obsts[im,i,10:].mean() if m > colorthreshold: colors[im,i] = 1 if m < -colorthreshold: colors[im,i] = -1 #plt.plot(hindts[y,:],color=color)#,label='Persis') #plt.title(f'{tag} hindcasts bias, init {Moninitname} (scr from hist, obs also scr)') #plt.savefig(f'{tag}_hindcasts_bias_init{Moninitname}.png') #plt.clf() #plt.show() plotting_4_lines(f'hindcast ensmean bias (scr) {tag}',f'bias_{tag}.png',hindts,-2.5,2.5,colors) def plot_hind_obs(): yb = 1985 ye = 2010 ny = ye-yb+1 hindts = np.empty([4,ny,13]) obsts = np.empty([4,ny,13]) remove_hist_clm12 = False colorthreshold = .3 colors = np.zeros_like(obsts[:,:,0]) for im,Moninit in enumerate(initMon): Moninitname = initmonName[initMon.index(Moninit)] #plt.ylim(-2.5,2.5) for i,y in enumerate(range(yb,ye+1)): hindts[im,i,:],obsts[im,i,:] = norcpm_hindcast_ensmean_region_ts_bias(path2,var='TS',component='atm',tag=tag,latbe=latbe,lonbe=lonbe,inityear=y,initmonth=Moninit,prefix='norcpm_ana_hindcast',ploting=False,remove_hist_clm12=remove_hist_clm12,difobs=False,withobs=True) #print(','.join([f'{i:.2f}' for i in a])) hindts[im,i,:] -= obsts[im,i,:] hindts[im,i,:] -= hindts[im,i,0] #m = obsts[im,i,10:].mean() #if m > colorthreshold: colors[im,i] = 1 #if m < -colorthreshold: colors[im,i] = -1 #plt.plot(hindts[y,:],color=color)#,label='Persis') #plt.title(f'{tag} hindcasts bias, init {Moninitname} (scr from hist, obs also scr)') #plt.savefig(f'{tag}_hindcasts_bias_init{Moninitname}.png') #plt.clf() #plt.show() title = f'hindcast ensmean bias {tag} (hind-obs)' if remove_hist_clm12: title += ' scr' plotting_4_lines(title,f'bias_testing_{tag}.png',hindts,-2.5,2.5,colors,show=True) def plot_hinds_obs_line1(initmonth,show=False): yb = 1985 ye = 2010 ny = ye-yb+1 hindts = np.empty([4,ny,13]) #initmonth = 1 ## 1,4,7,10 iim = initMon.index(initmonth) scr = True shift_init = True for im,Moninit in enumerate(initMon): Moninitname = initmonName[initMon.index(Moninit)] for i,y in enumerate(range(yb,ye+1)): hindts[im,i,:] = norcpm_hindcast_ensmean_region_ts_bias(path2,var='TS',component='atm',tag=tag,latbe=latbe,lonbe=lonbe,inityear=y,initmonth=Moninit,prefix='norcpm_ana_hindcast',ploting=False,remove_hist_clm12=scr,difobs=False,withobs=False) yyyymm = [] for y in range(yb,ye+1+1): yyyymm.extend([f'{y}{m:02}' for m in range(1,12+1)]) obs = data_hadisst.hadisst() obsv = obs.read_monthly_data(yyyymm=yyyymm) m2d = mask2d_by_latlon(obsv.dims[1]['value'],obsv.dims[2]['value'],latbe,lonbe) obsts = varTLL_ts(obsv,m2d) if scr: obsts = remove_seasonal_cycle(obsts,[ int(i[4:6]) for i in yyyymm]) fig = plt.figure(figsize=( 6.4 *4.0, 4.8 *1.2)) plt.plot(obsts,color='black') xtickvalue = list() xticklabels = list() shift = 0 for i,y in enumerate(range(yb,ye+1)): iib = yyyymm.index(f'{y}{initMon[iim]:02}') plt.plot(list(range(iib,iib+13)),hindts[iim,i,:],color='#808080') plt.plot(iib,hindts[iim,i,0],marker='o',color='#808080') if shift_init: shift = obsts[iib] - hindts[iim,i,0] ## shift init value to obs plt.plot(list(range(iib,iib+13)),hindts[iim,i,:]+shift,color='red') plt.plot(iib,hindts[iim,i,0]+shift,marker='o',color='red') xticklabels.append(f'{y}{initMon[iim]:02}') xtickvalue.append(iib) title = f'hindcasts ensmean {tag} initmonth={initmonth}' if scr: title += ' scr' if shift_init: title += ' init shifted' plt.title(title) plt.xticks(xtickvalue,xticklabels) scrs = '' if scr: scrs = '_scr' fn = f'hindcasts_ensmean{scrs}_{tag}_init{initmonth:02}.png' #plt.show() plt.savefig(fn,bbox_inches='tight',pad_inches=0) #plotting_4_lines(title,f'bias_testing2_{tag}.png',hindts,-2.5,2.5,colors,show=True) for i in initMon: plot_hinds_obs_line1(i) #r,p,rmse = norcpm_get_r_p() ## leadmonth, initmonth, norcpm 1&2 and norcpm1 sst (atm) ## p is two-side p-value #plotting_4_cor(title=f'Ensmean correlation, {tag}. dot: p-val<0.05',fn=f'cor_{tag}.png',r=r,p=p,ymin=-0.2,ymax=1.0) #plotting_4_cor(title=f'Ensmean RMSE, {tag}',fn=f'rmse_{tag}.png',r=rmse,ymin=0.0,ymax=1.0)