source: BOL/Multi_atlas/atlas/atlas_test_NC.py @ 3817

Last change on this file since 3817 was 3684, checked in by idelkadi, 5 years ago

Repository under svn of a first version of Multiatlas diagnostics for LMDZ. This version is adapted to be able to run a LMDZ multiatlas on an individual account on the ciclad machine of the IPSL. In this version, the parts to be modified are identified so as to subsequently adapt it to other machines.
This version is still under development.

File size: 11.2 KB
Line 
1# -*- coding: iso-8859-1 -*-
2# Created : S.Sénési - nov 2015
3
4#
5desc="\nCreation d'un atlas pour une simu, une grille et une liste de variables et de saisons \n"+\
6"  Exemples : \n"+\
7"  >>> python ./atlas.py -v tas,hfls -s NPv3.1ada_1982_1991\n"+\
8""
9# Avec CliMAF, cette etape est loin d'etre necessaire; on la réalise pour 'exposer' ces fichiers
10# dans une arborescence à laquelle sont habitues certains utilisateurs
11
12
13# Répertoire de base pour les entrées et les résultats
14#dir_default='/data/hourdin/LMDZ6/SE/ORIG'
15#dir_default='/prodigfs/fabric/LMDZ6/SE/ORIG'
16dir_default='/prodigfs/ipslfs/dods/fabric/lmdz/SE/ORIG'
17
18# Gestion des options et arguments d'appel
19from optparse import OptionParser
20parser = OptionParser(desc) ; parser.set_usage("%%prog [-h]\n%s" % desc)
21#parser.add_option("-i", "--input", help="repertoire des donnes d'entree(defaut : %s)"%dir_default,
22#                  action="store",default=dir_default)
23parser.add_option("-i", "--input", help="repertoire des donnes d'entree(defaut : %s)",
24                  action="store",default=None)
25parser.add_option("-g", "--grid", help="nom de grille (default: VLR)", action="store",default='VLR')
26parser.add_option("-r", "--region", help="nom de zone (default: GLOB)", action="store",default='GLOB')
27parser.add_option("-p", "--season", help="saison a traiter " "(eg : JJA, DJF, YEAR, defaut=%YEAR)", 
28                  action="store", default='YEAR')
29parser.add_option("-s", "--simulation", help="simulation+annees a traiter (sim_YYY1_YYY2) - laisser vide pour lister",
30                  action="store",default=None)
31parser.add_option("-t", "--reference", help="simulation de reference (sim_YYY1_YYY2, default=OBS) ",
32                  action="store",default='OBS')
33parser.add_option("-v", "--variables", help="liste des variables (separees par des virgules)", action="store",default=None)
34#parser.add_option("--root", help="Path to the root directory", action="store",default=None)
35parser.add_option("--root", help="Path to the root directory", action="store",default=None)
36parser.add_option("-f", "--force", help="force le recalcul de champs existants", 
37                  action="store_true",default=None)
38parser.add_option("-o", "--pdf", help="nom du pdf de sortie (default: atlas_<SIMU>_<SAISON>.pdf)", action="store")
39(opts, args) = parser.parse_args()
40
41#---------------------------------------------------------------------------------
42import math
43from climaf.api import *
44from climaf.html import * 
45# La description de l'organisation des données SE et des alias et rescalings
46# est partagée dans une micro-librairie :
47from lmdz_SE import * # svsg, all_SE_simulations
48from plot_params import plot_params
49#---------------------------------------------------------------------------------
50#
51def apply_scale_offset(dat,scale,offset):
52    return ccdo(ccdo(dat,operator='mulc,'+str(float(scale))),operator='addc,'+str(float(offset)))
53#
54#craz()
55if opts.simulation is None:
56    print "Available simulations at %s are : "%opts.input,
57    for s in all_SE_simulations() : print s,
58    exit(0)
59#
60lvars=opts.variables
61if lvars is not None : lvars=lvars.split(',')
62else : lvars=variables_list
63#
64# Preparons une commande pour assembler les sorties Pdf
65if opts.pdf : pdffile=opts.pdf
66else: pdffile="atlas_"+opts.simulation+"_"+opts.season+".pdf"
67pdfargs=["pdfjam","--landscape","-o ",pdffile]
68#
69# Initialisation de l'index html
70index= header("LMDZ Atlas for "+opts.simulation+ " versus "+opts.reference+" ("+opts.season+")") 
71index += cell('PDF',pdffile)
72index += section("2d vars", level=4)
73index += open_table()
74#
75# Titres de colonnes
76ref=opts.reference ; 
77if (ref == 'OBS' ) : text_diff='bias'
78else:                text_diff='diff'
79index+=open_line('VARIABLE')+cell('bias')+cell('rmse')+cell('mean')+cell(ref)+cell(text_diff)+\
80        cell('zonal')+cell('all')+cell('pdf')+close_line()
81#
82# -- Declare the script ml2pl for vertical interpolation
83cscript("ml2pl", "/home/jservon/Evaluation/CliMAF/Atlas_LMDz/ml2pl.sh -p ${var_2} -m ${var_1} ${in_1} ${out} ${in_2}",
84    commuteWithTimeConcatenation=True, commuteWithSpaceConcatenation=True)
85# -- Vertical levels for the vertical interpolation
86fixed_fields("ml2pl",("press_levels.txt","/home/fabric/LMDZ/atlas/press_levels.txt"))
87#
88for variable  in lvars :
89    # Get the model and the reference
90    if opts.root:
91       simu=svsg(opts.simulation,variable,opts.season,opts.grid, root=opts.root)
92    else:
93       simu=svsg(opts.simulation,variable,opts.season,opts.grid)
94    print 'variable = ',variable
95    reff=svsg(opts.reference,variable,opts.season,opts.grid)
96    #
97    # If the variable is a 3D field:
98    #  - interpolate the variable on the standard pressure levels with ml2pl (L. Guez)
99    #  - Compute the difference model-ref with diff_zonmean (computes the zonal mean lat/pressure fields,
100    #    interpolates the model on the ref, both vertically and horizontally, and returns the difference)
101    if is3d(variable) :
102       simu_pres = svsg(opts.simulation,'pres',opts.season,opts.grid)
103       simu = ml2pl(simu,simu_pres)
104       simu = zonmean(simu)
105       reff = zonmean(reff)
106       diff = diff_zonmean(simu,reff)
107    else:
108        if (opts.grid == '' ) : reff=regrid(reff,simu)
109        diff=minus(simu,reff)
110
111    pparams = plot_params(variable,'full_field')
112    vertical_interval = 'trYMaxF=1000|trYMinF=1'
113    stringFontHeight=0.018
114    if is3d(variable):
115        pparams.update({'options':vertical_interval})
116        stringFontHeight=0.03
117    # Map for simulation
118    simu_fig=plot(simu,title="",
119                  gsnLeftString=variable,
120                  gsnCenterString=opts.simulation,
121                  gsnRightString=opts.season,
122                  gsnStringFontHeightF=stringFontHeight,
123                  mpCenterLonF=0,
124                  **pparams)
125    simu_avg=cvalue(space_average(simu))
126    #
127    # Map for reference
128    ref_fig=plot(reff,title="",
129                 gsnLeftString=variable,
130                 gsnCenterString=ref,
131                 gsnRightString=opts.season,
132                 gsnStringFontHeightF=stringFontHeight,
133                 mpCenterLonF=0,
134                 **pparams)
135    ref_avg=cvalue(space_average(reff))
136    #
137    # Bias (or difference between simulations) map
138    if (ref == 'OBS' ) : p=plot_params(variable,'bias')
139    else:                p=plot_params(variable,'model_model')
140    tmp_aux_params = plot_params(variable,'full_field')
141    scale = 1.0 ; offset = 0.0
142    if 'offset' in tmp_aux_params or 'scale' in tmp_aux_params:
143       if 'offset' in tmp_aux_params:
144          offset = tmp_aux_params['offset']
145       else:
146          offset=0.0
147       if 'scale' in tmp_aux_params:
148          scale = tmp_aux_params['scale']
149       else:
150          scale=1.0
151       wreff = apply_scale_offset(reff,scale,offset)
152       wsimu = apply_scale_offset(simu,scale,offset)
153    else:
154       wreff = reff
155       wsimu = simu
156    #
157    if is3d(variable):
158        p.update({'options':vertical_interval})
159    if variable in ['ua','va','ta','hus','hur']:
160        tmp_levs = tmp_aux_params['colors']
161        p.update({'contours':tmp_levs})
162        diff_fig=plot(diff,wreff,title="", format='png', mpCenterLonF=0,
163                  gsnLeftString=variable,
164                  gsnCenterString=opts.simulation+' - '+ref,
165                  gsnRightString=opts.season,
166                  gsnStringFontHeightF=stringFontHeight,
167                  aux_options='cnLineThicknessF=2|cnLineLabelsOn=True', **p)
168    else:
169        p.update({'contours':1})
170        diff_fig=plot(diff,title="", format='png', mpCenterLonF=0,
171                  gsnLeftString=variable,
172                  gsnCenterString=opts.simulation+' - '+ref,
173                  gsnRightString=opts.season,
174                  gsnStringFontHeightF=stringFontHeight,
175                  **p)
176
177    #
178    # Bias mean value, and RMSD/RMSE
179    diff_avg=cvalue(space_average(diff))
180    rmsd=math.sqrt(cvalue(space_average(ccdo(diff,operator='-b F64 sqr'))))
181    #
182    # Zonal means
183    if not is3d(variable):
184        # -- apply a mask corresponding to the reference
185        mask = fdiv(reff,reff)
186        msimu = fmul(wsimu,mask)
187        # -- Compute the zonal mean
188        zmean=ccdo(msimu, operator='zonmean')
189        ref_zmean=ccdo(wreff, operator='zonmean')
190        #
191        sim=opts.simulation
192        #if variable in ['zg500']:
193        #   ref_zmean = ccdo(ref_zmean,operator='-b F32 mulc,1')
194        #   zmean = ccdo(zmean,operator='-b F32 mulc,1')
195        #zmean_fig=curves(cens([sim,ref],zmean,ref_zmean),
196        zmean_fig=curves(cens({sim:zmean,ref:ref_zmean}),
197                         title="",
198                         lgcols=3,
199                         options=#'tiYAxisString=""|'+\
200                                 #'+\'+\
201                                 'tmYROn=True|'+\
202                                 'tmYRBorderOn=True|'+\
203                                 'tmYLOn=False|'+\
204                                 'tmYUseRight=True|'+\
205                                 'vpXF=0|'+\
206                                 'vpWidthF=0.66|'+\
207                                 'vpHeightF=0.33|'+\
208                                 'tmYRLabelsOn=True|'+\
209                                 'tmXBLabelFontHeightF=0.018|'+\
210                                 'tmYLLabelFontHeightF=0.016|'+\
211                                 'lgLabelFontHeightF=0.018|'+\
212                                 #'pmLegendSide=Bottom|'+\
213                                 'pmLegendOrthogonalPosF=-0.32|'+\
214                                 'pmLegendParallelPosF=1.0|'+\
215                                 'tmXMajorGrid=True|'+\
216                                 'tmYMajorGrid=True|'+\
217                                 'tmXMajorGridLineDashPattern=2|'+\
218                                 'tmYMajorGridLineDashPattern=2|'+\
219                                 'xyLineThicknessF=8|'+\
220                                 'gsnLeftString='+variable+'|'+\
221                                 'gsnCenterString='+opts.simulation+' vs '+ref+'|'+\
222                                 'gsnRightString='+opts.season+'|'+\
223                                 'gsnStringFontHeightF='+str(stringFontHeight))
224    #    # Composite figure
225    if is3d(variable):
226        page=cpage([[simu_fig,ref_fig],[diff_fig,None]],orientation='landscape', page_trim=True, fig_trim=True)
227        pdf_page=cpage([[simu_fig,ref_fig],[diff_fig,None]],orientation='landscape',
228                   page_trim=True, fig_trim=True, format='pdf', title=variable+" "+opts.simulation)
229    else:
230        page=cpage([[simu_fig,ref_fig],[diff_fig,zmean_fig]],orientation='landscape', page_trim=True, fig_trim=True)
231        pdf_page=cpage([[simu_fig,ref_fig],[diff_fig,zmean_fig]],orientation='landscape', 
232                   page_trim=True, fig_trim=True, format='pdf', title=variable+" "+opts.simulation)
233    pdfargs.append(cfile(pdf_page))
234    #
235    thumbnail_size = 200
236    if is3d(variable):
237            index+=open_line(varlongname(variable)+' ('+variable+')')+\
238                    cell("%.2g"%diff_avg,cfile(diff_fig))+\
239                    cell("%.2g"%rmsd,cfile(diff_fig))+\
240                    cell(simu,cfile(simu_fig),thumbnail=thumbnail_size,hover=False)+\
241                    cell(ref,cfile(ref_fig),thumbnail=thumbnail_size,hover=False)+\
242                    cell(text_diff,cfile(diff_fig),thumbnail=thumbnail_size,hover=False)+\
243                    ' '+\
244                    cell('page',cfile(page),thumbnail=thumbnail_size,hover=False)+\
245                    cell('Pdf',cfile(pdf_page))
246            close_line()
247    else:
248            index+=open_line(varlongname(variable)+' ('+variable+')')+\
249                    cell("%.2g"%diff_avg,cfile(diff_fig))+\
250                    cell("%.2g"%rmsd,cfile(diff_fig))+\
251                    cell(simu,cfile(simu_fig),thumbnail=thumbnail_size,hover=False)+\
252                   cell(ref,cfile(ref_fig),thumbnail=thumbnail_size,hover=False)+\
253                    cell(text_diff,cfile(diff_fig),thumbnail=thumbnail_size,hover=False)+\
254                    cell('zonal mean',cfile(zmean_fig),thumbnail=thumbnail_size,hover=False)+\
255                    cell('page',cfile(page),thumbnail=thumbnail_size,hover=False)+\
256                    cell('Pdf',cfile(pdf_page))
257            close_line()
258#
259# Finalisons l'index html
260index += close_table()
261index += trailer()
262#out="index_example.html"
263out="index_example_"+opts.season+"_"+opts.simulation+".html"
264with open(out,"w") as filout : filout.write(index)
265#
266# Creation du Pdf multi-pages
267comm=subprocess.Popen(pdfargs, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
268#
269import os,os.path ; 
270# print("Attendez un bon peu : lancemement de firefox sur Ciclad....")
271# os.system("firefox file://"+os.path.abspath(os.path.curdir)+"/"+out+"&")
Note: See TracBrowser for help on using the repository browser.