Changeset 184


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
Files:
3 edited

Legend:

Unmodified
Added
Removed

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

    r181 r184  
    3535        return wlon,wlat
    3636
     37def getproj (nc):
     38    map_proj = getattr(nc, 'MAP_PROJ')
     39    cen_lat  = getattr(nc, 'CEN_LAT')
     40    if map_proj == 2:
     41        if cen_lat > 10.:   
     42            proj="npstere"
     43            print "NP stereographic polar domain"
     44        else:           
     45            proj="spstere"
     46            print "SP stereographic polar domain"
     47    elif map_proj == 1:
     48        print "lambert projection domain"
     49        proj="lcc"
     50    elif map_proj == 3:
     51        print "mercator projection"
     52        proj="merc"
     53    return proj   
     54
    3755def ptitle (name):
    3856    from matplotlib.pyplot import title
    … …  
    4361    import numpy as np
    4462    return [[np.min(lon2d),np.max(lon2d)],[np.min(lat2d),np.max(lat2d)]]
     63
     64def wrfinterv (lon2d,lat2d):
     65    nx = len(lon2d[0,:])-1
     66    ny = len(lon2d[:,0])-1
     67    return [[lon2d[0,0],lon2d[nx,ny]],[lat2d[0,0],lat2d[nx,ny]]]
    4568
    4669def makeplotpngres (filename,res,pad_inches_value=0.25,folder='',disp=True):
    … …  
    5881    return
    5982
    60 def getcoord2d (nc,nlat='XLAT',nlon='XLONG'):
    61         import numpy as np
    62         lat = nc.variables[nlat][0,:,:]
    63         lon = nc.variables[nlon][0,:,:]
    64         if np.array(lat).ndim != 2:     [lon2d,lat2d] = np.meshgrid(lon,lat)
    65         else:                           [lon2d,lat2d] = [lon,lat]
    66         return lon2d,lat2d
     83def dumpbdy (field):
     84    nx = len(field[0,:])-1
     85    ny = len(field[:,0])-1
     86    return field[5:ny-5,5:nx-5]
     87
     88def getcoord2d (nc,nlat='XLAT',nlon='XLONG',is1d=False):
     89    import numpy as np
     90    if is1d:
     91        lat = nc.variables[nlat][:]
     92        lon = nc.variables[nlon][:]
     93        [lon2d,lat2d] = np.meshgrid(lon,lat)
     94    else:
     95        lat = nc.variables[nlat][0,:,:]
     96        lon = nc.variables[nlon][0,:,:]
     97        [lon2d,lat2d] = [lon,lat]
     98    return lon2d,lat2d
    6799
    68100def smooth (field, coeff):
    … …  
    94126        return improc
    95127
    96 def vectorfield (u, v, x, y, stride=3, scale=15., factor=250., color='black', csmooth=1):
    97         ## scale regle la reference du vecteur
    98         ## factor regle toutes les longueurs (dont la reference). l'AUGMENTER pour raccourcir les vecteurs.
    99         import  matplotlib.pyplot               as plt
    100         import  numpy                           as np
    101         posx = np.max(x)*0.90
    102         posy = np.mean(y)
    103         u = smooth(u,csmooth)
    104         v = smooth(v,csmooth)
    105         q = plt.quiver( x[::stride,::stride],\
    106                         y[::stride,::stride],\
    107                         u[::stride,::stride],\
    108                         v[::stride,::stride],\
    109                         angles='xy',color=color,\
    110                         scale=factor,width=0.003 )
    111         if color=='white':      kcolor='black'
    112         elif color=='yellow':   kcolor=color
    113         else:                   kcolor=color
    114         p = plt.quiverkey(q,posx,posy,scale,\
    115                         str(int(scale)),coordinates='data',color=kcolor)
    116         return p
     128def vectorfield (u, v, x, y, stride=3, scale=15., factor=250., color='black', csmooth=1, key=True):
     129    ## scale regle la reference du vecteur
     130    ## factor regle toutes les longueurs (dont la reference). l'AUGMENTER pour raccourcir les vecteurs.
     131    import  matplotlib.pyplot               as plt
     132    import  numpy                           as np
     133    posx = np.max(x) + np.std(x) / 3.  ## pb pour les domaines globaux ...
     134    posy = np.mean(y)
     135    #posx = np.min(x)
     136    #posy = np.max(x)
     137    u = smooth(u,csmooth)
     138    v = smooth(v,csmooth)
     139    widthvec = 0.005 #0.003
     140    q = plt.quiver( x[::stride,::stride],\
     141                    y[::stride,::stride],\
     142                    u[::stride,::stride],\
     143                    v[::stride,::stride],\
     144                    angles='xy',color=color,\
     145                    scale=factor,width=widthvec )
     146    if color=='white':     kcolor='black'
     147    elif color=='yellow':  kcolor=color
     148    else:                  kcolor=color
     149    if key: p = plt.quiverkey(q,posx,posy,scale,\
     150                   str(int(scale)),coordinates='data',color=kcolor)
     151    return
    117152
    118153def display (name):
    119         from    os                              import system
    120         system("display "+name+" > /dev/null 2> /dev/null &")
    121         return name
     154    from os import system
     155    system("display "+name+" > /dev/null 2> /dev/null &")
     156    return name
    122157
    123158def findstep (wlon):
    124     steplon = int((wlon[1]-wlon[0])/3.)
    125     step = 60.
    126     if steplon < 60.: step = 30.
    127     if steplon < 30.: step = 15.
    128     if steplon < 15.: step = 10.
    129     if steplon < 10.: step = 5.
    130     if steplon < 5.:  step = 1.
     159    steplon = int((wlon[1]-wlon[0])/4.)  #3
     160    step = 120.
     161    while step > steplon and step > 15. :       step = step / 2.
     162    if step <= 15.:
     163        while step > steplon and step > 5.  :   step = step - 5.
     164    if step <= 5.:
     165        while step > steplon and step > 1.  :   step = step - 1.
     166    if step <= 1.:
     167        step = 1.
    131168    return step
    132169
    … …  
    137174    meanlon = 0.5*(wlon[0]+wlon[1])
    138175    meanlat = 0.5*(wlat[0]+wlat[1])
     176    if   wlat[0] >= 80.:   blat =  40.
     177    elif wlat[0] <= -80.:  blat = -40.
     178    else:                  blat = wlat[0]
    139179    h = 2000.
    140     if   char == "cyl":     m = Basemap(projection='cyl',llcrnrlat=wlat[0],urcrnrlat=wlat[1],llcrnrlon=wlon[0],urcrnrlon=wlon[1])
    141     elif char == "moll":    m = Basemap(projection='moll',lon_0=meanlon)
    142     elif char == "ortho":   m = Basemap(projection='ortho',lon_0=meanlon,lat_0=meanlat)
    143     elif char == "lcc":     m = Basemap(projection='lcc',lat_1=meanlat,lat_0=meanlat,lon_0=meanlon,\
     180    radius = 3397200
     181    if   char == "cyl":     m = Basemap(rsphere=radius,projection='cyl',\
    144182                              llcrnrlat=wlat[0],urcrnrlat=wlat[1],llcrnrlon=wlon[0],urcrnrlon=wlon[1])
    145     elif char == "npstere": m = Basemap(projection='npstere', boundinglat=wlat[0], lon_0=0.)
    146     elif char == "spstere": m = Basemap(projection='spstere', boundinglat=wlat[0], lon_0=0.)
    147     elif char == "nsper":   m = Basemap(projection='nsper',lon_0=meanlon,lat_0=meanlat,satellite_height=h*1000.)
    148     fontsizemer = int(mpl.rcParams['font.size']*2./3.)
    149     if char in ["cyl","lcc"]:   step = findstep(wlon)
    150     else:                       step = 10.
     183    elif char == "moll":    m = Basemap(rsphere=radius,projection='moll',lon_0=meanlon)
     184    elif char == "ortho":   m = Basemap(rsphere=radius,projection='ortho',lon_0=meanlon,lat_0=meanlat)
     185    elif char == "lcc":     m = Basemap(rsphere=radius,projection='lcc',lat_1=meanlat,lat_0=meanlat,lon_0=meanlon,\
     186                              llcrnrlat=wlat[0],urcrnrlat=wlat[1],llcrnrlon=wlon[0],urcrnrlon=wlon[1])
     187    elif char == "npstere": m = Basemap(rsphere=radius,projection='npstere', boundinglat=blat, lon_0=0.)
     188    elif char == "spstere": m = Basemap(rsphere=radius,projection='spstere', boundinglat=blat, lon_0=0.)
     189    elif char == "nsper":   m = Basemap(rsphere=radius,projection='nsper',lon_0=meanlon,lat_0=meanlat,satellite_height=h*1000.)
     190    elif char == "merc":    m = Basemap(rsphere=radius,projection='merc',lat_ts=0.,\
     191                              llcrnrlat=wlat[0],urcrnrlat=wlat[1],llcrnrlon=wlon[0],urcrnrlon=wlon[1])
     192    fontsizemer = int(mpl.rcParams['font.size']*3./4.)
     193    if char in ["cyl","lcc","merc"]:   step = findstep(wlon)
     194    else:                              step = 10.
    151195    m.drawmeridians(np.r_[-180.:180.:step*2.], labels=[0,0,0,1], color='grey', fontsize=fontsizemer)
    152196    m.drawparallels(np.r_[-90.:90.:step], labels=[1,0,0,0], color='grey', fontsize=fontsizemer)

  • 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.