#!/usr/bin/env python ###------------------------------------- ### Choix Type Image (none, png, eps) ###------------------------------------- img='png' # img= none,png,eps,pdf #domain_list=['PERU','CALIFORNIA','SENEGAL','BENGUELA'] #domain_list=['PERU','CALIFORNIA'] domain_list=['PERU'] #domain_list=['CALIFORNIA'] #domain_list=['SENEGAL'] #domain_list=['BENGUELA'] #depth_list=['1','54','87','29_114','120','181'] #depth_list=['1','54','87','120','181'] #depth_list=['1','29_114','120'] #depth_list=['1','54','120'] #depth_list=['1','87'] #depth_list=['1','120'] #depth_list=['1'] depth_list=['1','29_114','181'] #depth_list=['29_114'] # 85W n'a pas de sens car la detection de l'EUC ne fonctionne pas => valeurs de temperature absurdes # section_list_pacific=['98.5W','85W','110.5W','120.25W','130W','140.5W','170.5W'] # section_list_pacific=['85W'] #section_list_pacific=['98.5W','110.5W','120.25W','130W','140.5W','170.5W'] #section_list_pacific=['98.5W','110.5W','120.25W','170.5W'] #section_list_pacific=['98.5W','110.5W','120.25W'] #section_list_pacific=['98.5W','110.5W'] #section_list_pacific=['98.5W'] #section_list_pacific=['110.5W'] #section_list_pacific=['120.25W'] #section_list_pacific=['130W'] #section_list_pacific=['140.5W'] section_list_pacific=['170.5W'] # 10W n'a pas de sens car la detection de l'EUC ne fonctionne pas => valeurs de temperature absurdes # section_list_atlantic=['10W','20_5W','30_25W'] # section_list_atlantic=['10W'] #section_list_atlantic=['20_5W','30_25W'] #section_list_atlantic=['30_25W'] section_list_atlantic=['20_5W'] area_name_idem=True #area_name_list_idem=['T_upA','T_upB','T_upBas'] #area_name_list_idem=['T_upA'] #area_name_list_idem=['T_upB'] area_name_list_idem=['T_upBas'] #seuil_uo='EUC_uo_max' seuil_uo='EUC_uo_seuil_20pct' #seuil_uo='30%_of_uo' # LEGENDS: color_dots=True # False: meme couleur pour toutes les simus # True : une couleur par type (erai, quick, cmip5...) identifier=False # True: plot les identifiants de chaque point (chaque simu), # => genere un fichier texte avec les legendes plot_legend_simu=False # plot les legendes de chaque simu sur le plot # bien quand il n'y a que quelques simulations. label_plus=True # pour chaque droite de regression, indique la pente et le coefficiant linear_regression_line_label=True size2='14' size3='18' #xlim1=14; xlim2=23 # pour tous les plots! xlim1=13; xlim2=23 #ylim1=5; ylim2=30 # pour tous les plots! ylim1=5; ylim2=28 # pour tous les plots! #ylim1=10; ylim2=28 # Peru #ylim1=5; ylim2=23 # California #ylim1=10; ylim2=30 # Benguela #ylim1=9; ylim2=30 # Senegal if domain_list[0] == 'PERU': ylim1=10; ylim2=29 elif domain_list[0] == 'CALIFORNIA': ylim1=7; ylim2=25 elif domain_list[0] == 'SENEGAL': ylim1=10; ylim2=29 elif domain_list[0] == 'BENGUELA': ylim1=8; ylim2=30 yplus=ylim2-1 #yplus=24 # Peru T_upA #yplus=27 # Peru T_upB et T_upBas #figsize1=(15,15) figsize1=(xlim2-xlim1,ylim2-ylim1) xfigsize1=xlim2-xlim1 yfigsize1=ylim2-ylim1 ms1=5 #import t2 y_obs1=''; y_obs2='' id_obs1=''; id_obs2='' #------------------------------------------------------------------------------- # python Modules import #------------------------------------------------------------------------------- import matplotlib as mpl print_image=True if img=='png' : mpl.use('Agg') elif img=='eps' : mpl.use('ps') elif img=='pdf' : mpl.use('pdf') else : print_image=False import matplotlib.pyplot as plt from netCDF4 import Dataset as netcdf import numpy as np import numpy.ma as ma import sys import os from matplotlib.font_manager import FontProperties import datetime import shutil import matplotlib.offsetbox as offsetbox local_time1=datetime.datetime.strftime(datetime.datetime.now(), '%Y-%m-%d') # '2018-09-19' local_time2=datetime.datetime.strftime(datetime.datetime.now(), '%d-%b-%Y') # '19-Sep-2018' fig1, ax1 = plt.subplots() fig1.set_figwidth(xfigsize1) fig1.set_figheight(yfigsize1) ax2 = ax1.twinx() domain_text_incr=0 area_label='' domain_str='' kk=0 label1=[] label2=[] for domain in domain_list : #------------------------------------------------------------------------------------ # read data #------------------------------------------------------------------------------------ # PERU (EUC_uo_seuil_20pct) data_file='list_'+domain+'.txt' print 'reading data file '+data_file+'...' f1 = open(data_file, 'r') # list_PERU.txt in one line # str1=f1.read() # exec(str1) # list_PERU.txt with a line for a variable for line in f1.readlines(): exec(line.strip()) # line without '\n' with .strip() f1.close() print 'cmip5_res=%3s , cmip5_2D_interp_method=%s , F_hfds=%s'%(cmip5_res,cmip5_2D_interp_method,F_hfds) label1.append('%s : '%domain) label2.append(domain+'_') domain_str=domain+'_'+domain_str print '\n',domain,'\n' #------------------------------------------------------------------------------------ # Sections ATLANTIC / PACIFIC #------------------------------------------------------------------------------------ if (domain == 'PERU') | (domain == 'CALIFORNIA') : section_list=section_list_pacific pacific_str1='section(' pacific_str2='' for section in section_list : pacific_str1=pacific_str1+section+' ' pacific_str2=pacific_str2+section+'_' else : section_list=section_list_atlantic atlantic_str1='section(' atlantic_str2='' for section in section_list : atlantic_str1=atlantic_str1+section+' ' atlantic_str2=atlantic_str2+section+'_' if area_name_idem == True : #------------------------------------------------------------------------------------ # choix area par domain #------------------------------------------------------------------------------------ area_name_list=area_name_list_idem #area_name_list=['T_upBas'] else : if domain == 'PERU' : area_name_list=['T_upA','T_upBas'] elif domain == 'CALIFORNIA' : area_name_list=['T_upA'] elif domain == 'SENEGAL' : area_name_list=['T_upBas'] elif domain == 'BENGUELA' : area_name_list=['T_upBas'] #------------------------------------------------------------------------------------ if seuil_uo == 'EUC_uo_max' : #------------------------------------------------------------------------------------ if domain == 'PERU' : # PERU (EUC_uo_max) exp_list=['tr075n_newflx','tr075n_n02nogm','tr075_quik','tr075n_quik','tr075_erainocrt','tr075n_erainocrt','tr075_cpl12L60_sol093','tr075_cpl12L60_sol094','tr075now_11L60m_swrad04','tr075now_12L60_sol094obcokalb3','tr075now_44L60_sol099obcokalb3zhao3nn','tr025_cpl12L60_sol094','tr025n_gm00','tr025n_gm01','tr025_quik','tr025_quikubs','tr025_quiktau','tr025_erainocrt','tr12_quik','tr12_erainocrt','tr12n_n02','tr12n_n03'] elif domain == 'CALIFORNIA' : # CALIFORNIA (EUC_uo_max) exp_list=['tr075n_newflx','tr075n_n02nogm','tr075_quik','tr075n_quik','tr075_erainocrt','tr075n_erainocrt','tr075_cpl12L60_sol093','tr075_cpl12L60_sol094','tr075now_11L60m_swrad04','tr075now_12L60_sol094obcokalb3','tr075now_44L60_sol099obcokalb3zhao3nn','tr025_cpl12L60_sol094','tr025n_gm00','tr025n_gm01','tr025_quik','tr025_quikubs','tr025_quiktau','tr025_erainocrt','tr12_quik','tr12_erainocrt','tr12n_n02','tr12n_n03'] elif domain == 'SENEGAL' : # SENEGAL (EUC_uo_max) exp_list=['tr075n_newflx','tr075n_n02nogm','tr075_quik','tr075n_quik','tr075_erainocrt','tr075n_erainocrt','tr075_cpl12L60_sol093','tr075_cpl12L60_sol094','tr075now_11L60m_swrad04','tr075now_12L60_sol094obcokalb3','tr075now_44L60_sol099obcokalb3zhao3nn','tr025_cpl12L60_sol094','tr025n_gm00','tr025n_gm01','tr025_quik','tr025_quikubs','tr025_quiktau','tr025_erainocrt','tr12_quik','tr12_erainocrt','tr12n_n02','tr12n_n03'] elif domain == 'BENGUELA' : # BENGUELA (EUC_uo_max) exp_list=['tr075n_newflx','tr075n_n02nogm','tr075_quik','tr075n_quik','tr075_erainocrt','tr075n_erainocrt','tr075_cpl12L60_sol093','tr075_cpl12L60_sol094','tr075now_11L60m_swrad04','tr075now_12L60_sol094obcokalb3','tr075now_44L60_sol099obcokalb3zhao3nn','tr025_cpl12L60_sol094','tr025n_gm00','tr025n_gm01','tr025_quik','tr025_quikubs','tr025_quiktau','tr025_erainocrt','tr12_quik','tr12_erainocrt','tr12n_n02','tr12n_n03'] #=============================================================================== elif seuil_uo == 'EUC_uo_seuil_20pct' : #=============================================================================== if domain == 'PERU' : domain_color='blue' elif domain == 'CALIFORNIA' : domain_color='magenta' elif domain == 'SENEGAL' : domain_color='red' elif domain == 'BENGUELA' : domain_color='green' #------------------------------------------------------------------------------- area_name_str='' area_text_incr=0 label1[kk]=label1[kk]+'area(' for area_name in area_name_list : area_name_str=area_name+'_'+area_name_str if area_name=='T_upA': area_title='UpA area' area_style=':' # dotted segment_size=1.0 area_style2='....' elif area_name=='T_upB': area_title='UpB area' area_style='-' # solid area_style2='____' segment_size=0.5 elif area_name=='T_upBas': area_title='upB area (Along Shore) ' area_style='-.' #dashdot area_style2='-.-.-.' segment_size=0.7 label1[kk]=label1[kk]+area_style2+' '+area_name.replace('T_','')+' ' label2[kk]=label2[kk]+area_name.replace('T_','')+'_' # depth_str1='depth( ' depth_str2='' label_plus_incr=0 if F_hfds : HFDS_name=area_name.replace('T_','HFDS_') HFDS_list=eval(HFDS_name+'_list') decompte=len(depth_list) for depth in depth_list : # print 'depth=', depth depth_str1=depth_str1+depth+'m ' depth_str2=depth_str2+depth+'m_' # T_area_str=area_name+'_'+depth+'m' # print T_area_str T_area_list=eval(T_area_str+'_list') mm=0 label_plus_incr=label_plus_incr+0.3 for jj in range(len(section_list)) : section=section_list[jj] T_EUC_list=eval('S_EUC_'+section.replace('.','_')+'_list') # print 'T_EUC_list=', T_EUC_list # print 'section_list=', section_list # print 'section=', section print jj,T_EUC_list legend1='' T_area_mm=np.zeros(len(T_upBas_120m_list)) T_EUC_mm=np.zeros(len(T_upBas_120m_list)) for ii in range(len(T_upBas_120m_list)) : T_area=T_area_list[ii] if F_hfds : HFDS=HFDS_list[ii] T_EUC=T_EUC_list[ii] exp=exp_list[ii] id=id_list[ii] if F_hfds : print ii,' ',id,' ',exp,' ',T_EUC,' ',T_area,' ',HFDS else : print ii,' ',id,' ',exp,' ',T_EUC,' ',T_area # if id[0:2] == '75' : marker1='x' elif id[0:2] == '25' : marker1='4' elif id[0:2] == '13' : marker1='+' # trop 12 elif id[0:2] == '12' : marker1='o' # BDY 12 elif id[0:2] == 'c5' : marker1='^' # cmip5 elif id[0:2] == 'gl' : marker1='x' # glorys # if (id=='12q1') : T_area_12q1=T_area T_EUC_12q1=T_EUC color1='cyan' if (id=='12q2') : T_area_12q2=T_area T_EUC_12q2=T_EUC color1='cyan' if (id=='12q3') : T_area_12q3=T_area T_EUC_12q3=T_EUC color1='yellow' if (id=='12q4') : T_area_12q4=T_area T_EUC_12q4=T_EUC color1='yellow' # if color_dots == True : ms1=5 if id[0:2] == 'c5' : color1='magenta' elif id == 'gl01' : ii_obs1=ii id_obs1=id T_area_obs1=T_area T_EUC_obs1=T_EUC color1_obs1='magenta' color1=color1_obs1 marker1='^' ms1=10 elif id == 'gl02' : ii_obs2=ii id_obs2=id T_area_obs2=T_area T_EUC_obs2=T_EUC color1_obs2='green' color1=color1_obs2 marker1='^' ms1=10 else : if id[2] == 'e' : # color1='red' color1='green' ms1=10 elif id[2] == 'q' : # color1='blue' color1='green' ms1=10 elif id[2] == 'c' : # color1='green' color1='magenta' ms1=10 color2='k' if ('GISS' in exp ) : color1='black' ms1=10 else : color1=domain_color color2=domain_color # legend1=legend1+id+'='+'%-40s'%(exp)+'\n' T_area_mm[ii]=T_area T_EUC_mm[ii]=T_EUC ax1.plot(T_EUC,T_area,marker1,color=color1, markeredgewidth=1, markersize=ms1) if F_hfds : ax2.plot(T_EUC,HFDS,marker1,color='green', markeredgewidth=1, markersize=ms1) if identifier : if 'p12n_bdy' in exp : ax2.text(T_EUC,HFDS,id,color='green') if identifier : if 'p12n_bdy' in exp : ax1.text(T_EUC,T_area,id) xmin=T_EUC_mm.min() xmax=T_EUC_mm.max() imin=np.where(T_EUC_mm==T_EUC_mm.min())[0][0] imax=np.where(T_EUC_mm==T_EUC_mm.max())[0][0] ymin=T_area_mm.min() ymax=T_area_mm.max() # Droite de regression from sklearn import linear_model from sklearn.metrics import r2_score regr1 = linear_model.LinearRegression() regr2 = linear_model.LinearRegression() # dim=len(id_list) xx=np.zeros(dim) yy=np.zeros(dim) hf=np.zeros(dim) yy_pred=np.zeros(dim) hf_pred=np.zeros(dim) for ii in range(dim) : # print ii,T_EUC_list[ii] xx[ii]=T_EUC_list[ii] for ii in range(dim) : # print ii,T_EUC_list[ii] yy[ii]=T_area_list[ii] if F_hfds : hf[ii]=HFDS_list[ii] if (id_obs1 == 'gl01') : y_obs1=yy[ii] # print 'y_obs1=',y_obs1 if (id_obs2 == 'gl02') : y_obs2=yy[ii] # print 'y_obs2=',y_obs2 regr1.fit(xx[:,np.newaxis],yy) if F_hfds : regr2.fit(xx[:,np.newaxis],hf) for ii in range(dim) : yy_pred[ii]=regr1.predict(xx[ii]) if F_hfds : hf_pred[ii]=regr2.predict(xx[ii]) ax1.plot([xx[ii],xx[ii]], [yy[ii],yy_pred[ii]], color=domain_color, linewidth=segment_size, linestyle=area_style) # ax1.plot(xx[ii],yy_pred[ii],marker1,color=color1, mew=2, ms=1) if F_hfds : ax2.plot([xx[ii],xx[ii]], [hf[ii],hf_pred[ii]], color='green', linewidth=segment_size, linestyle=area_style) #x_test = np.linspace(np.min(xx), np.max(xx), 100) xx_test = np.linspace(xmin, xmax, 100) if area_name_list_idem==['T_upA','T_upB','T_upBas']: ax1.plot(xx_test, regr1.predict(xx_test[:,np.newaxis]), color=domain_color, linewidth=2, linestyle=area_style) ax2.plot(xx_test, regr2.predict(xx_test[:,np.newaxis]), color='green', linewidth=2, linestyle=area_style) else : ax1.plot(xx_test, regr1.predict(xx_test[:,np.newaxis]), color=domain_color, linewidth=2, linestyle='-') if F_hfds : ax2.plot(xx_test, regr2.predict(xx_test[:,np.newaxis]), color='green', linewidth=2, linestyle='-') if (id_obs1 == 'gl01') : print 'y_obs1=',y_obs1 y_glorys=np.full((xx_test.shape), T_area_obs1) ax1.plot(xx_test, y_glorys, color=color1_obs1, linewidth=1, linestyle='-') if (id_obs2 == 'gl02') : print 'y_obs2=',y_obs2 y_glorys=np.full((xx_test.shape), T_area_obs2) ax1.plot(xx_test, y_glorys, color=color1_obs2, linewidth=1, linestyle='-') # ax1.plot(xx, yy, color=domain_color, linewidth=1, linestyle=area_style) # ax1.plot(xx, yy_pred, color=domain_color, linewidth=1, linestyle=area_style) r2=r2_score(yy,yy_pred) if F_hfds : r2_hfds=r2_score(hf,hf_pred) label3='%sm %s'%(depth,section) label4='%sm %s (a=%.2f r2=%.2f)'%(depth,section,regr1.coef_[0],r2) if F_hfds : label5='HFDS %s (a=%.2f r2=%.2f)'%(section,regr2.coef_[0],r2_hfds) # print 'label4=',label4 if linear_regression_line_label : if mm%2 == 0 : ax1.text(xx[imax],yy_pred[imax],label3,verticalalignment='bottom',horizontalalignment='left',color=color2) else : ax1.text(xx[imin],yy_pred[imin],label3,verticalalignment='top',horizontalalignment='right',color=color2) mm=mm+1 if label_plus : ax1.text(xlim2,yplus-label_plus_incr,label4,verticalalignment='bottom',horizontalalignment='right',color=color2, size=size2) if decompte == 1 : if F_hfds : ax1.text(xlim2,yplus-label_plus_incr*2,label5,verticalalignment='bottom',horizontalalignment='right',color='green', size=size2) label_plus_incr=label_plus_incr + 0.3 del section # # print ii, "xmax, ymax",xmax,ymax # ax1.text(xmax,ymax,area_name+'('+depth+'m)',verticalalignment='top',horizontalalignment='left',color=domain_color) if plot_legend_simu : # legend des simulations ax1.text(xmin,ymax,legend1,verticalalignment='top',horizontalalignment='left',multialignment='left') # del depth decompte=decompte-1 print 'decompte ', decompte # ax1.text(xlim2,ylim2-1-area_text_incr,area_style2+' '+area_title,verticalalignment='top',horizontalalignment='right',color='k') area_text_incr=area_text_incr + 0.5 # if (domain == 'PERU') | (domain == 'CALIFORNIA') : label1[kk]=label1[kk]+') '+pacific_str1+')' label2[kk]=label2[kk] else : label1[kk]=label1[kk]+') '+atlantic_str1+')' label2[kk]=label2[kk] # label1[kk]=project+' '+label1[kk] print 'label1 (title):',label1[kk] ax1.text(xlim1,ylim2-domain_text_incr,label1[kk],verticalalignment='top',horizontalalignment='left',color=domain_color,size=size3) # domain_text_incr=domain_text_incr + 0.3 kk=kk+1 # label2_str='' for ll in range(kk) : label2_str=label2_str+label2[ll] if (domain == 'PERU') | (domain == 'CALIFORNIA') : label2_str=label2_str+'_'+depth_str2+'_'+pacific_str2 if (domain == 'SENEGAL') | (domain == 'BENGUELA') : label2_str=label2_str+'_'+depth_str2+'_'+atlantic_str2 #label2_str=label2_str+'_'+seuil_uo+'_'+project label2_str=label2_str+'_'+seuil_uo print 'label2_str (figname):',label2_str ax1.grid() ax1.set_xlabel('EUC temperature ') ax1.set_ylabel('Upwelling area Temperature') ax2.set_ylabel('Surface flux (W/m2)') ax1.set_xticks(np.arange(xlim1, xlim2, 1)) ax1.set_yticks(np.arange(ylim1, ylim2, 1)) if F_hfds : title1=seuil_uo+' (cmip5 '+cmip5_2D_interp_method+' '+cmip5_res+') - '+F_Flux else : title1=seuil_uo+' (cmip5 '+cmip5_2D_interp_method+' '+cmip5_res+')' ax1.text(xlim1,ylim2-domain_text_incr,' '+depth_str1+')',verticalalignment='top',horizontalalignment='left',color='k') #ax1.text(xlim2,ylim2,depth_str1,verticalalignment='top',horizontalalignment='right',color='k') ax1.set_title(title1 , fontsize=18,fontweight='bold' ) ax1.set_xlim(xlim1,xlim2) ax1.set_ylim(ylim1,ylim2) ax2.set_ylim(-600,150) #ax2.set_ylim(0,600) ax1.text(xlim2,ylim1,local_time2,color='g',verticalalignment='bottom',horizontalalignment='right') plt.show() #domain_str='' #for domain in domain_list : # domain_str=domain+'_'+domain_str #area_name_str='' #for area_name in area_name_list : # area_name_str=area_name+'_'+area_name_str #depth_str1='' #for depth in depth_list : # depth_str1=depth+'m_'+depth_str1 #plt.savefig(domain_str+area_name_str+depth_str2+seuil_uo+'.'+img) #plt.savefig(label2_str+'_'+cmip5_2D_interp_method+'_'+cmip5_res+'.'+img) if F_hfds : label3_str=label2_str+'_'+F_Flux else : label3_str=label2_str+'_'+section_list[0] if identifier : plt.savefig(label3_str+'_identifier.'+img) flegend = open(label3_str+'_legend.txt', 'w') flegend.write(legend1) flegend.close() else : plt.savefig(label3_str+'.'+img) scriptname=sys.argv[0] shutil.copy2(scriptname,label3_str+'.py') plt.close(fig1)