source: BOL/Multi_atlas/atlas/atlas_none_1.0.3.py @ 4995

Last change on this file since 4995 was 3684, checked in by idelkadi, 4 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: 12.9 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
49crm(pattern='pme')
50crm(pattern='ua850')
51#---------------------------------------------------------------------------------
52#
53def apply_scale_offset(dat,scale,offset):
54    return ccdo(ccdo(dat,operator='mulc,'+str(float(scale))),operator='addc,'+str(float(offset)))
55#
56#craz()
57if opts.simulation is None:
58    print "Available simulations at %s are : "%opts.input,
59    for s in all_SE_simulations() : print s,
60    exit(0)
61#
62lvars=opts.variables
63if lvars is not None : lvars=lvars.split(',')
64else : lvars=variables_list
65#
66# Preparons une commande pour assembler les sorties Pdf
67if opts.pdf : pdffile=opts.pdf
68else: pdffile="atlas_"+opts.simulation+"_"+opts.season+".pdf"
69pdfargs=["pdfjam","--landscape","-o ",pdffile]
70#
71# Initialisation de l'index html
72index= header("LMDZ Atlas for "+opts.simulation+ " versus "+opts.reference+" ("+opts.season+")") 
73index += cell('PDF',pdffile)
74index += section("2d vars", level=4)
75index += open_table()
76#
77# Titres de colonnes
78ref=opts.reference ; 
79if (ref == 'OBS' ) : text_diff='bias'
80else:                text_diff='diff'
81index+=open_line('VARIABLE')+cell('bias')+cell('rmse')+cell('mean')+cell(ref)+cell(text_diff)+\
82        cell('zonal')+cell('all')+cell('pdf')+close_line()
83#
84# -- Declare the script ml2pl for vertical interpolation
85cscript("ml2pl", "/home/jservon/Evaluation/CliMAF/Atlas_LMDz/ml2pl.sh -p ${var_2} -m ${var_1} ${in_1} ${out} ${in_2}",
86    commuteWithTimeConcatenation=True, commuteWithSpaceConcatenation=True)
87# -- Vertical levels for the vertical interpolation
88fixed_fields("ml2pl",("press_levels.txt","/home/fabric/LMDZ/atlas/press_levels.txt"))
89#
90for variable  in lvars :
91    #
92    print 'variable = ',variable
93    # -- Interpolation sur les niveaux verticaux a partir des fichiers interpoles avec ml2pl
94    dum3dvars = ['ua', 'va', 'ta', 'hur', 'hus', 'zg']
95    # --> On cherche un pattern du genre ua200, ta850 (et on veut éviter hurs, tas...)
96    var_interp = False
97    for dum3dvar in dum3dvars:
98        # -- si on le trouve...
99        if dum3dvar in variable:
100           # -- On enleve le nom de la variable 3d du nom de variable
101           tmplev = str.replace(variable,dum3dvar,'')
102           # -- Et si ce qu'il reste n'est pas 's' (comme dans tas ou hurs)
103           if tmplev not in ['s', '']:
104              var_interp = dum3dvar
105    # -- Si variable est une variable qui contient un nom de variable 3d avec un niveau vertical:
106    if var_interp:
107       if opts.root:
108          simu=svsg(opts.simulation,var_interp,opts.season,opts.grid, root=opts.root)
109       else:
110          simu=svsg(opts.simulation,var_interp,opts.season,opts.grid)
111       # -- On peut selectionner le niveau vertical dans le fichier interpole
112       simu_pres = svsg(opts.simulation,'pres',opts.season,opts.grid)
113       simu_interp = ml2pl(simu,simu_pres)
114       simu = ccdo(simu_interp, operator='intlevel,'+tmplev)
115       print "cfile(simu) = ",cfile(simu)
116       reff=svsg(opts.reference,var_interp,opts.season,opts.grid)
117       reff = ccdo(reff, operator='intlevel,'+tmplev)
118       print "cfile(reff) = ",cfile(reff)
119       # -- Pas besoin d'interpoler pour les OBS => fait en amont avec derive() sur les niveaux les plus courants
120       #
121    # -- Si ce n'est pas le cas
122    else:
123       # Get the model and the reference
124       if opts.root:
125          simu=svsg(opts.simulation,variable,opts.season,opts.grid, root=opts.root)
126       else:
127          simu=svsg(opts.simulation,variable,opts.season,opts.grid)
128       reff=svsg(opts.reference,variable,opts.season,opts.grid)
129    #
130    # If the variable is a 3D field:
131    #  - interpolate the variable on the standard pressure levels with ml2pl (L. Guez)
132    #  - Compute the difference model-ref with diff_zonmean (computes the zonal mean lat/pressure fields,
133    #    interpolates the model on the ref, both vertically and horizontally, and returns the difference)
134    if is3d(variable) :
135       simu_pres = svsg(opts.simulation,'pres',opts.season,opts.grid)
136       simu = ml2pl(simu,simu_pres)
137       simu = zonmean(simu)
138       reff = zonmean(reff)
139       diff = diff_zonmean(simu,reff)
140    else:
141        if (opts.grid == '' ) : reff=regrid(reff,simu)
142        diff=minus(simu,reff)
143    #
144             
145    pparams = plot_params(variable,'full_field')
146    vertical_interval = 'trYMaxF=1000|trYMinF=1'
147    stringFontHeight=0.018
148    if is3d(variable):
149        pparams.update({'options':vertical_interval, 'y':'log'})
150        stringFontHeight=0.03
151    # Map for simulation
152    simu_fig=plot(simu,title="",
153                  gsnLeftString=variable,
154                  gsnCenterString=opts.simulation,
155                  gsnRightString=opts.season,
156                  gsnStringFontHeightF=stringFontHeight,
157                  mpCenterLonF=0,
158                  **pparams)
159    simu_avg=cvalue(space_average(simu))
160    #
161    # Map for reference
162    ref_fig=plot(reff,title="",
163                 gsnLeftString=variable,
164                 gsnCenterString=ref,
165                 gsnRightString=opts.season,
166                 gsnStringFontHeightF=stringFontHeight,
167                 mpCenterLonF=0,
168                 **pparams)
169    ref_avg=cvalue(space_average(reff))
170    #
171    # Bias (or difference between simulations) map
172    if (ref == 'OBS' ) : p=plot_params(variable,'bias')
173    else:                p=plot_params(variable,'model_model')
174    tmp_aux_params = plot_params(variable,'full_field')
175    scale = 1.0 ; offset = 0.0
176    if 'offset' in tmp_aux_params or 'scale' in tmp_aux_params:
177       if 'offset' in tmp_aux_params:
178          offset = tmp_aux_params['offset']
179       else:
180          offset=0.0
181       if 'scale' in tmp_aux_params:
182          scale = tmp_aux_params['scale']
183       else:
184          scale=1.0
185       wreff = apply_scale_offset(reff,scale,offset)
186       wsimu = apply_scale_offset(simu,scale,offset)
187    else:
188       wreff = reff
189       wsimu = simu
190    #
191    if is3d(variable):
192        p.update({'options':vertical_interval})
193    if variable in ['ua','va','ta','hus','hur']:
194        tmp_levs = tmp_aux_params['colors']
195        p.update({'contours':tmp_levs, 'y':'log'})
196        diff_fig=plot(diff,wreff,title="", format='png', mpCenterLonF=0,
197                  gsnLeftString=variable,
198                  gsnCenterString=opts.simulation+' - '+ref,
199                  gsnRightString=opts.season,
200                  gsnStringFontHeightF=stringFontHeight,
201                  aux_options='cnLineThicknessF=2|cnLineLabelsOn=True', **p)
202    else:
203        p.update({'contours':1})
204        diff_fig=plot(diff,title="", format='png', mpCenterLonF=0,
205                  gsnLeftString=variable,
206                  gsnCenterString=opts.simulation+' - '+ref,
207                  gsnRightString=opts.season,
208                  gsnStringFontHeightF=stringFontHeight,
209                  **p)
210
211    #
212    # Bias mean value, and RMSD/RMSE
213    diff_avg=cvalue(space_average(diff))
214    rmsd=math.sqrt(cvalue(space_average(ccdo(diff,operator='-b F64 sqr'))))
215    #
216    # Zonal means
217    if not is3d(variable):
218        # -- apply a mask corresponding to the reference
219        mask = fdiv(reff,reff)
220        msimu = fmul(wsimu,mask)
221        # -- Compute the zonal mean
222        zmean=ccdo(msimu, operator='zonmean')
223        ref_zmean=ccdo(wreff, operator='zonmean')
224        #
225        sim=opts.simulation
226        #if variable in ['zg500']:
227        #   ref_zmean = ccdo(ref_zmean,operator='-b F32 mulc,1')
228        #   zmean = ccdo(zmean,operator='-b F32 mulc,1')
229        #zmean_fig=curves(cens([sim,ref],zmean,ref_zmean),
230        zmean_fig=curves(cens({sim:zmean,ref:ref_zmean}),
231                         title="",
232                         lgcols=3,
233                         options=#'tiYAxisString=""|'+\
234                                 #'+\'+\
235                                 'tmYROn=True|'+\
236                                 'tmYRBorderOn=True|'+\
237                                 'tmYLOn=False|'+\
238                                 'tmYUseRight=True|'+\
239                                 'vpXF=0|'+\
240                                 'vpWidthF=0.66|'+\
241                                 'vpHeightF=0.33|'+\
242                                 'tmYRLabelsOn=True|'+\
243                                 'tmXBLabelFontHeightF=0.018|'+\
244                                 'tmYLLabelFontHeightF=0.016|'+\
245                                 'lgLabelFontHeightF=0.018|'+\
246                                 #'pmLegendSide=Bottom|'+\
247                                 'pmLegendOrthogonalPosF=-0.32|'+\
248                                 'pmLegendParallelPosF=1.0|'+\
249                                 'tmXMajorGrid=True|'+\
250                                 'tmYMajorGrid=True|'+\
251                                 'tmXMajorGridLineDashPattern=2|'+\
252                                 'tmYMajorGridLineDashPattern=2|'+\
253                                 'xyLineThicknessF=8|'+\
254                                 'gsnLeftString='+variable+'|'+\
255                                 'gsnCenterString='+opts.simulation+' vs '+ref+'|'+\
256                                 'gsnRightString='+opts.season+'|'+\
257                                 'gsnStringFontHeightF='+str(stringFontHeight))
258    #    # Composite figure
259    if is3d(variable):
260        page=cpage([[simu_fig,ref_fig],[diff_fig,None]],orientation='landscape', page_trim=True, fig_trim=True)
261        pdf_page=cpage([[simu_fig,ref_fig],[diff_fig,None]],orientation='landscape',
262                   page_trim=True, fig_trim=True, format='pdf', title=variable+" "+opts.simulation)
263    else:
264        page=cpage([[simu_fig,ref_fig],[diff_fig,zmean_fig]],orientation='landscape', page_trim=True, fig_trim=True)
265        pdf_page=cpage([[simu_fig,ref_fig],[diff_fig,zmean_fig]],orientation='landscape', 
266                   page_trim=True, fig_trim=True, format='pdf', title=variable+" "+opts.simulation)
267    pdfargs.append(cfile(pdf_page))
268    #
269    thumbnail_size = 200
270    if is3d(variable):
271            index+=open_line(varlongname(variable)+' ('+variable+')')+\
272                    cell("%.2g"%diff_avg,cfile(diff_fig))+\
273                    cell("%.2g"%rmsd,cfile(diff_fig))+\
274                    cell(simu,cfile(simu_fig),thumbnail=thumbnail_size,hover=False)+\
275                    cell(ref,cfile(ref_fig),thumbnail=thumbnail_size,hover=False)+\
276                    cell(text_diff,cfile(diff_fig),thumbnail=thumbnail_size,hover=False)+\
277                    ' '+\
278                    cell('page',cfile(page),thumbnail=thumbnail_size,hover=False)+\
279                    cell('Pdf',cfile(pdf_page))
280            close_line()
281    else:
282            index+=open_line(varlongname(variable)+' ('+variable+')')+\
283                    cell("%.2g"%diff_avg,cfile(diff_fig))+\
284                    cell("%.2g"%rmsd,cfile(diff_fig))+\
285                    cell(simu,cfile(simu_fig),thumbnail=thumbnail_size,hover=False)+\
286                   cell(ref,cfile(ref_fig),thumbnail=thumbnail_size,hover=False)+\
287                    cell(text_diff,cfile(diff_fig),thumbnail=thumbnail_size,hover=False)+\
288                    cell('zonal mean',cfile(zmean_fig),thumbnail=thumbnail_size,hover=False)+\
289                    cell('page',cfile(page),thumbnail=thumbnail_size,hover=False)+\
290                    cell('Pdf',cfile(pdf_page))
291            close_line()
292#
293# Finalisons l'index html
294index += close_table()
295index += trailer()
296#out="index_example.html"
297out="index_example_"+opts.season+"_"+opts.simulation+".html"
298with open(out,"w") as filout : filout.write(index)
299#
300# Creation du Pdf multi-pages
301comm=subprocess.Popen(pdfargs, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
302#
303import os,os.path ; 
304# print("Attendez un bon peu : lancemement de firefox sur Ciclad....")
305# os.system("firefox file://"+os.path.abspath(os.path.curdir)+"/"+out+"&")
Note: See TracBrowser for help on using the repository browser.