Subversion Repositories lagranto.arpege

Rev

Details | Last modification | View Log | RSS feed

Rev Author Line No. Line
3 michaesp 1
#!/usr/bin/python
2
 
3
# =============================================================================
4
# Quick Viewer for meteorological fields from ECMWF
5
# Michael Sprenger, Sebastian Limbach [ based on IWAL], Winter 2012/13  
6
# =============================================================================
7
 
8
#  Imports
9
import os
10
import sys
11
from   optparse import OptionParser
12
import subprocess
13
import datetime
14
 
15
# The absolute path of the test-scripts
16
SCRIPT_PATH = os.getenv('DYN_TOOLS') + '/quickview/'
17
 
18
# ----------------------------------------------------------------------------
19
# Help page / show colorbar only
20
# ----------------------------------------------------------------------------
21
 
22
# Write the help page on demand
23
if ( sys.argv[1] == 'help' ):
24
 
25
    print ( '--------------------------------------------------------------------------------------------------------' )
26
    print ( 'Usage: [1] qv help                                                                                      ' ) 
27
    print ( '       [2] qv data field [options]       -> horizontal and vertical cross sections                      ' )
28
    print ( '       [3] qv data [options]             -> skew-T/log-P diagram and satellite imagery                  ' )
29
    print ( '       [4] qv colorbar                   -> show colorbar only                                          ' )
30
    print ( '                                                                                                        ' )
31
    print ( 'Arguments:                                                                                              ' )
32
    print ( '                                                                                                        ' )
33
    print ( '      data = P{date}           : explicit filename                                                      ' )
34
    print ( '      data = an.{date}         : operational ECMWF analysis                                             ' )
35
    print ( '      data = fc.{date}         : most recent operational ECMWF analysis for {date}                      ' )
36
    print ( '      data = fc.{date0}.{date} : operational ECMWF initialized at {date0} analysis for {date}           ' )
37
    print ( '      data = ir|wv.date        : IR and WV satellite imagery                                            ' )
38
    print ( '                                                                                                        ' )
39
    print ( '      field = Q@200hPa         : e.g. specific humidity (Q) at 200 hPa                                  ' )
40
    print ( '      field = PV@320K          : e.g. potential vorticity (PV) at 320 K                                 ' )
41
    print ( '      field = T@15             : e.g. temperature (T) at model level 15                                 ' )
42
    print ( '      field = SLP              : e.. reduced surface pressure                                           ' )
43
    print ( '                                                                                                        ' )
44
    print ( 'Options:                                                                                                ' )
45
    print ('                                                                                                         ' )
46
    print ('        -o  figname                              : name of output figure                                 ' )
47
    print ('        -d  lonmin,lonmax,latmin,latmax|region   : horizontal domain (region:europe,natlantic,nhem,shem) ' )
48
    print ('        -g  np|sp|latlon                         : geographic projection (np:north pole, sp:south pole, latlon (default)')      
49
    print ('        -v  lon1,lat1,lon2,lat2[,lev1,lev2,type] : vertical cross section (pype:gc,linear,ns,ew)         ' )
50
    print ('        -p  lon,lat[,lev1,lev2]                  : vertical profile (skew-T/log-P)                       ' )
51
    print ('        -c  cmin,cmax[,mode,#cols|colname]       : specification of colortable                           ' )
52
    print ('        -n  npixel                               : pixel size of figure (default=512)                    ' )
53
    print ('        -s                                       : show figure                                           ' )
54
    print ('        -r                                       : remove figure after showing it                        ' )
55
    print ('        -l                                       : show label/color bar                                  ' )
56
    print ( '--------------------------------------------------------------------------------------------------------' )
57
 
58
    quit()
59
 
60
elif ( sys.argv[1] == 'colorbar' ):
61
 
62
    figname = 'quickview.color'
63
    figpath = os.getcwd() + '/'
64
 
65
    print figname, figpath
66
 
67
    wsname  = 'testws'
68
    npixel  = 512
69
    os.system('python ' + SCRIPT_PATH + '/ext_scripts/colorbar.py ' +
70
          ' '.join([str(npixel), wsname, figpath, figname])) 
71
 
72
    cmd = 'xv ' + figpath + '/quickview.color.png;'
73
    subprocess.Popen(cmd,stdout=subprocess.PIPE,shell=True)
74
 
75
    quit()
76
 
77
# ----------------------------------------------------------------------------
78
# Set ERA-Interim dictionary
79
# ----------------------------------------------------------------------------
80
 
81
erai = { 'T'       : [ 'P', 'dansgaard'],
82
         'OMEGA'   : [ 'P', 'dansgaard'],
83
         'U'       : [ 'P', 'dansgaard'],
84
         'V'       : [ 'P', 'dansgaard'],
85
         'SLP'     : [ 'P', 'dansgaard'],
86
         'PS'      : [ 'P', 'dansgaard'],
87
         'Q'       : [ 'P', 'dansgaard'],
88
         'LWC'     : [ 'P', 'dansgaard'],
89
         'IWC'     : [ 'P', 'dansgaard'],
90
         'U10M'    : [ 'B', 'dansgaard'],
91
         'V10M'    : [ 'B', 'dansgaard'],
92
         'TD2M'    : [ 'P', 'dansgaard'],
93
         'T2M'     : [ 'P', 'dansgaard'],
94
         'SKT'     : [ 'B', 'dansgaard'],
95
         'G10M'    : [ 'G', 'dansgaard'],
96
         'T2MAX'   : [ 'G', 'dansgaard'],
97
         'T2MIN'   : [ 'G', 'dansgaard'],
98
         'RTOT'    : [ 'R', 'dansgaard'],
99
         'LSP'     : [ 'R', 'dansgaard'],
100
         'CP'      : [ 'R', 'dansgaard'],
101
         'SF'      : [ 'R', 'dansgaard'],
102
         'SSHF'    : [ 'R', 'dansgaard'],
103
         'O3'      : [ 'O', 'litho'    ],
104
         'TROPO'   : [ 'L', 'litho'    ],
105
         'LABEL'   : [ 'L', 'litho'    ],
106
         'CAPE_MU' : [ 'C', 'litho'    ],
107
         'CIN_MU'  : [ 'C', 'litho'    ],
108
         'CAPE_ML' : [ 'C', 'litho'    ],
109
         'CIN_ML'  : [ 'C', 'litho'    ],
110
         'Z'       : [ 'Z', 'dansgaard'],
111
         'TH'      : [ 'S', 'thermo'   ],
112
         'PV'      : [ 'S', 'thermo'   ],
113
         'RH'      : [ 'S', 'thermo'   ],
114
         'THE'     : [ 'S', 'thermo'   ]
115
       }    
116
 
117
# ----------------------------------------------------------------------------
118
# Get the arguments
119
# ----------------------------------------------------------------------------
120
 
121
# Parse optional parameters
122
parser = OptionParser()
123
parser.add_option("-o",  "--outfile",    type='string',       dest="outfile",    default = 'NONE'  )
124
parser.add_option("-d",  "--domain",     type='string',       dest="domain",     default = 'NONE'  )
125
parser.add_option("-p",  "--profile",    type='string',       dest="profile",    default = 'NONE'  )
126
parser.add_option("-v",  "--vertical",   type='string',       dest="vertical",   default = 'NONE'  )
127
parser.add_option("-c",  "--color",      type='string',       dest="color",      default = 'NONE'  )
128
parser.add_option("-n",  "--npixel",     type='int',          dest="npixel",     default = 512     )
129
parser.add_option("-s",  "--show",       action="store_true", dest="show",       default = False   )
130
parser.add_option("-r",  "--showremove", action="store_true", dest="showremove", default = False   )
131
parser.add_option("-l",  "--labelbar",   action="store_true", dest="labelbar",   default = False   ) 
132
parser.add_option("-g",  "--geoproj",    type='string',       dest="projection", default = 'latlon')
133
 
134
(options, args) = parser.parse_args()
135
 
136
# Get mandatory arguments (inpfile / field) 
137
inpfile = args[0]
138
if ( len(args) == 2 ) :
139
    field   = args[1]
140
else:
141
    field = 'T'
142
 
143
# Adjusts defaults
144
if options.outfile   == 'NONE' :
145
    outfile = inpfile
146
else:
147
    outfile = options.outfile
148
 
149
# Set some internal parameters 
150
wsname     = 'testws'
151
 
152
# Set name of output figures
153
figname = outfile
154
figpath = os.getcwd()+'/'
155
 
156
# Decide about the plotting mode
157
if any( elem == '-d' for elem in sys.argv ):
158
    print "Mode : horizontal cross scection"
159
    mode = 'horizontal'
160
 
161
elif any( elem == '-v' for elem in sys.argv ):
162
    print "Mode : vertical cross scection"
163
    mode = 'vertical'
164
 
165
elif any( elem == '-p' for elem in sys.argv ):
166
    print "Mode : vertical profile"
167
    mode = 'profile'
168
 
169
else:
170
    print "Mode : horizontal cross section (default)"
171
    mode = 'horizontal'
172
 
173
# Check that the number of pixel is reasonable
174
if ( options.npixel > 2000 ):
175
    sys.exit('number of pixels too large');
176
 
177
# ----------------------------------------------------------------------------
178
# Set the meteorological field 
179
# ----------------------------------------------------------------------------
180
 
181
field_split = field.split('@')
182
 
183
if len(field_split) == 1 :
184
    field_fieldname = field_split[0]
185
    field_level     = 0
186
    field_unit      = 'nil'
187
 
188
elif len(field_split) == 2 :
189
    field_fieldname = field_split[0]
190
    field_level     = field_split[1]
191
    i0 = field_level.find('hPa')
192
    i1 = field_level.find('K'  )
193
    if ( i0 != -1 ):
194
        field_unit  = 'hPa'
195
        field_level = field_level[0:i0]
196
    elif ( i1 != -1 ):
197
        field_unit = 'K'
198
        field_level = field_level[0:i1]
199
    else:
200
        field_unit = 'ML'
201
        field_level = field_level
202
 
203
print field_unit
204
 
205
# ----------------------------------------------------------------------------
206
# Set the data source + check whether field is available
207
# ----------------------------------------------------------------------------
208
 
209
inpfile_split = inpfile.split('.')
210
 
211
if len(inpfile_split) == 1 :
212
    inpfile_filename = os.path.split(inpfile_split[0])[1]
213
    inpfile_folder   = os.path.split(inpfile_split[0])[0]
214
    if inpfile_folder == '':
215
        inpfile_folder = os.getcwd()
216
    inpfile_folder = inpfile_folder + '/'
217
 
218
elif ( inpfile_split[0] == 'an' ) & ( len(inpfile_split) == 2 ):
219
    inpfile_filename = 'P'+inpfile_split[1]
220
    inpfile_folder   = '/net/litho/atmosdyn/ec.analysis/'
221
 
222
elif  ( inpfile_split[0] == 'fc' ) & ( len(inpfile_split) == 2 ):
223
    inpfile_filename = 'P'+inpfile_split[1]
224
    today            = datetime.date.today()
225
    weekday          = today.weekday()
226
    if ( weekday == 0 ):
227
        weekday = 6
228
    else:
229
        weekday = weekday -1
230
    daynames         = ( 'Mon','Tue','Wed','Thu','Fri','Sat','Sun' )
231
    inpfile_folder   = '/net/litho/atmosdyn/ec.forecast/' + daynames[weekday] + '/'
232
 
233
elif  ( inpfile_split[0] == 'fc' ) & ( len(inpfile_split) == 3 ):
234
    inpfile_filename = 'P'+inpfile_split[2]
235
    inpfile_folder   = '/net/litho/atmosdyn/ec.forecast/' + inpfile_split[1] + '/'
236
 
237
elif ( inpfile_split[0] == 'ei' ) & ( len(inpfile_split) == 2 ):
238
 
239
    inpfile_filename = erai[field_fieldname][0]+inpfile_split[1]
240
    yyyy             = inpfile_filename[1:5]
241
    mm               = inpfile_filename[5:7]
242
    inpfile_folder   = '/net/'+erai[field_fieldname][1]+'/atmosdyn/erainterim/cdf/'+yyyy+'/'+mm+'/'
243
 
244
else:
245
    sys.exit('Unsupported input file ' + inpfile);
246
 
247
# Check whether the input file is available
248
if os.path.exists(inpfile_folder+inpfile_filename):
249
    print 'input file found: ' + inpfile_folder + inpfile_filename
250
else:
251
    sys.exit( 'input file is missing : ' + inpfile_folder + inpfile_filename )
252
 
253
# Check whether field is available on input file
254
cmd        = 'getvars ' + inpfile_folder + '/' + inpfile_filename
255
proc       = subprocess.Popen(cmd,stdout=subprocess.PIPE,shell=True)
256
(out, err) = proc.communicate()  
257
ok         = any( (elem.strip() == field_fieldname) for elem in out.split('\n') ) 
258
 
259
# If not available, try to calculate it with p2s
260
if ( ok == False ):
261
    print 'Try to calculate '+field_fieldname+' with p2s'
262
    cmd = 'p2s ' + inpfile_folder + '/' + inpfile_filename
263
    print 'NOT YET IMPLEMTED'
264
 
265
if ok == False:
266
    sys.exit('variable '+field_fieldname+' missing on '+inpfile_filename);
267
 
268
# ----------------------------------------------------------------------------
269
# Set the color table
270
# ----------------------------------------------------------------------------
271
 
272
color_split = options.color.split(',')
273
 
274
if len(color_split) == 5:
275
    color_cmin = color_split[0]
276
    color_cmax = color_split[1]
277
    color_mode = color_split[2]
278
    color_ncol = color_split[3]
279
    color_name = color_split[4]
280
 
281
elif len(color_split) == 4:
282
    color_cmin = color_split[0]
283
    color_cmax = color_split[1]
284
    color_mode = color_split[2]
285
    color_ncol = color_split[3]
286
    color_name = 'adjust'
287
 
288
elif len(color_split) == 3:
289
    color_cmin = color_split[0]
290
    color_cmax = color_split[1]
291
    color_mode = color_split[2]
292
    color_ncol = 18
293
    color_name = 'hotcold_18lev'
294
 
295
elif len(color_split) == 2:
296
    color_cmin = color_split[0]
297
    color_cmax = color_split[1]
298
    color_mode = 'fill'
299
    color_ncol = '18'
300
    color_name = 'hotcold_18lev'
301
 
302
elif (color_split[0] == 'PV.default') | ( color_split[0] == 'NONE' ) & ( field_fieldname == 'PV' ):
303
    color_cmin = '-0.5'
304
    color_cmax = '10.'
305
    color_mode = 'fill'
306
    color_ncol = 'adjust'
307
    color_name = 'PV.default'
308
 
309
elif color_split[0] == 'NONE':
310
    color_cmin = 'p95'
311
    color_cmax = 'p05'
312
    color_mode = 'fill'
313
    color_ncol = '18'
314
    color_name = 'hotcold_18lev'
315
 
316
# ----------------------------------------------------------------------------
317
# Set the horizontal domain
318
# ----------------------------------------------------------------------------
319
 
320
# Set the geographical projection 
321
projection = options.projection
322
print projection      
323
 
324
domain_split = options.domain.split(',')
325
 
326
print domain_split
327
 
328
if len(domain_split) == 4:
329
    domain_xmin = domain_split[0]
330
    domain_xmax = domain_split[1]
331
    domain_ymin = domain_split[2]
332
    domain_ymax = domain_split[3]
333
 
334
elif ( len(domain_split) == 1 ) & (options.domain == 'europe' ):
335
    domain_xmin = -20
336
    domain_xmax = 40
337
    domain_ymin = 30
338
    domain_ymax = 75
339
 
340
elif ( len(domain_split) == 1 ) & (options.domain == 'natlantic' ):
341
    domain_xmin = -100
342
    domain_xmax = 10
343
    domain_ymin = 10
344
    domain_ymax = 75
345
 
346
elif ( len(domain_split) == 1 ) & (options.domain == 'global' ):
347
    domain_xmin = -180
348
    domain_xmax = 180
349
    domain_ymin = -90
350
    domain_ymax = 90
351
 
352
elif ( len(domain_split) == 1 ) & (options.domain == 'nhem' ):
353
    domain_xmin = -180
354
    domain_xmax = 180
355
    domain_ymin = 0
356
    domain_ymax = 90
357
 
358
elif ( len(domain_split) == 1 ) & (options.domain == 'shem' ):
359
    domain_xmin = -180
360
    domain_xmax = 180
361
    domain_ymin = -90
362
    domain_ymax = 0
363
 
364
elif ( len(domain_split) == 2 ):
365
    domain_xmin = -180
366
    domain_xmax = 180
367
    domain_ymin = domain_split[0]
368
    domain_ymax = domain_split[1]
369
 
370
elif ( projection == 'np' ):
371
    domain_xmin = -180
372
    domain_xmax = 180
373
    domain_ymin = 30
374
    domain_ymax = 90
375
 
376
elif ( projection == 'sp' ):
377
    domain_xmin = -180
378
    domain_xmax = 180
379
    domain_ymin = -90
380
    domain_ymax = -30
381
 
382
else:
383
    domain_xmin = 'adjust'
384
    domain_xmax = 'adjust'
385
    domain_ymin = 'adjust'
386
    domain_ymax = 'adjust'
387
 
388
# Check that the geographical domain is reasonable
389
if ( domain_xmin != 'adjust' ):
390
    if (domain_xmin >= domain_xmax) | (domain_ymin >= domain_ymax):
391
        sys.exit('domain problem:')
392
 
393
# ----------------------------------------------------------------------------
394
# Set the settings for vertical cross section
395
# ----------------------------------------------------------------------------
396
 
397
vertical_split = options.vertical.split(',')
398
 
399
if len(vertical_split) == 7:
400
    vertical_lon1 = vertical_split[0]
401
    vertical_lat1 = vertical_split[1]
402
    vertical_lon2 = vertical_split[2]
403
    vertical_lat2 = vertical_split[3]
404
    vertical_lev1 = vertical_split[4]
405
    vertical_lev2 = vertical_split[5]
406
    vertical_type = vertical_split[6]
407
 
408
if len(vertical_split) == 6:
409
    vertical_lon1 = vertical_split[0]
410
    vertical_lat1 = vertical_split[1]
411
    vertical_lon2 = vertical_split[2]
412
    vertical_lat2 = vertical_split[3]
413
    vertical_lev1 = vertical_split[4]
414
    vertical_lev2 = vertical_split[5]
415
    vertical_type = 'linear'
416
 
417
if len(vertical_split) == 4:
418
    vertical_lon1 = vertical_split[0]
419
    vertical_lat1 = vertical_split[1]
420
    vertical_lon2 = vertical_split[2]
421
    vertical_lat2 = vertical_split[3]
422
    vertical_lev1 = 1030
423
    vertical_lev2 = 100
424
    vertical_type = 'linear'
425
 
426
elif ( len(vertical_split) == 1 ) & (options.vertical == 'zurich.ns' ):
427
    vertical_lon1 = 8.538
428
    vertical_lat1 = 47.3690 - 7.5
429
    vertical_lon2 = 8.538
430
    vertical_lat2 = 47.3690 + 7.5
431
    vertical_lev1 = 1030
432
    vertical_lev2 = 100
433
    vertical_type = 'linear'
434
 
435
# ----------------------------------------------------------------------------
436
# Set the profile settings
437
# ----------------------------------------------------------------------------
438
 
439
profile_split = options.profile.split(',')
440
 
441
if len(profile_split) == 4:
442
    profile_lon  = profile_split[0]
443
    profile_lat  = profile_split[1]
444
    profile_lev1 = profile_split[2]
445
    profile_lev2 = profile_split[3]
446
 
447
elif len(profile_split) == 2:
448
    profile_lon  = profile_split[0]
449
    profile_lat  = profile_split[1]
450
    profile_lev1 = 1030
451
    profile_lev2 = 100
452
 
453
# ----------------------------------------------------------------------------
454
# Horizontal cross section 
455
# ----------------------------------------------------------------------------
456
 
457
if ( mode == 'horizontal' ):
458
 
459
    os.system('python ' + SCRIPT_PATH + '/ext_scripts/horizontal.py ' + 
460
        ' '.join([inpfile_filename, 
461
                  inpfile_folder, 
462
                  field_fieldname, 
463
                  str(field_level),
464
                  str(domain_xmin),str(domain_xmax),str(domain_ymin),str(domain_ymax), 
465
                  str(options.npixel), 
466
                  wsname, 
467
                  figpath,figname, 
468
                  str(color_cmin), str(color_ncol), str(color_cmax), color_name, 
469
                  field_unit,projection,color_mode]))
470
 
471
 
472
 
473
# ----------------------------------------------------------------------------
474
# Vertical cross section 
475
# ----------------------------------------------------------------------------
476
 
477
if ( mode == 'vertical' ) :
478
 
479
    os.system('python ' + SCRIPT_PATH + '/ext_scripts/vertical.py ' +
480
          ' '.join([inpfile_filename, 
481
                    inpfile_folder, 
482
                    field_fieldname,
483
                    str(vertical_lon1), str(vertical_lon2), 
484
                    str(vertical_lat1), str(vertical_lat2), 
485
                    str(vertical_lev1), str(vertical_lev2),
486
                    vertical_type, 
487
                    str(options.npixel), 
488
                    wsname, 
489
                    figpath, figname,
490
                    str(color_cmin), str(color_ncol), str(color_cmax), 
491
                    color_name,color_mode ]))
492
 
493
# ----------------------------------------------------------------------------
494
# Skew-T/log-P (skewtlogp)
495
# ----------------------------------------------------------------------------
496
 
497
if ( mode == 'profile' ):
498
 
499
    os.system('python ' + SCRIPT_PATH + '/ext_scripts/skewtlogp.py ' +
500
          ' '.join([inpfile_filename,inpfile_folder,
501
          str(profile_lon), str(profile_lat), str(profile_lev1), str(profile_lev2),
502
          str(options.npixel), wsname, figpath, figname ]))
503
 
504
 
505
# ----------------------------------------------------------------------------
506
# Show the figure if requested
507
# ----------------------------------------------------------------------------
508
 
509
if options.show:
510
    cmd = 'xv ' + figpath + '/' + figname + '.png'
511
    subprocess.Popen(cmd,stdout=subprocess.PIPE,shell=True)
512
 
513
if options.showremove:
514
    cmd = 'xv ' + figpath + '/' + figname + '.png; rm ' + figpath + '/' + figname + '.png'
515
    subprocess.Popen(cmd,stdout=subprocess.PIPE,shell=True)
516
 
517
if options.labelbar | ( mode == 'colorbar' ):
518
 
519
    figname = 'quickview.color'
520
    os.system('python ' + SCRIPT_PATH + '/ext_scripts/colorbar.py ' +
521
          ' '.join([str(options.npixel), wsname, figpath, figname])) 
522
 
523
    if options.showremove:
524
        cmd = 'xv ' + figpath + '/quickview.color.png;' + 'rm ' + figpath + '/quickview.color.png'
525
        subprocess.Popen(cmd,stdout=subprocess.PIPE,shell=True)
526
    else:
527
        cmd = 'xv ' + figpath + '/quickview.color.png;'
528
        subprocess.Popen(cmd,stdout=subprocess.PIPE,shell=True)