Changeset 233 for trunk/MESOSCALE_DEV


Ignore:
Timestamp:
Jul 19, 2011, 2:22:19 AM (13 years ago)
Author:
aslmd
Message:

MESOSCALE: improved and clean winds.py script to make fast and powerful plots. preparing the routine to be easily adapted to sections. also: suppressed the old directory mars_lmd_new.bak [<< ancienne nouvelle physique >>]

Location:
trunk/MESOSCALE_DEV/PLOT/PYTHON
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/MESOSCALE_DEV/PLOT/PYTHON/mylib/myplot.py

    r232 r233  
     1def errormess(text):
     2    print text
     3    exit()
     4    return
     5
     6def whatkindfile (nc):
     7    if 'controle' in nc.variables:   typefile = 'gcm'
     8    elif 'vert' in nc.variables:     typefile = 'mesoapi'
     9    elif 'U' in nc.variables:        typefile = 'meso'
     10    elif 'HGT_M' in nc.variables:    typefile = 'geo'
     11    else:                            errormess("whatkindfile: typefile not supported.")
     12    return typefile
     13
     14def getfield (nc,var):
     15    ## this allows to get much faster (than simply referring to nc.variables[var])
     16    dimension = len(nc.variables[var].dimensions)
     17    if dimension == 2:    field = nc.variables[var][:,:]
     18    elif dimension == 3:  field = nc.variables[var][:,:,:]
     19    elif dimension == 4:  field = nc.variables[var][:,:,:,:]
     20    return field
     21
     22def reducefield (input,d4=None,d3=None,d2=None,d1=None):
     23    ### do it the reverse way to be compliant with netcdf "t z y x" or "t y x" or "y x"
     24    ### it would be actually better to name d4 d3 d2 d1 as t z y x
     25    import numpy as np
     26    dimension = np.array(input).ndim
     27    shape = np.array(input).shape
     28    print 'dim,shape: ',dimension,shape
     29    output = input
     30    error = False
     31    if dimension == 2:
     32        if   d2 >= shape[0]: error = True
     33        elif d1 >= shape[1]: error = True
     34        elif d1 is not None and d2 is not None:  output = input[d2,d1]
     35        elif d1 is not None:         output = input[:,d1]
     36        elif d2 is not None:         output = input[d2,:]
     37    elif dimension == 3:
     38        if   d4 >= shape[0]: error = True
     39        elif d2 >= shape[1]: error = True
     40        elif d1 >= shape[2]: error = True
     41        elif d4 is not None and d2 is not None and d1 is not None:  output = input[d4,d2,d1]
     42        elif d4 is not None and d2 is not None:         output = input[d4,d2,:]
     43        elif d4 is not None and d1 is not None:         output = input[d4,:,d1]
     44        elif d2 is not None and d1 is not None:         output = input[:,d2,d1]
     45        elif d1 is not None:                output = input[:,:,d1]
     46        elif d2 is not None:                output = input[:,d2,:]
     47        elif d4 is not None:                output = input[d4,:,:]
     48    elif dimension == 4:
     49        if   d4 >= shape[0]: error = True
     50        elif d3 >= shape[1]: error = True
     51        elif d2 >= shape[2]: error = True
     52        elif d1 >= shape[3]: error = True
     53        elif d4 is not None and d3 is not None and d2 is not None and d1 is not None:  output = input[d4,d3,d2,d1]
     54        elif d4 is not None and d3 is not None and d2 is not None:         output = input[d4,d3,d2,:]
     55        elif d4 is not None and d3 is not None and d1 is not None:         output = input[d4,d3,:,d1]
     56        elif d4 is not None and d2 is not None and d1 is not None:         output = input[d4,:,d2,d1]
     57        elif d3 is not None and d2 is not None and d1 is not None:         output = input[:,d3,d2,d1]
     58        elif d4 is not None and d3 is not None:                output = input[d4,d3,:,:]
     59        elif d4 is not None and d2 is not None:                output = input[d4,:,d2,:]
     60        elif d4 is not None and d1 is not None:                output = input[d4,:,:,d1]
     61        elif d3 is not None and d2 is not None:                output = input[:,d3,d2,:]
     62        elif d3 is not None and d1 is not None:                output = input[:,d3,:,d1]
     63        elif d2 is not None and d1 is not None:                output = input[:,:,d2,d1]
     64        elif d1 is not None:                       output = input[:,:,:,d1]
     65        elif d2 is not None:                       output = input[:,:,d2,:]
     66        elif d3 is not None:                       output = input[:,d3,:,:]
     67        elif d4 is not None:                       output = input[d4,:,:,:]
     68    dimension = np.array(output).ndim
     69    shape = np.array(output).shape
     70    print 'dim,shape: ',dimension,shape
     71    return output, error
     72
    173def latinterv (area):
    274        if   area == "Europe":
    … …  
    35107        return wlon,wlat
    36108
     109def definesubplot ( numplot, fig ):
     110    from matplotlib.pyplot import rcParams
     111    rcParams['font.size'] = 12. ## default (important for multiple calls)
     112    if   numplot == 4:
     113        sub = 221
     114        fig.subplots_adjust(wspace = 0.3, hspace = 0.3)
     115        rcParams['font.size'] = int( rcParams['font.size'] * 2. / 3. )
     116    elif numplot == 2:
     117        sub = 121
     118        fig.subplots_adjust(wspace = 0.35)
     119        rcParams['font.size'] = int( rcParams['font.size'] * 3. / 4. )
     120    elif numplot == 3:
     121        sub = 131
     122        fig.subplots_adjust(wspace = 0.5)
     123        rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
     124    elif numplot == 6:
     125        sub = 231
     126        fig.subplots_adjust(wspace = 0.4, hspace = 0.0)
     127        rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
     128    elif numplot == 8:
     129        sub = 331 #241
     130        fig.subplots_adjust(wspace = 0.3, hspace = 0.3)
     131        rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
     132    elif numplot == 9:
     133        sub = 331
     134        fig.subplots_adjust(wspace = 0.3, hspace = 0.3)
     135        rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
     136    elif numplot == 1:
     137        sub = 99999
     138    elif numplot < 0:
     139        sub = 99999
     140    else:
     141        print "supported: 1,2,3,4,6,8,9"
     142        exit()
     143    return sub
     144
     145def getstralt(nc,nvert):
     146    typefile = whatkindfile(nc)
     147    if typefile is 'meso':                     
     148        stralt = "_lvl" + str(nvert)
     149    elif typefile is 'mesoapi':
     150        zelevel = int(nc.variables['vert'][nvert])
     151        if abs(zelevel) < 10000.:   strheight=str(zelevel)+"m"
     152        else:                       strheight=str(int(zelevel/1000.))+"km"
     153        if 'altitude'       in nc.dimensions:   stralt = "_"+strheight+"-AMR"
     154        elif 'altitude_abg' in nc.dimensions:   stralt = "_"+strheight+"-ALS"
     155        elif 'bottom_top'   in nc.dimensions:   stralt = "_"+strheight
     156        elif 'pressure'     in nc.dimensions:   stralt = "_"+str(zelevel)+"Pa"
     157        else:                                   stralt = ""
     158    else:
     159        stralt = ""
     160    return stralt
     161
    37162def getlschar ( namefile ):
    38     #### strangely enough this does not work for api or ncrcat results!
    39163    from netCDF4 import Dataset
    40164    from timestuff import sol2ls
     165    from numpy import array
    41166    nc  = Dataset(namefile)
    42     if namefile[0]+namefile[1]+namefile[2] != "geo" \
    43        and 'Times' in nc.variables \
     167    if 'Times' in nc.variables:
     168        zetime = nc.variables['Times'][0]
     169        shape = array(nc.variables['Times']).shape
     170        if shape[0] < 2: zetime = None
     171    if zetime is not None \
    44172       and 'vert' not in nc.variables:
    45         zetime = nc.variables['Times'][0]
     173       #### strangely enough this does not work for api or ncrcat results!
    46174        zetimestart = getattr(nc, 'START_DATE')
    47175        zeday = int(zetime[8]+zetime[9]) - int(zetimestart[8]+zetimestart[9])
    … …  
    67195
    68196def getproj (nc):
    69     map_proj = getattr(nc, 'MAP_PROJ')
    70     cen_lat  = getattr(nc, 'CEN_LAT')
    71     if map_proj == 2:
    72         if cen_lat > 10.:   
    73             proj="npstere"
    74             print "NP stereographic polar domain"
    75         else:           
    76             proj="spstere"
    77             print "SP stereographic polar domain"
    78     elif map_proj == 1:
    79         print "lambert projection domain"
    80         proj="lcc"
    81     elif map_proj == 3:
    82         print "mercator projection"
    83         proj="merc"
    84     else:
    85         proj="merc"
     197    typefile = whatkindfile(nc)
     198    if typefile in ['mesoapi','meso','geo']:
     199        ### (il faudrait passer CEN_LON dans la projection ?)
     200        map_proj = getattr(nc, 'MAP_PROJ')
     201        cen_lat  = getattr(nc, 'CEN_LAT')
     202        if map_proj == 2:
     203            if cen_lat > 10.:   
     204                proj="npstere"
     205                print "NP stereographic polar domain"
     206            else:           
     207                proj="spstere"
     208                print "SP stereographic polar domain"
     209        elif map_proj == 1:
     210            print "lambert projection domain"
     211            proj="lcc"
     212        elif map_proj == 3:
     213            print "mercator projection"
     214            proj="merc"
     215        else:
     216            proj="merc"
     217    elif typefile in ['gcm']:  proj="cyl"  ## pb avec les autres (de trace derriere la sphere ?)
     218    else:                      proj="ortho"
    86219    return proj   
    87220
    … …  
    102235    lat1 = lat2d[0,0]
    103236    lat2 = lat2d[nx,ny]
    104     wider = 0.5 * (abs(lon1)+abs(lon2)) * 0.1
    105     if lon1 < lon2:  wlon = [lon1, lon2 + wider]  ## a tester en normal
     237    if abs(0.5*(lat1+lat2)) > 60.:  wider = 0.5 * (abs(lon1)+abs(lon2)) * 0.1
     238    else:                           wider = 0.
     239    if lon1 < lon2:  wlon = [lon1, lon2 + wider] 
    106240    else:            wlon = [lon2, lon1 + wider]
    107241    if lat1 < lat2:  wlat = [lat1, lat2]
    … …  
    118252    return
    119253
    120 def makeplotpng (filename,pad_inches_value=0.25,minres=100.,folder=''):
    121     makeplotpngres(filename,minres,     pad_inches_value=pad_inches_value,folder=folder)
     254def makeplotpng (filename,pad_inches_value=0.25,minres=100.,folder='',disp=True):
     255    makeplotpngres(filename,minres,     pad_inches_value=pad_inches_value,folder=folder,disp=disp)
    122256    makeplotpngres(filename,minres+200.,pad_inches_value=pad_inches_value,folder=folder,disp=False)
    123257    return
    124258
    125 def dumpbdy (field):
     259def dumpbdy (field,stag=None):
    126260    nx = len(field[0,:])-1
    127261    ny = len(field[:,0])-1
     262    if stag == 'U': nx = nx-1
     263    if stag == 'V': ny = ny-1
    128264    return field[5:ny-5,5:nx-5]
     265
     266def getcoorddef ( nc ):   
     267    ## getcoord2d for predefined types
     268    typefile = whatkindfile(nc)
     269    if typefile in ['mesoapi','meso']:
     270        [lon2d,lat2d] = getcoord2d(nc)
     271        lon2d = dumpbdy(lon2d)
     272        lat2d = dumpbdy(lat2d)
     273    elif typefile in ['gcm']:
     274        [lon2d,lat2d] = getcoord2d(nc,nlat="latitude",nlon="longitude",is1d=True)
     275    elif typefile in ['geo']:
     276        [lon2d,lat2d] = getcoord2d(nc,nlat='XLAT_M',nlon='XLONG_M')
     277    return lon2d,lat2d   
    129278
    130279def getcoord2d (nc,nlat='XLAT',nlon='XLONG',is1d=False):
    … …  
    168317        return improc
    169318
    170 def getwinds (nc,charu='Um',charv='Vm'):
    171     import numpy as np
    172     u = nc.variables[charu]
    173     v = nc.variables[charv]
    174     if charu == 'U': u = u[:, :, :, 0:len(u[0,0,0,:])-1]
    175     if charv == 'V': v = v[:, :, 0:len(v[0,0,:,0])-1, :]
    176                      ### ou alors prendre les coordonnees speciales
    177     return u,v
     319def getwinddef (nc):   
     320    ## getwinds for predefined types
     321    typefile = whatkindfile(nc)
     322    ###
     323    if typefile is 'mesoapi':    [uchar,vchar] = ['Um','Vm']
     324    elif typefile is 'gcm':      [uchar,vchar] = ['u','v']
     325    elif typefile is 'meso':     [uchar,vchar] = ['U','V']
     326    else:                        [uchar,vchar] = ['not found','not found']
     327    ###
     328    if typefile in ['meso']:     metwind = False ## geometrical (wrt grid)
     329    else:                        metwind = True  ## meteorological (zon/mer)
     330    if metwind is False:         print "Not using meteorological winds. You trust numerical grid as being (x,y)"
     331    ###
     332    return uchar,vchar,metwind
    178333
    179334def vectorfield (u, v, x, y, stride=3, scale=15., factor=250., color='black', csmooth=1, key=True):
    … …  
    216371    return step
    217372
    218 def define_proj (char,wlon,wlat,back="."):
     373def define_proj (char,wlon,wlat,back=None):
    219374    from    mpl_toolkits.basemap            import Basemap
    220375    import  numpy                           as np
    … …  
    247402    m.drawmeridians(np.r_[-180.:180.:step*2.], labels=[0,0,0,1], color='grey', fontsize=fontsizemer)
    248403    m.drawparallels(np.r_[-90.:90.:step], labels=[1,0,0,0], color='grey', fontsize=fontsizemer)
    249     if back == ".":      m.warpimage(marsmap(),scale=0.75)
    250     elif back == None:   pass
    251     else:                m.warpimage(marsmap(back),scale=0.75)
     404    if back: m.warpimage(marsmap(back),scale=0.75)
     405            #if not back:
     406            #    if not var:                                        back = "mola"    ## if no var:         draw mola
     407            #    elif typefile in ['mesoapi','meso','geo'] \
     408            #       and proj not in ['merc','lcc','nsper','laea']:  back = "molabw"  ## if var but meso:   draw molabw
     409            #    else:                                              pass             ## else:              draw None
    252410    return m
    253411
    … …  
    270428    return
    271429
     430def calculate_bounds(field,vmin=None,vmax=None):
     431    import numpy as np
     432    from mymath import max,min,mean
     433    ind = np.where(field < 9e+35)
     434    fieldcalc = field[ ind ] # la syntaxe compacte ne marche si field est un tuple
     435    ###
     436    dev = np.std(fieldcalc)*3.0
     437    ###
     438    if vmin is None:
     439        zevmin = mean(fieldcalc) - dev
     440    else:             zevmin = vmin
     441    ###
     442    if vmax is None:  zevmax = mean(fieldcalc) + dev
     443    else:             zevmax = vmax
     444    if vmin == vmax:
     445                      zevmin = mean(fieldcalc) - dev  ### for continuity
     446                      zevmax = mean(fieldcalc) + dev  ### for continuity           
     447    ###
     448    if zevmin < 0. and min(fieldcalc) > 0.: zevmin = 0.
     449    print "field ", min(fieldcalc), max(fieldcalc)
     450    print "bounds ", zevmin, zevmax
     451    return zevmin, zevmax
     452
     453def bounds(what_I_plot,zevmin,zevmax):
     454    ### might be convenient to add the missing value in arguments
     455    what_I_plot[ what_I_plot < zevmin ] = zevmin*(1. + 1.e-7)
     456    what_I_plot[ what_I_plot > 9e+35  ] = -9e+35
     457    what_I_plot[ what_I_plot > zevmax ] = zevmax*(1. - 1.e-7)
     458    return what_I_plot
     459
     460def zoomset (wlon,wlat,zoom):
     461    dlon = abs(wlon[1]-wlon[0])/2.
     462    dlat = abs(wlat[1]-wlat[0])/2.
     463    [wlon,wlat] = [ [wlon[0]+zoom*dlon/100.,wlon[1]-zoom*dlon/100.],\
     464                    [wlat[0]+zoom*dlat/100.,wlat[1]-zoom*dlat/100.] ]
     465    print "zoom %",zoom,wlon,wlat
     466    return wlon,wlat
    272467
    273468def fmtvar (whichvar="def"):
    … …  
    286481    return fmtvar[whichvar]
    287482
     483####################################################################################################################
     484### Colorbars http://www.scipy.org/Cookbook/Matplotlib/Show_colormaps?action=AttachFile&do=get&target=colormaps3.png
    288485def defcolorb (whichone="def"):
    289486    whichcolorb =    { \
    … …  
    319516
    320517def marsmap (whichone="vishires"):
     518        from os import uname
     519        mymachine = uname()[1]
     520        ### not sure about speed-up with this method... looks the same
     521        if "lmd.jussieu.fr" in mymachine: domain = "/u/aslmd/WWW/maps/"
     522        else:                             domain = "http://www.lmd.jussieu.fr/~aslmd/maps/"
    321523        whichlink =     { \
    322                 "vis":          "http://maps.jpl.nasa.gov/pix/mar0kuu2.jpg",\
    323                 "vishires":     "http://dl.dropbox.com/u/11078310/MarsMap_2500x1250.jpg",\
    324                 "geolocal":     "http://dl.dropbox.com/u/11078310/geolocal.jpg",\
    325                 "mola":         "http://www.lns.cornell.edu/~seb/celestia/mars-mola-2k.jpg",\
    326                 "molabw":       "http://dl.dropbox.com/u/11078310/MarsElevation_2500x1250.jpg",\
     524                #"vis":         "http://maps.jpl.nasa.gov/pix/mar0kuu2.jpg",\
     525                #"vishires":    "http://www.lmd.jussieu.fr/~aslmd/maps/MarsMap_2500x1250.jpg",\
     526                #"geolocal":    "http://dl.dropbox.com/u/11078310/geolocal.jpg",\
     527                #"mola":        "http://www.lns.cornell.edu/~seb/celestia/mars-mola-2k.jpg",\
     528                #"molabw":      "http://dl.dropbox.com/u/11078310/MarsElevation_2500x1250.jpg",\
     529                "vis":         domain+"mar0kuu2.jpg",\
     530                "vishires":    domain+"MarsMap_2500x1250.jpg",\
     531                "geolocal":    domain+"geolocal.jpg",\
     532                "mola":        domain+"mars-mola-2k.jpg",\
     533                "molabw":      domain+"MarsElevation_2500x1250.jpg",\
    327534                        }
    328535        if whichone not in whichlink:

  • trunk/MESOSCALE_DEV/PLOT/PYTHON/scripts/makefig.sh

    r225 r233  
    55 file2="$fold/POLAR_THOMAS_ls20/wrfout_d03_2024-01-42_06:00:00"
    66 file3="$fold/POLAR_THOMAS_ls0/wrfout_d03_2024-01-06_06:00:00"
     7 file4="$fold/POLAR_THOMAS_ls10/wrfout_d03_2024-01-19_06:00:00"
    78 dest="/u/aslmd/WWW/antichambre/thomas/"
    89
    … …  
    1213           -v USTM -m   0. -M 0.8 \
    1314           -v HGT  -m   0  -M 0 \
    14            -f $file1 -f $file2 -f $file3
     15           -f $file1 -f $file2 -f $file3 -f $file4 \
     16           -d

  • trunk/MESOSCALE_DEV/PLOT/PYTHON/scripts/winds.py

    r232 r233  
    2323           vmax=None,\
    2424           tile=False,\
    25            zoom=None):
     25           zoom=None,\
     26           display=True,\
     27           itstep=None):
    2628
    2729    ####################################################################################################################
    … …  
    3234    from netCDF4 import Dataset
    3335    from myplot import getcoord2d,define_proj,makeplotpng,simplinterv,vectorfield,ptitle,latinterv,getproj,wrfinterv,dumpbdy,\
    34                        fmtvar,definecolorvec,getwinds,defcolorb,getprefix,putpoints
     36                       fmtvar,definecolorvec,defcolorb,getprefix,putpoints,calculate_bounds,errormess,definesubplot,\
     37                       zoomset,getcoorddef,getwinddef,whatkindfile,reducefield,bounds,getstralt,getfield
    3538    from mymath import deg,max,min,mean
    3639    from matplotlib.pyplot import contour,contourf, subplot, figure, rcParams, savefig, colorbar, pcolor
    … …  
    3942    from numpy.core.defchararray import find
    4043
    41     ###
    42     #rcParams['text.usetex'] = True
    43     #rcParams['cairo.format'] = 'svg'
    44 
    4544    ######################
    4645    ### Load NETCDF object
    4746    nc  = Dataset(namefile)
    48 
    49     ###################################
    50     ### Recognize predefined file types
    51     if 'controle' in nc.variables:   typefile = 'gcm'
    52     elif 'vert' in nc.variables:     typefile = 'mesoapi'
    53     elif 'U' in nc.variables:        typefile = 'meso'
    54     elif 'HGT_M' in nc.variables:    typefile = 'geo'
    55     else:                           
    56         print "typefile not supported."
    57         print nc.variables
    58         exit()
    59 
    60     ##############################################################
    61     ### Try to guess the projection from wrfout if not set by user
    62     if typefile in ['mesoapi','meso','geo']:
    63         if proj == None:       proj = getproj(nc)
    64                                     ### (il faudrait passer CEN_LON dans la projection ?)
    65     elif typefile in ['gcm']:
    66         if proj == None:       proj = "cyl"   
    67                                     ## pb avec les autres (de trace derriere la sphere ?)
    68 
    69     ############################################
    70     #### Choose underlying topography by default
    71     if not back:
    72         if not var:                                        back = "mola"    ## if no var:         draw mola
    73         elif typefile in ['mesoapi','meso','geo'] \
    74            and proj not in ['merc','lcc','nsper','laea']:  back = "molabw"  ## if var but meso:   draw molabw
    75         else:                                              pass             ## else:              draw None
    76 
    77     ####################################################
    78     ### Get geographical coordinates and plot boundaries
    79     if typefile in ['mesoapi','meso']:
    80         [lon2d,lat2d] = getcoord2d(nc)
    81         lon2d = dumpbdy(lon2d)
    82         lat2d = dumpbdy(lat2d)
    83     elif typefile in ['gcm']:
    84         [lon2d,lat2d] = getcoord2d(nc,nlat="latitude",nlon="longitude",is1d=True)
    85     elif typefile in ['geo']:
    86         [lon2d,lat2d] = getcoord2d(nc,nlat='XLAT_M',nlon='XLONG_M')
     47   
     48    ##################################
     49    ### Initial checks and definitions
     50    typefile = whatkindfile(nc)                                  ## TYPEFILE
     51    if var not in nc.variables: var = False                      ## VAR
     52    if winds: [uchar,vchar,metwind] = getwinddef(nc)             ## WINDS
     53    if uchar == 'not found': winds = False
     54    [lon2d,lat2d] = getcoorddef(nc)                              ## COORDINATES, could be moved below
     55    if proj == None:   proj = getproj(nc)                        ## PROJECTION
     56
     57    ##########################
     58    ### Define plot boundaries
    8759    if proj == "npstere":             [wlon,wlat] = latinterv("North_Pole")
    8860    elif proj in ["lcc","laea"]:      [wlon,wlat] = wrfinterv(lon2d,lat2d)
    8961    else:                             [wlon,wlat] = simplinterv(lon2d,lat2d)
    90     print wlon
    91     print wlat
    92     if zoom: 
    93         dlon = abs(wlon[1]-wlon[0])/2.
    94         dlat = abs(wlat[1]-wlat[0])/2.
    95         [wlon,wlat] = [ [wlon[0]+zoom*dlon/100.,wlon[1]-zoom*dlon/100.],\
    96                         [wlat[0]+zoom*dlat/100.,wlat[1]-zoom*dlat/100.] ]
    97         print "zoom %",zoom,wlon,wlat
    98 
    99     ##############################################################################
    100     ### Get winds and know if those are meteorological winds (ie. zon, mer) or not
    101     if winds:
    102         if typefile is 'mesoapi':
    103             [u,v] = getwinds(nc)
    104             metwind = True  ## meteorological (zon/mer)
    105         elif typefile is 'gcm':
    106             [u,v] = getwinds(nc,charu='u',charv='v')
    107             metwind = True  ## meteorological (zon/mer)
    108         elif typefile is 'meso':
    109             [u,v] = getwinds(nc,charu='U',charv='V')
    110             metwind = False ## geometrical (wrt grid)
    111             print "Beware ! Not using meteorological winds. You trust numerical grid as being (x,y)."
    112         elif typefile is 'geo':
    113             winds = None
    114 
    115     #####################################################
    116     ### Load the chosen variables, whether it is 2D or 3D
    117     if var:
    118         if var not in nc.variables:
    119             print "not found in file:",var
    120             exit()
    121         else:   
    122             dimension = np.array(nc.variables[var]).ndim
    123             if dimension == 2:     field = nc.variables[var][:,:]
    124             elif dimension == 3:   field = nc.variables[var][:,:,:]
    125             elif dimension == 4:   field = nc.variables[var][:,nvert,:,:] 
    126         fieldcalc = field[ field < 9e+35 ]
    127         dev = np.std(fieldcalc)*2.0
    128         if vmin is None:  zevmin = mean(fieldcalc) - dev
    129         else:             zevmin = vmin
    130         if vmax is None:  zevmax = mean(fieldcalc) + dev
    131         else:             zevmax = vmax
    132         if vmin == vmax:   
    133                           zevmin = mean(fieldcalc) - dev  ### for continuity
    134                           zevmax = mean(fieldcalc) + dev  ### for continuity           
    135         print "field ", min(fieldcalc), max(fieldcalc)
    136         print "bounds ", zevmin, zevmax
    137         ### some already defined colormaps
    138         if colorb is True:    colorb = defcolorb(var)
    139     else:
    140         dimension = 0
    141  
    142     ###########################
    143     ### Get length of time axis
    144     if winds:               nt = len(u[:,0,0,0])
    145     elif var: 
    146         if dimension == 2:  nt = 1
    147         else             :  nt = len(field[:,0,0])
     62    if zoom:                          [wlon,wlat] = zoomset(wlon,wlat,zoom)
    14863
    14964    #########################################
    … …  
    15267    elif var:             basename = var
    15368    elif winds:           basename = 'UV'
    154     else:                 exit()
    155     ###
    156     if dimension == 4 or winds:
    157         if typefile is 'meso':                      stralt = "_lvl" + str(nvert)
    158         elif typefile is 'mesoapi': 
    159             zelevel = int(nc.variables['vert'][nvert])
    160             if abs(zelevel) < 10000.:   strheight=str(zelevel)+"m"
    161             else:                       strheight=str(int(zelevel/1000.))+"km"
    162             if 'altitude'       in nc.dimensions:   stralt = "_"+strheight+"-AMR"
    163             elif 'altitude_abg' in nc.dimensions:   stralt = "_"+strheight+"-ALS"
    164             elif 'bottom_top'   in nc.dimensions:   stralt = "_"+strheight
    165             elif 'pressure'     in nc.dimensions:   stralt = "_"+str(zelevel)+"Pa"
    166         else:                                       stralt = ""         
    16769    else:
    168         stralt = ""
    169     ###
    170     basename = basename + stralt
     70        print nc.variables                 
     71        errormess("please set at least winds or var")
     72    basename = basename + getstralt(nc,nvert)  ## can be moved elsewhere for a more generic routine
    17173
    17274    ##################################
    17375    ### Open a figure and set subplots
    17476    fig = figure()
    175     rcParams['font.size'] = 12.
    176     if   typefile in ['geo']:  numplot = 1
    177     if   numplot > 0:   
    178         if   numplot == 4:
    179             sub = 221
    180             fig.subplots_adjust(wspace = 0.3, hspace = 0.3)
    181             rcParams['font.size'] = int( rcParams['font.size'] * 2. / 3. )
    182         elif numplot == 2:
    183             sub = 121
    184             fig.subplots_adjust(wspace = 0.35)
    185             rcParams['font.size'] = int( rcParams['font.size'] * 3. / 4. )
    186         elif numplot == 3:
    187             sub = 131
    188             fig.subplots_adjust(wspace = 0.5)
    189             rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
    190         elif numplot == 6:
    191             sub = 231
    192             fig.subplots_adjust(wspace = 0.4, hspace = 0.0)
    193             rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
    194         elif numplot == 8:
    195             sub = 331 #241
    196             fig.subplots_adjust(wspace = 0.3, hspace = 0.3)
    197             rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
    198         elif numplot == 9:
    199             sub = 331
    200             fig.subplots_adjust(wspace = 0.3, hspace = 0.3)
    201             rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
    202         elif numplot == 1:
    203             sub = 99999
    204         else:
    205             print "supported: 1,2,3,4,6,8,9"
    206             exit()
    207         ### Prepare time loop
    208         if nt <= numplot or numplot == 1: 
    209             tabrange = [0]
    210             numplot = 1
    211         else:                         
    212             tabrange = range(0,nt,int(nt/numplot))  #nt-1
    213             tabrange = tabrange[0:numplot]
    214     else:
    215         tabrange = range(0,nt,1)
    216         sub = 99999
    217     print tabrange
    218 
     77    sub = definesubplot( numplot, fig )
     78 
    21979    #################################
    22080    ### Time loop for plotting device
    22181    found_lct = False
    222     for i in tabrange:
     82    itime = 0  ## could be an argument
     83    nplot = 1
     84    error = False
     85    if itstep is None: itstep = int(24./numplot)
     86    while error is False:
    22387
    22488       ### Which local time ?
    225        ltst = ( interv[0] + 0.5*(wlon[0]+wlon[1])/15.) + i*interv[1]
     89       print interv[0], interv[1], itime
     90       ltst = ( interv[0] + 0.5*(wlon[0]+wlon[1])/15.) + itime*interv[1]
    22691       ltst = int (ltst * 10) / 10.
    22792       ltst = ltst % 24
    22893
    22994       ### General plot settings
    230        if numplot > 1:
    231            subplot(sub)
     95       #print itime, int(ltst), numplot, nplot
     96       if numplot >= 1:
     97           if nplot > numplot: break
     98           if numplot > 1:     
     99               if typefile not in ['geo']:  subplot(sub+nplot-1)
     100
    232101           found_lct = True
    233        elif numplot == 1:
    234            found_lct = True
    235             ### If only one local time is requested (numplot < 0)
     102       ### If only one local time is requested (numplot < 0)
    236103       elif numplot <= 0:
    237            if int(ltst) + numplot != 0:         continue
    238            else:                                found_lct = True
     104           if int(ltst) + numplot != 0:         
     105                 itime += 1
     106                 if found_lct is True: break     ## because it means LT was found at previous iteration
     107                 else:                 continue  ## continue to iterate to find the correct LT
     108           else:
     109                 found_lct = True
    239110
    240111       ### Map projection
    … …  
    244115       #### Contour plot
    245116       if var:
    246            if typefile in ['mesoapi','meso']:    what_I_plot = dumpbdy(field[i,:,:])
    247            elif typefile in ['geo']:             what_I_plot = field[0,:,:]
    248            elif typefile in ['gcm']:             
    249                if dimension == 2:                what_I_plot = field[:,:]
    250                elif dimension == 3:              what_I_plot = field[i,:,:]
    251            palette = get_cmap(name=colorb)
    252            #palette.set_over('b', 1.0)
    253            #print np.array(x).shape
    254            #print np.array(y).shape
    255            #print np.array(what_I_plot).shape
    256            if not tile:
    257                zelevels = np.linspace(zevmin,zevmax)
    258                hole = True
    259                if not hole:
    260                    what_I_plot[ what_I_plot < zevmin ] = zevmin*(1. + 1.e-7)
    261                    what_I_plot[ what_I_plot > 9e+35  ] = -9e+35
    262                    what_I_plot[ what_I_plot > zevmax ] = zevmax*(1. - 1.e-7)
    263                contourf( x, y, what_I_plot, 10, cmap = palette, levels = zelevels )
    264            else:   
    265                pcolor( x, y, what_I_plot, cmap = palette, vmin=zevmin, vmax=zevmax )
    266            #putpoints(m,fig)
    267            if var in ['HGT']: pass
    268            elif colorb:             
    269                               ndiv = 10
    270                               colorbar(fraction=0.05,pad=0.1,format=fmtvar(var),\
    271                                        ticks=np.linspace(zevmin,zevmax,ndiv+1),\
    272                                        extend='max',spacing='proportional')
    273                                        # both min max neither
     117           what_I_plot, error = reducefield( getfield(nc,var), d4=itime, d3=nvert )
     118           if not error:
     119               if typefile in ['mesoapi','meso']:    what_I_plot = dumpbdy(what_I_plot)
     120               zevmin, zevmax = calculate_bounds(what_I_plot,vmin=vmin,vmax=vmax)
     121               if colorb is True: colorb = defcolorb(var)
     122               palette = get_cmap(name=colorb)
     123               if not tile:
     124                   hole = True
     125                   hole = None
     126                   if not hole:  what_I_plot = bounds(what_I_plot,zevmin,zevmax)
     127                   zelevels = np.linspace(zevmin,zevmax)
     128                   contourf( x, y, what_I_plot, zelevels, cmap = palette )
     129               else:   
     130                   pcolor( x, y, what_I_plot, cmap = palette, vmin=zevmin, vmax=zevmax )
     131               #putpoints(m,fig)
     132               if var in ['HGT']: pass
     133               elif colorb:             
     134                         ndiv = 10
     135                         colorbar(fraction=0.05,pad=0.1,format=fmtvar(var),\
     136                                           ticks=np.linspace(zevmin,zevmax,ndiv+1),\
     137                                           extend='max',spacing='proportional')
     138                                           # both min max neither
     139                 
    274140       ### Vector plot
    275141       if winds:
    276            if   typefile in ['mesoapi','meso']:   
    277                [vecx,vecy] = [dumpbdy(u[i,nvert,:,:]), dumpbdy(v[i,nvert,:,:])]
    278                key = True
    279            elif typefile in ['gcm']:               
    280                [vecx,vecy] = [        u[i,nvert,:,:] ,         v[i,nvert,:,:] ]
    281                key = False
    282            if metwind:  [vecx,vecy] = m.rotate_vector(vecx, vecy, lon2d, lat2d)
    283            if var == None:        colorvec = definecolorvec(back)
    284            else:                  colorvec = definecolorvec(colorb)
    285            vectorfield(vecx, vecy,\
    286                       x, y, stride=stride, csmooth=2,\
    287                       scale=15., factor=300., color=colorvec, key=key)
    288                                         #200.         ## or csmooth=stride
    289        
     142           vecx, error = reducefield( getfield(nc,uchar), d4=itime, d3=nvert )
     143           vecy, error = reducefield( getfield(nc,vchar), d4=itime, d3=nvert )
     144           if not error:
     145               if typefile in ['mesoapi','meso']:   
     146                   [vecx,vecy] = [dumpbdy(vecx,stag=uchar), dumpbdy(vecy,stag=vchar)]
     147                   key = True
     148               elif typefile in ['gcm']:           
     149                   key = False
     150               if metwind:  [vecx,vecy] = m.rotate_vector(vecx, vecy, lon2d, lat2d)
     151               if var == False:       colorvec = definecolorvec(back)
     152               else:                  colorvec = definecolorvec(colorb)
     153               vectorfield(vecx, vecy,\
     154                          x, y, stride=stride, csmooth=2,\
     155                          scale=15., factor=300., color=colorvec, key=key)
     156                                            #200.         ## or csmooth=stride
     157               
    290158       ### Next subplot
    291159       plottitle = basename
    292160       if typefile in ['mesoapi','meso']:
    293            if addchar:  plottitle = plottitle + addchar + "_LT"+str(ltst)
    294            else:        plottitle = plottitle + "_LT"+str(ltst)
     161            if addchar:  plottitle = plottitle + addchar + "_LT"+str(ltst)
     162            else:        plottitle = plottitle + "_LT"+str(ltst)
    295163       ptitle( plottitle )
    296        sub += 1
     164       itime += itstep
     165       nplot += 1
    297166
    298167    ##########################################################################
    … …  
    309178    else:               zeplot = target + "/" + zeplot 
    310179    ###
    311     if found_lct:     makeplotpng(zeplot,pad_inches_value=0.35)   
     180    if found_lct:     makeplotpng(zeplot,pad_inches_value=0.35,disp=display)   
    312181    else:             print "Local time not found"
    313182
    … …  
    315184    ### Now the end
    316185    return zeplot
     186
     187##############################
     188### A specific stuff for below
     189def adjust_length (tab, zelen):
     190    from numpy import ones
     191    if tab is None:
     192        outtab = ones(zelen) * -999999
     193    else:
     194        print zelen, len(tab)
     195        if zelen != len(tab):
     196            print "not enough or too much values... setting same values all variables"
     197            outtab = ones(zelen) * tab[0]
     198        else:
     199            outtab = tab
     200    return outtab
    317201
    318202###########################################################################################
    … …  
    325209    from netCDF4 import Dataset
    326210    from myplot import getlschar
    327     from numpy import ones
    328211
    329212    #############################
    330213    ### Get options and variables
    331214    parser = OptionParser()
    332     parser.add_option('-f', action='append', dest='namefile',    type="string",  default=None,  help='[NEEDED] name of WRF file (append)')
    333     parser.add_option('-l', action='store', dest='nvert',       type="float",   default=0,     help='vertical level (def=0)(-i 2: p,mbar)(-i 3,4: z,km)')
    334     parser.add_option('-p', action='store', dest='proj',        type="string",  default=None,  help='projection')
    335     parser.add_option('-b', action='store', dest='back',        type="string",  default=None,  help='background')
    336     parser.add_option('-t', action='store', dest='target',      type="string",  default=None,  help='destination folder')
    337     parser.add_option('-s', action='store', dest='stride',      type="int",     default=3,     help='stride vectors (def=3)')
    338     parser.add_option('-v', action='append', dest='var',         type="string",  default=None,  help='variable contoured (append)')
    339     parser.add_option('-n', action='store', dest='numplot',     type="int",     default=2,     help='number of plots (def=1)(<0: 1 plot of LT -*numplot*)')
    340     parser.add_option('-i', action='store', dest='interp',      type="int",     default=None,  help='interpolation method (2: press, 3: z-amr, 4:z-als)')
    341     parser.add_option('-c', action='store', dest='colorb',      type="string",  default=True,  help='change colormap')
    342     parser.add_option('-x', action='store_false', dest='winds',                 default=True,  help='no wind vectors')
    343     parser.add_option('-m', action='append', dest='vmin',        type="float",   default=None,  help='bounding minimum value (append)')   
    344     parser.add_option('-M', action='append', dest='vmax',        type="float",   default=None,  help='bounding maximum value (append)')
    345     parser.add_option('-T', action='store_true', dest='tile',                   default=False, help='draw a tiled plot (no blank zone)')
    346     parser.add_option('-z', action='store', dest='zoom',        type="float",   default=None,  help='zoom factor in %')
    347     parser.add_option('-N', action='store_true', dest='nocall',                 default=False, help='do not recreate api file')
     215    parser.add_option('-f', action='append',dest='namefile',    type="string",  default=None,  help='[NEEDED] name of WRF file (append)')
     216    parser.add_option('-l', action='store',dest='nvert',        type="float",   default=0,     help='vertical level (def=0)(-i 2: p,mbar)(-i 3,4: z,km)')
     217    parser.add_option('-p', action='store',dest='proj',         type="string",  default=None,  help='projection')
     218    parser.add_option('-b', action='store',dest='back',         type="string",  default=None,  help='background image (def: None)')
     219    parser.add_option('-t', action='store',dest='target',       type="string",  default=None,  help='destination folder')
     220    parser.add_option('-s', action='store',dest='stride',       type="int",     default=3,     help='stride vectors (def=3)')
     221    parser.add_option('-v', action='append',dest='var',         type="string",  default=None,  help='variable contoured (append)')
     222    parser.add_option('-n', action='store',dest='numplot',      type="int",     default=2,     help='number of plots (def=1)(<0: 1 plot of LT -*numplot*)')
     223    parser.add_option('-i', action='store',dest='interp',       type="int",     default=None,  help='interpolation method (2: press, 3: z-amr, 4:z-als)')
     224    parser.add_option('-c', action='store',dest='colorb',       type="string",  default=True,  help='change colormap')
     225    parser.add_option('-x', action='store_false',dest='winds',                  default=True,  help='no wind vectors')
     226    parser.add_option('-m', action='append',dest='vmin',        type="float",   default=None,  help='bounding minimum value (append)')   
     227    parser.add_option('-M', action='append',dest='vmax',        type="float",   default=None,  help='bounding maximum value (append)')
     228    parser.add_option('-T', action='store_true',dest='tile',                    default=False, help='draw a tiled plot (no blank zone)')
     229    parser.add_option('-z', action='store',dest='zoom',         type="float",   default=None,  help='zoom factor in %')
     230    parser.add_option('-N', action='store_true',dest='nocall',                  default=False, help='do not recreate api file')
     231    parser.add_option('-d', action='store_false',dest='display',                default=True,  help='do not pop up created images')
     232    parser.add_option('-e', action='store',dest='itstep',       type="int",     default=None,  help='stride time (def=4)')
    348233    #parser.add_option('-V', action='store', dest='comb',        type="float",   default=None,  help='a defined combination of variables')
    349234    (opt,args) = parser.parse_args()
    … …  
    352237        exit()
    353238    print "Options:", opt
    354    
     239
     240    listvar = ''
    355241    if opt.var is None:
    356         pass
     242        zerange = [-999999]
    357243    else:
    358         listvar = ''
    359244        zelen = len(opt.var)
    360         if zelen == 1: listvar = opt.var[0] + ','
    361         else         :
    362             for jjj in range(zelen): listvar += opt.var[jjj] + ','
     245        zerange = range(zelen)
     246        #if zelen == 1: listvar = opt.var[0] + ','
     247        #else         :
     248        for jjj in zerange: listvar += opt.var[jjj] + ','
    363249        listvar = listvar[0:len(listvar)-1]
    364         if opt.vmin:
    365             if zelen != len(opt.vmin):
    366                 print "not enough or too much vmin values... setting same values all variables"
    367                 vmintab = ones(zelen) * opt.vmin[0]
    368             else:
    369                 vmintab = opt.vmin
    370         else:
    371             vmintab = None
    372         if opt.vmax:
    373             if zelen != len(opt.vmax):
    374                 print "not enough or too much vmax values... setting same values all variables"
    375                 vmaxtab = ones(zelen) * opt.vmax[0]
    376             else:
    377                 vmaxtab = opt.vmax
    378         else:
    379             vmaxtab = None
     250        vmintab = adjust_length (opt.vmin, zelen) 
     251        vmaxtab = adjust_length (opt.vmax, zelen)
    380252
    381253    for i in range(len(opt.namefile)):
    … …  
    395267            else                    :  zefields = ''
    396268            ### var or no var
    397             if opt.var is None      :  pass
    398             elif zefields == ''     :  zefields = listvar
     269            #if opt.var is None      :  pass
     270            if zefields == ''       :  zefields = listvar
    399271            else                    :  zefields = zefields + "," + listvar
    400272            print zefields
    … …  
    409281            zelevel = 0 ## so that zelevel could play again the role of nvert
    410282
    411         if opt.var is None: 
     283        if opt.var is None: zerange = [-999999]
     284        else:               zerange = range(zelen)
     285        for jjj in zerange:
     286            if jjj == -999999:
     287                argvar  = None
     288                argvmin = None
     289                argvmax = None
     290            else:
     291                argvar = opt.var[jjj]
     292                if vmintab[jjj] != -999999:  argvmin = vmintab[jjj]
     293                else:                        argvmin = None
     294                if vmaxtab[jjj] != -999999:  argvmax = vmaxtab[jjj]
     295                else:                        argvmax = None
    412296            #############
    413297            ### Main call
    414298            name = winds (zefile,int(zelevel),\
    415                proj=opt.proj,back=opt.back,target=opt.target,stride=opt.stride,var=opt.var,numplot=opt.numplot,colorb=opt.colorb,winds=opt.winds,\
    416                addchar=lschar,interv=[zehour,zehourin],vmin=opt.vmin,vmax=opt.vmax,tile=opt.tile,zoom=opt.zoom)
     299                proj=opt.proj,back=opt.back,target=opt.target,stride=opt.stride,var=argvar,\
     300                numplot=opt.numplot,colorb=opt.colorb,winds=opt.winds,\
     301                addchar=lschar,interv=[zehour,zehourin],vmin=argvmin,vmax=argvmax,\
     302                tile=opt.tile,zoom=opt.zoom,display=opt.display,\
     303                itstep=opt.itstep)
    417304            print 'Done: '+name
    418         else:       
    419             for jjj in range(len(opt.var)): 
    420                 if vmintab: argvmin = vmintab[jjj]
    421                 else:       argvmin = None
    422                 if vmaxtab: argvmax = vmaxtab[jjj]
    423                 else:       argvmax = None
    424                 #############
    425                 ### Main call
    426                 name = winds (zefile,int(zelevel),\
    427                    proj=opt.proj,back=opt.back,target=opt.target,stride=opt.stride,var=opt.var[jjj],numplot=opt.numplot,colorb=opt.colorb,winds=opt.winds,\
    428                    addchar=lschar,interv=[zehour,zehourin],vmin=argvmin,vmax=argvmax,tile=opt.tile,zoom=opt.zoom)
    429                 print 'Done: '+name
    430305   
    431306        #########################################################
Note: See TracChangeset for help on using the changeset viewer.