Ignore:
Timestamp:
Jul 1, 2011, 5:14:38 AM (14 years ago)
Author:
aslmd
Message:

MESOSCALE: much improved python scripts for winds. lots of automatic settings.

Location:
trunk/MESOSCALE/PLOT/PYTHON/scripts
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/MESOSCALE/PLOT/PYTHON/scripts/domain.py

    r181 r184  
    11#!/usr/bin/env python
    22
    3 ### A. Spiga LMD 29/05/2011
     3### A. Spiga -- LMD -- 30/05/2011
    44
    5 def usage():
    6     print 'USAGE:  plot.py file (target)'
    7     print 'file  : name of input netcdf file.'
    8     print '(target) : a directory with write rights.'
    9     print 'Example:  domain.py /d5/aslmd/LMD_MM_MARS_SIMUS/OM/OM6_TI85/wrfout_d01_2024-06-43_06:00:00 ~/'
    10 
    11 def domain (namefile,proj="ortho",back="molabw",target=None): #vishires
     5def domain (namefile,proj="ortho",back="molabw",target=None,var='HGT'):
    126    from netCDF4 import Dataset
    137    from myplot import getcoord2d,define_proj,makeplotpng,simplinterv
    … …  
    1812    m = define_proj(proj,wlon,wlat,back=back)
    1913    x, y = m(lon2d, lat2d)
    20     contourf(x, y, nc.variables['HGT'][0,:,:] / 1000., 50)
     14    contourf(x, y, nc.variables[var][0,:,:], 40)
    2115    if not target:   zeplot = namefile+".domain"
    22     else:            zeplot = target+"domain"
     16    else:            zeplot = target+"/domain"
    2317    makeplotpng(zeplot,pad_inches_value=0.35)
    2418
    2519if __name__ == "__main__":
    2620    import sys
    27     if (len(sys.argv)) == 2:     domain( str(sys.argv[1]) )
    28     elif (len(sys.argv)) == 3:   domain( str(sys.argv[1]) , target = str(sys.argv[2]) )
     21    ### to be replaced by argparse
     22    from optparse import OptionParser
     23    parser = OptionParser()
     24    parser.add_option('-f', action='store', dest='namefile',    type="string",  default=None,     help='name of WRF file [NEEDED]')
     25    parser.add_option('-p', action='store', dest='proj',        type="string",  default="ortho",  help='projection')
     26    parser.add_option('-b', action='store', dest='back',        type="string",  default="molabw", help='background')
     27    parser.add_option('-t', action='store', dest='target',      type="string",  default=None,     help='destination folder')
     28    parser.add_option('-v', action='store', dest='var',         type="string",  default='HGT',    help='variable contoured')
     29    (opt,args) = parser.parse_args()
     30    if opt.namefile is None:
     31        print "I want to eat one file at least ! Use domain.py -f name_of_my_file. Or type domain.py -h"
     32        exit()
     33    print "Options:", opt
     34    domain (opt.namefile,proj=opt.proj,back=opt.back,target=opt.target,var=opt.var)

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

    r183 r184  
    11#!/usr/bin/env python
    22
    3 ### A. Spiga and A. Colaitis -- LMD -- 30/05/2011
    4 
     3### A. Spiga -- LMD -- 30/05/2011
     4### Thanks to A. Colaitis for the parser trick
     5
     6
     7####################################
     8####################################
     9### The main program to plot vectors
    510def winds (namefile,\
    611           nvert,\
    7            proj="cyl",\
     12           proj=None,\
    813           back=None,\
    914           target=None,
    10            stride=5,\
    11            var='HGT'):
     15           stride=3,\
     16           numplot=4,\
     17           var=None):
     18
     19    #################################
     20    ### Load librairies and functions
    1221    from netCDF4 import Dataset
    13     from myplot import getcoord2d,define_proj,makeplotpng,simplinterv,vectorfield,ptitle
    14     from matplotlib.pyplot import contourf, subplot, figure
     22    from myplot import getcoord2d,define_proj,makeplotpng,simplinterv,vectorfield,ptitle,latinterv,getproj,wrfinterv,dumpbdy
     23    from matplotlib.pyplot import contourf, subplot, figure, rcParams, savefig
    1524    import numpy as np
     25
     26    #############################
     27    ### Lower a bit the font size
     28    rcParams['font.size'] = int( rcParams['font.size'] * 2. / 3. )
     29
     30    ######################
     31    ### Load NETCDF object
    1632    nc  = Dataset(namefile)
    17     [lon2d,lat2d] = getcoord2d(nc)
    18     [wlon,wlat] = simplinterv(lon2d,lat2d)
    19     [u,v] = getwinds(nc)
     33
     34    ###################################
     35    ### Recognize predefined file types
     36    if 'controle' in nc.variables:   typefile = 'gcm'
     37    elif 'Um' in nc.variables:       typefile = 'mesoapi'
     38    elif 'U' in nc.variables:        typefile = 'meso'
     39    else:                           
     40        print "typefile not supported."
     41        exit()
     42
     43    ##############################################################
     44    ### Try to guess the projection from wrfout if not set by user
     45    if typefile in ['mesoapi','meso']:
     46        if proj == None:       proj = getproj(nc)
     47                                    ### (il faudrait passer CEN_LON dans la projection ?)
     48    elif typefile in ['gcm']:
     49        if proj == None:       proj = "cyl"   
     50                                    ## pb avec les autres (de trace derriere la sphere ?)
     51
     52    ###################################################################
     53    ### For mesoscale results plot the underlying topography by default
     54    if typefile in ['mesoapi','meso']:
     55        if var == None: var = 'HGT'
     56
     57    ####################################################
     58    ### Get geographical coordinates and plot boundaries
     59    if typefile in ['mesoapi','meso']:
     60        [lon2d,lat2d] = getcoord2d(nc)
     61        lon2d = dumpbdy(lon2d)
     62        lat2d = dumpbdy(lat2d)
     63    elif typefile in ['gcm']:
     64        [lon2d,lat2d] = getcoord2d(nc,nlat="latitude",nlon="longitude",is1d=True)
     65    if proj == "npstere":    [wlon,wlat] = latinterv("North_Pole")
     66    elif proj == "lcc":      [wlon,wlat] = wrfinterv(lon2d,lat2d)
     67    else:                    [wlon,wlat] = simplinterv(lon2d,lat2d)
     68
     69    ##################
     70    ### Get local time
     71    if typefile in ['mesoapi','meso']:  ltst = int(getattr(nc, 'GMT') + 0.5*(wlon[0]+wlon[1])/15.)
     72    elif typefile in ['gcm']:           ltst = 0
     73
     74    ##############################################################################
     75    ### Get winds and know if those are meteorological winds (ie. zon, mer) or not
     76    if typefile is 'mesoapi':
     77        [u,v] = getwinds(nc)
     78        metwind = True  ## meteorological (zon/mer)
     79    elif typefile is 'gcm':
     80        [u,v] = getwinds(nc,charu='u',charv='v')
     81        metwind = True  ## meteorological (zon/mer)
     82    elif typefile is 'meso':
     83        [u,v] = getwinds(nc,charu='U',charv='V')
     84        metwind = False ## geometrical (wrt grid)
     85
     86    ####################################################
     87    ### Load the chosen variables, whether it is 2D or 3D
     88    if var:
     89        if var not in nc.variables:
     90            print "not found in file:",var
     91            exit()
     92        else:   
     93            dimension = np.array(nc.variables[var]).ndim
     94            if dimension == 2:     field = nc.variables[var][:,:]
     95            elif dimension == 3:   field = nc.variables[var][:,:,:]
     96            elif dimension == 4:   field = nc.variables[var][:,nvert,:,:] 
     97
     98    ##################################
     99    ### Open a figure and set subplots
     100    fig = figure()
     101    if   numplot == 4:
     102        sub = 221
     103        fig.subplots_adjust(wspace = 0.1, hspace = 0.3)
     104    elif numplot == 2:
     105        sub = 121
     106        fig.subplots_adjust(wspace = 0.3)
     107    elif numplot == 3:
     108        sub = 131
     109        fig.subplots_adjust(wspace = 0.3)
     110    elif numplot == 6:
     111        sub = 231
     112        fig.subplots_adjust(wspace = 0.4, hspace = 0.0)
     113    elif numplot == 8:
     114        sub = 331 #241
     115        fig.subplots_adjust(wspace = 0.1, hspace = 0.3)
     116    else:
     117        print "supported: 1,2,3,4,6,8"
     118        exit()
     119
     120    #####################
     121    ### Prepare time loop
    20122    nt = len(u[:,0,0,0])
    21     if var not in nc.variables:
    22         print "not found in file:",var
    23         exit()
    24     else:   
    25         dimension = np.array(nc.variables[var]).ndim
    26         if dimension == 3:     what_I_plot = nc.variables[var][:,:,:]
    27         elif dimension == 4:   what_I_plot = nc.variables[var][:,nvert,:,:] 
    28     fig = figure()
    29     fig.subplots_adjust(wspace = 0.0, hspace = 0.3)
    30     sub = 221
    31     for i in range(0,nt-1,int(nt/4.)):
    32        subplot(sub)
    33        ptitle("winds+"+var+" t"+str(i)+" l"+str(nvert))
     123    if nt <= numplot or numplot == 1: 
     124        print "I am plotting only one map ",nt,numplot
     125        tabrange = [0]
     126    else:                         
     127        tabrange = range(0,nt-1,int(nt/numplot))
     128
     129    #################################
     130    ### Time loop for plotting device
     131    for i in tabrange:
     132
     133       print i
     134
     135       ### General plot settings
     136       if tabrange != [0]: subplot(sub)
     137       zetitle = "WINDS" + "_"
     138
     139       ### Map projection
    34140       m = define_proj(proj,wlon,wlat,back=back)
    35141       x, y = m(lon2d, lat2d)
    36        contourf(x, y, what_I_plot[i,:,:], 30)
    37        vectorfield(u[i,nvert,:,:], v[i,nvert,:,:],\
     142
     143       #### Contour plot
     144       if var:
     145           zetitle = zetitle + var + "_"
     146           if typefile in ['mesoapi','meso']:    what_I_plot = dumpbdy(field[i,:,:])
     147           elif typefile in ['gcm']:             
     148               if dimension == 2:                what_I_plot = field[:,:]
     149               elif dimension == 3:              what_I_plot = field[i,:,:]
     150           contourf(x, y, what_I_plot, 30)
     151
     152       ### Vector plot
     153       if   typefile in ['mesoapi','meso']:   
     154           [vecx,vecy] = [dumpbdy(u[i,nvert,:,:]), dumpbdy(v[i,nvert,:,:])]
     155           key = True
     156       elif typefile in ['gcm']:               
     157           [vecx,vecy] = [        u[i,nvert,:,:] ,         v[i,nvert,:,:] ]
     158           key = False
     159       if metwind:  [vecx,vecy] = m.rotate_vector(vecx, vecy, lon2d, lat2d)
     160       else:        print "Beware ! Not using meteorological winds. You trust numerical grid as being (x,y)."
     161       vectorfield(vecx, vecy,\
    38162                      x, y, stride=stride, csmooth=2,\
    39                       scale=20., factor=300., color='k')
     163                      scale=15., factor=200., color='k', key=key)
     164       
     165       ### Next subplot
     166       zetitle = zetitle + "LT"+str((ltst+i)%24)
     167       ptitle(zetitle)
    40168       sub += 1
    41     if not target:   zeplot = namefile+".winds"+var+str(nvert)
    42     else:            zeplot = target+"/winds"+var+str(nvert)
     169
     170    ##########################################################################
     171    ### Save the figure in a file in the data folder or an user-defined folder
     172    if not target:   zeplot = namefile+zetitle
     173    else:            zeplot = target+"/"+zetitle
    43174    makeplotpng(zeplot,pad_inches_value=0.35)   
    44175
    45 def getwinds (nc,charu='U',charv='V'):
     176
     177
     178####################################################
     179####################################################
     180### A simple program to get wind vectors' components
     181def getwinds (nc,charu='Um',charv='Vm'):
    46182    import numpy as np
    47183    u = nc.variables[charu]
    … …  
    49185    if charu == 'U': u = u[:, :, :, 0:len(u[0,0,0,:])-1]
    50186    if charv == 'V': v = v[:, :, 0:len(v[0,0,:,0])-1, :]
    51     #print np.array(u).shape, np.array(v).shape
     187                     ### ou alors prendre les coordonnees speciales
    52188    return u,v
    53189
     190
     191
     192###########################################################################################
     193###########################################################################################
     194### What is below relate to running the file as a command line executable (very convenient)
    54195if __name__ == "__main__":
    55196    import sys
    56     ### to be replaced by argparse
    57     from optparse import OptionParser
     197    from optparse import OptionParser    ### to be replaced by argparse
    58198    parser = OptionParser()
    59199    parser.add_option('-f', action='store', dest='namefile',    type="string",  default=None,  help='name of WRF file [NEEDED]')
    60200    parser.add_option('-l', action='store', dest='nvert',       type="int",     default=0,     help='subscript for vertical level')
    61     parser.add_option('-p', action='store', dest='proj',        type="string",  default='cyl', help='projection')
     201    parser.add_option('-p', action='store', dest='proj',        type="string",  default=None, help='projection')
    62202    parser.add_option('-b', action='store', dest='back',        type="string",  default=None,  help='background')
    63203    parser.add_option('-t', action='store', dest='target',      type="string",  default=None,  help='destination folder')
    64     parser.add_option('-s', action='store', dest='stride',      type="int",     default=5,     help='stride vectors')
    65     parser.add_option('-v', action='store', dest='var',         type="string",  default='HGT', help='variable contoured')
     204    parser.add_option('-s', action='store', dest='stride',      type="int",     default=3,     help='stride vectors')
     205    parser.add_option('-v', action='store', dest='var',         type="string",  default=None,  help='variable contoured')
     206    parser.add_option('-n', action='store', dest='numplot',     type="int",     default=4,     help='number of plots')
    66207    (opt,args) = parser.parse_args()
    67208    if opt.namefile is None:
    … …  
    69210        exit()
    70211    print "Options:", opt
    71     winds (opt.namefile,opt.nvert,proj=opt.proj,back=opt.back,target=opt.target,stride=opt.stride,var=opt.var)
     212    winds (opt.namefile,opt.nvert,proj=opt.proj,back=opt.back,target=opt.target,stride=opt.stride,var=opt.var,numplot=opt.numplot)
     213
     214#    if typefile in ['gcm']:
     215#        if var == 'HGT':    var = 'phisinit'  ## default choice for GCM
     216
     217
Note: See TracChangeset for help on using the changeset viewer.