Subversion Repositories lagranto.ecmwf

Rev

Rev 35 | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 35 Rev 36
Line 1... Line 1...
1 
#!/usr/bin/env python
 1 
#!/usr/bin/python
2 
# -*- coding:utf-8 -*-
 - 
 
3 
 
 2 
 
4 
# ---------------------------------------------------------------------
 3 
# =============================================================================
5 
# import modules
 4 
# Quick Viewer for meteorological fields from ECMWF
6 
# ---------------------------------------------------------------------
 5 
# Michael Sprenger, Sebastian Limbach [ based on IWAL], Winter 2012/13
- 
 
 6 
# =============================================================================
7 
 
 7 
 
- 
 
 8 
#  Imports
8 
from __future__ import print_function
 9 
import os
9 
import sys
 10 
import sys
10 
import argparse
 - 
 
11 
from dypy.lagranto import Tra
 - 
 
12 
from dypy.plotting import Mapfigure
 11 
from   optparse import OptionParser
13 
import matplotlib.pyplot as plt
 12 
import subprocess
14 
import matplotlib
 13 
import datetime
- 
 
 14 
 
15 
matplotlib.rcParams['backend'] = "Qt4Agg"
 15 
# The absolute path of the test-scripts
- 
 
 16 
SCRIPT_PATH = os.getenv('DYN_TOOLS') + '/quickview/'
- 
 
 17 
 
16 
# ---------------------------------------------------------------------
 18 
# ----------------------------------------------------------------------------
17 
# main function
 19 
# Help page / show colorbar only
18 
# ---------------------------------------------------------------------
 20 
# ----------------------------------------------------------------------------
- 
 
 21 
 
- 
 
 22 
# Write the help page on demand
- 
 
 23 
if ( sys.argv[1] == 'help' ):
19 
 
 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 ( '--------------------------------------------------------------------------------------------------------' )
20 
 
 57 
 
21 
def main(argv):
 58 
    quit()
- 
 
 59 
 
- 
 
 60 
elif ( sys.argv[1] == 'colorbar' ):
- 
 
 61 
 
- 
 
 62 
    figname = 'quickview.color'
- 
 
 63 
    figpath = os.getcwd() + '/'
22 
 
 64 
 
23 
    # parameters
 65 
    print figname, figpath
24 
 
 66 
 
- 
 
 67 
    wsname  = 'testws'
25 
    domain = None
 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 
 
26 
    var = ''
 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'],
27 
    eu = [-15, 30, 30, 75]
 91 
         'V10M'    : [ 'B', 'dansgaard'],
- 
 
 92 
         'TD2M'    : [ 'P', 'dansgaard'],
- 
 
 93 
         'T2M'     : [ 'P', 'dansgaard'],
- 
 
 94 
         'SKT'     : [ 'B', 'dansgaard'],
28 
    nh = [-180, 180, 0, 90]
 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'   ]
29 
    filename = ""
 115 
       }
- 
 
 116 
 
- 
 
 117 
# ----------------------------------------------------------------------------
- 
 
 118 
# Get the arguments
- 
 
 119 
# ----------------------------------------------------------------------------
30 
 
 120 
 
- 
 
 121 
# Parse optional parameters
31 
    parser = argparse.ArgumentParser(
 122 
parser = OptionParser()
32 
        description="Display trajectories on world map")
 123 
parser.add_option("-o",  "--outfile",    type='string',       dest="outfile",    default = 'NONE'  )
33 
    parser.add_argument("-d", "--domain", type=str,
 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'  )
34 
                        help="Choose domain to plot nh (north hemisphere) \
 127 
parser.add_option("-c",  "--color",      type='string',       dest="color",      default = 'NONE'  )
35 
                        / eu (europe) / lon1,lon2,lat1,lat2")
 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   )
36 
    parser.add_argument("-v", "--variable", type=str, help="Choose a variable")
 130 
parser.add_option("-r",  "--showremove", action="store_true", dest="showremove", default = False   )
37 
    parser.add_argument("-s", "--save", type=str, help="Save to file")
 131 
parser.add_option("-l",  "--labelbar",   action="store_true", dest="labelbar",   default = False   )
38 
    parser.add_argument("file", type=str, help="filename")
 132 
parser.add_option("-g",  "--geoproj",    type='string',       dest="projection", default = 'latlon')
- 
 
 133 
 
39 
    args = parser.parse_args()
 134 
(options, args) = parser.parse_args()
40 
 
 135 
 
- 
 
 136 
# Get mandatory arguments (inpfile / field)
- 
 
 137 
inpfile = args[0]
- 
 
 138 
if ( len(args) == 2 ) :
- 
 
 139 
    field   = args[1]
- 
 
 140 
else:
41 
    # read options
 141 
    field = 'T'
42 
 
 142 
 
- 
 
 143 
# Adjusts defaults
- 
 
 144 
if options.outfile   == 'NONE' :
43 
    filename = args.file
 145 
    outfile = inpfile
- 
 
 146 
else:
44 
    if args.domain is None:
 147 
    outfile = options.outfile
- 
 
 148 
 
- 
 
 149 
# Set some internal parameters
45 
        pass
 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'  )
46 
    elif args.domain == "nh":
 193 
    if ( i0 != -1 ):
47 
        domain = nh
 194 
        field_unit  = 'hPa'
- 
 
 195 
        field_level = field_level[0:i0]
48 
    elif args.domain == "eu":
 196 
    elif ( i1 != -1 ):
49 
        domain = eu
 197 
        field_unit = 'K'
- 
 
 198 
        field_level = field_level[0:i1]
50 
    else:
 199 
    else:
- 
 
 200 
        field_unit = 'ML'
51 
        domain = [float(coo) for coo in args.domain.split(',')]
 201 
        field_level = field_level
52 
 
 202 
 
53 
    # -----------------------------------------
 - 
 
54 
    # read the trajectories
 - 
 
55 
    # -----------------------------------------
 - 
 
56 
 
 - 
 
57 
    trajs = Tra(filename)
 - 
 
58 
 
 - 
 
59 
    print(trajs.shape)
 203 
print field_unit
60 
    if trajs.ntra > 500:
 - 
 
61 
        it = trajs.ntra / 500
 - 
 
62 
        trajs.set_array(trajs[::it, :])
 - 
 
63 
        printstring = 'Only {} trajectories are plotted (every {})'
 - 
 
64 
        print(printstring.format(trajs.ntra, it))
 - 
 
65 
 
 - 
 
66 
    if domain is None:
 - 
 
67 
        domain = [trajs['lon'].min(), trajs['lon'].max(),
 - 
 
68 
                  trajs['lat'].min(), trajs['lat'].max()]
 - 
 
69 
    # set default parameters
 - 
 
70 
    if args.variable is None:
 - 
 
71 
        var = trajs.dtype.names[-1]
 - 
 
72 
        printstring = "Fields available for plotting: \n {0}"
 - 
 
73 
        print(printstring.format(trajs.dtype.names))
 - 
 
74 
    else:
 - 
 
75 
        var = args.variable
 - 
 
76 
 
 - 
 
77 
    # plot trajectories
 - 
 
78 
 
 - 
 
79 
    m = Mapfigure(domain=domain, resolution='l')
 - 
 
80 
 
 - 
 
81 
    fig = plt.figure()
 - 
 
82 
    m.ax = fig.add_subplot(111)
 - 
 
83 
    m.drawmap()
 - 
 
84 
    lc = m.plot_traj(trajs, var, rasterized=True)
 - 
 
85 
    m.colorbar(lc, 'right', label=var)
 - 
 
86 
 
 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()
87 
    if args.save is None:
 226 
    if ( weekday == 0 ):
88 
        plt.show()
 227 
        weekday = 6
89 
    else:
 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):
90 
        plt.savefig(args.save, bbox_inches='tight')
 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 ]))
91 
 
 492 
 
- 
 
 493 
# ----------------------------------------------------------------------------
- 
 
 494 
# Skew-T/log-P (skewtlogp)
- 
 
 495 
# ----------------------------------------------------------------------------
92 
 
 496 
 
93 
if __name__ == '__main__':
 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 
 
94 
    main(sys.argv[1:])
 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)