Ignore:
Timestamp:
Jul 23, 2011, 1:38:02 AM (13 years ago)
Author:
aslmd
Message:

MESOSCALE: python graphics: contour plots + various small but convenient improvements. fixed broken domain.py

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

Legend:

Unmodified
Added
Removed

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

    r240 r241  
    7171    return output, error
    7272
    73 def latinterv (area):
    74         if   area == "Europe":
    75                 wlat = [20.,80.]
    76                 wlon = [-50.,50.]
    77         elif area == "Central_America":
    78                 wlat = [-10.,40.]
    79                 wlon = [230.,300.]
    80         elif area == "Africa":
    81                 wlat = [-20.,50.]
    82                 wlon = [-50.,50.]
    83         elif area == "Whole":
    84                 wlat = [-90.,90.]
    85                 wlon = [-180.,180.]
    86         elif area == "Southern_Hemisphere":
    87                 wlat = [-90.,60.]
    88                 wlon = [-180.,180.]
    89         elif area == "Northern_Hemisphere":
    90                 wlat = [-60.,90.]
    91                 wlon = [-180.,180.]
    92         elif area == "Tharsis":
    93                 wlat = [-30.,60.]
    94                 wlon = [-170.,-10.]
    95         elif area == "Whole_No_High":
    96                 wlat = [-60.,60.]
    97                 wlon = [-180.,180.]
    98         elif area == "Chryse":
    99                 wlat = [-60.,60.]
    100                 wlon = [-60.,60.]
    101         elif area == "North_Pole":
    102                 wlat = [50.,90.]
    103                 wlon = [-180.,180.]
    104         elif area == "Close_North_Pole":
    105                 wlat = [75.,90.]
    106                 wlon = [-180.,180.]
    107         elif area == "South_Pole":
    108                 wlat = [-90.,-50.]
    109                 wlon = [-180.,180.]
    110         elif area == "Close_South_Pole":
    111                 wlat = [-90.,-75.]
    112                 wlon = [-180.,180.]
    113         return wlon,wlat
     73#def latinterv (area):
     74#       if   area == "Europe":                [wlat,wlon] = [[ 20., 80.],[- 50.,  50.]]
     75#       elif area == "Central_America":       [wlat,wlon] = [[-10., 40.],[ 230., 300.]]
     76#       elif area == "Africa":                [wlat,wlon] = [[-20., 50.],[- 50.,  50.]]
     77#       elif area == "Whole":                 [wlat,wlon] = [[-90., 90.],[-180.,-180.]]
     78#       elif area == "Southern_Hemisphere":   [wlat,wlon] = [[-90., 60.],[-180.,-180.]]
     79#        elif area == "Northern_Hemisphere":   [wlat,wlon] = [[-60., 90.],[-180.,-180.]]
     80#        elif area == "Tharsis":               [wlat,wlon] = [[-30., 60.],[-170.,- 10.]]
     81#       elif area == "Whole_No_High":         [wlat,wlon] = [[-60., 60.],[-180., 180.]]
     82#       elif area == "Chryse":                [wlat,wlon] = [[-60., 60.],[- 60.,  60.]]
     83#        elif area == "North_Pole":            [wlat,wlon] = [[ 50., 90.],[-180., 180.]]
     84#        elif area == "Close_North_Pole":      [wlat,wlon] = [[ 75., 90.],[-180., 180.]]
     85#        elif area == "Far_South_Pole":        [wlat,wlon] = [[-90.,-40.],[-180., 180.]]
     86#        elif area == "South_Pole":            [wlat,wlon] = [[-90.,-50.],[-180., 180.]]
     87#        elif area == "Close_South_Pole":      [wlat,wlon] = [[-90.,-75.],[-180., 180.]]
     88#       return wlon,wlat
    11489
    11590def definesubplot ( numplot, fig ):
    … …  
    142117    elif numplot == 1:
    143118        sub = 99999
    144     elif numplot < 0:
     119    elif numplot <= 0:
    145120        sub = 99999
    146121    else:
    … …  
    182157        zeday = int(zetime[8]+zetime[9]) - int(zetimestart[8]+zetimestart[9])
    183158        if zeday < 0:    lschar=""  ## might have crossed a month... fix soon
    184         else:            lschar="_Ls"+str( int( sol2ls ( getattr( nc, 'JULDAY' ) + zeday ) ) )
     159        else:            lschar="_Ls"+str( int( 10. * sol2ls ( getattr( nc, 'JULDAY' ) + zeday ) ) / 10. )
    185160        ###
    186161        zetime2 = nc.variables['Times'][1]
    … …  
    472447    return what_I_plot
    473448
     449def nolow(what_I_plot):
     450    from mymath import max,min
     451    lim = 0.15*0.5*(abs(max(what_I_plot))+abs(min(what_I_plot)))
     452    print "on vire en dessous de ", lim
     453    what_I_plot [ abs(what_I_plot) < lim ] = 1.e40
     454    return what_I_plot
     455
    474456def zoomset (wlon,wlat,zoom):
    475457    dlon = abs(wlon[1]-wlon[0])/2.
    … …  
    492474             "TAU_ICE":      "%.2f",\
    493475             "anomaly":      "%.1f",\
     476             "W":            "%.1f",\
    494477                    }
    495478    if whichvar not in fmtvar:
    … …  
    509492             "ICETOT":       "YlGnBu",\
    510493             "TAU_ICE":      "Blues",\
     494             "W":            "jet",\
    511495             "anomaly":      "RdBu_r",\
    512496                     }
     497#W --> spectral ou jet
    513498#spectral BrBG RdBu_r
     499    print "predefined colorbars"
    514500    if whichone not in whichcolorb:
    515501        whichone = "def"
    … …  
    567553        return whichlink
    568554
     555def latinterv (area="Whole"):
     556    list =    { \
     557        "Europe":                [[ 20., 80.],[- 50.,  50.]],\
     558        "Central_America":       [[-10., 40.],[ 230., 300.]],\
     559        "Africa":                [[-20., 50.],[- 50.,  50.]],\
     560        "Whole":                 [[-90., 90.],[-180.,-180.]],\
     561        "Southern_Hemisphere":   [[-90., 60.],[-180.,-180.]],\
     562        "Northern_Hemisphere":   [[-60., 90.],[-180.,-180.]],\
     563        "Tharsis":               [[-30., 60.],[-170.,- 10.]],\
     564        "Whole_No_High":         [[-60., 60.],[-180., 180.]],\
     565        "Chryse":                [[-60., 60.],[- 60.,  60.]],\
     566        "North_Pole":            [[ 50., 90.],[-180., 180.]],\
     567        "Close_North_Pole":      [[ 75., 90.],[-180., 180.]],\
     568        "Far_South_Pole":        [[-90.,-40.],[-180., 180.]],\
     569        "South_Pole":            [[-90.,-50.],[-180., 180.]],\
     570        "Close_South_Pole":      [[-90.,-75.],[-180., 180.]],\
     571              }
     572    if area not in list:   area = "Whole"
     573    [olat,olon] = list[area]
     574    return olon,olat
     575

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

    r225 r241  
    55def domain (namefile,proj=None,back="vishires",target=None):
    66    from netCDF4 import Dataset
    7     from myplot import getcoord2d,define_proj,makeplotpng,simplinterv,getprefix,dumpbdy,getproj,latinterv,wrfinterv,simplinterv
    8     from matplotlib.pyplot import contourf,rcParams
     7    from myplot import getcoord2d,define_proj,makeplotres,simplinterv,getprefix,dumpbdy,getproj,latinterv,wrfinterv,simplinterv
     8    from mymath import max,min
     9    from matplotlib.pyplot import contourf,rcParams,pcolor
    910    from numpy.core.defchararray import find
     11    from numpy import arange
    1012    ###
    1113    nc  = Dataset(namefile)
    … …  
    2325        zeplot = getprefix(nc) + "domain"
    2426    ###
    25     lon2d = dumpbdy(lon2d)
    26     lat2d = dumpbdy(lat2d)
     27    lon2d = dumpbdy(lon2d,5)
     28    lat2d = dumpbdy(lat2d,5)
    2729    if proj == "npstere":             [wlon,wlat] = latinterv("North_Pole")
    2830    elif proj in ["lcc","laea"]:      [wlon,wlat] = wrfinterv(lon2d,lat2d)
    … …  
    3234    x, y = m(lon2d, lat2d)
    3335    ###
    34     what_I_plot = dumpbdy( nc.variables[var][0,:,:] )
    35     contourf(x, y, what_I_plot, 40)
     36    what_I_plot = dumpbdy(nc.variables[var][0,:,:], 5)
     37    #levinterv = 250.
     38    #zelevels = arange(min(what_I_plot)-levinterv,max(what_I_plot)+levinterv,levinterv)
     39    zelevels = 30
     40    contourf(x, y, what_I_plot, zelevels)
     41    #pcolor(x,y,what_I_plot)  ## on voit trop les lignes !
    3642    ###
    3743    if not target:   zeplot = namefile[0:find(namefile,'wrfout')] + zeplot
    3844    else:            zeplot = target + "/" + zeplot         
    3945    ###
    40     makeplotpng(zeplot,pad_inches_value=0.35)
     46    pad_inches_value = 0.35
     47    makeplotres(zeplot,res=100.,pad_inches_value=pad_inches_value) #,erase=True)  ## a miniature
     48    makeplotres(zeplot,res=200.,pad_inches_value=pad_inches_value,disp=False)
     49    #makeplotpng(zeplot,pad_inches_value=0.35)
    4150    #rcParams['savefig.facecolor'] = 'black'
    4251    #makeplotpng(zeplot+"b",pad_inches_value=0.35)

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

    r240 r241  
    1616           numplot=2,\
    1717           var=None,\
    18            colorb=True,\
     18           colorb="def",\
    1919           winds=True,\
    2020           addchar=None,\
    … …  
    2828           hole=False,\
    2929           save="gui",\
    30            anomaly=False):
     30           anomaly=False,\
     31           var2=None):
    3132
    3233    ####################################################################################################################
    … …  
    3839    from myplot import getcoord2d,define_proj,makeplotres,simplinterv,vectorfield,ptitle,latinterv,getproj,wrfinterv,dumpbdy,\
    3940                       fmtvar,definecolorvec,defcolorb,getprefix,putpoints,calculate_bounds,errormess,definesubplot,\
    40                        zoomset,getcoorddef,getwinddef,whatkindfile,reducefield,bounds,getstralt,getfield,smooth
     41                       zoomset,getcoorddef,getwinddef,whatkindfile,reducefield,bounds,getstralt,getfield,smooth,nolow
    4142    from mymath import deg,max,min,mean
    4243    from matplotlib.pyplot import contour,contourf, subplot, figure, rcParams, savefig, colorbar, pcolor, show
    … …  
    4546    from numpy.core.defchararray import find
    4647
    47     #rcParams['backend'] = 'PS'
     48    #if save == 'eps': rcParams['backend'] = 'PS'
    4849
    4950    ######################
    … …  
    6465    ### Define plot boundaries
    6566    if proj == "npstere":             [wlon,wlat] = latinterv("North_Pole")
    66     elif proj == "spstere":           [wlon,wlat] = latinterv("South_Pole")
     67    elif proj == "spstere":           [wlon,wlat] = latinterv("Far_South_Pole")
    6768    elif proj in ["lcc","laea"]:      [wlon,wlat] = wrfinterv(lon2d,lat2d)
    6869    else:                             [wlon,wlat] = simplinterv(lon2d,lat2d)
    … …  
    7778        print nc.variables                 
    7879        errormess("please set at least winds or var")
     80    if anomaly:           basename = 'd' + basename
    7981    basename = basename + getstralt(nc,nvert)  ## can be moved elsewhere for a more generic routine
    8082
    … …  
    9092    nplot = 1
    9193    error = False
    92     if itstep is None: itstep = int(24./numplot)
     94    if itstep is None and numplot > 0: itstep = int(24./numplot)
     95    else:                              itstep = 1
    9396    while error is False:
    9497
    … …  
    109112       ### If only one local time is requested (numplot < 0)
    110113       elif numplot <= 0:
    111            if int(ltst) + numplot != 0:        
     114           if int(ltst) + numplot != 0:
    112115                 itime += 1
    113116                 if found_lct is True: break     ## because it means LT was found at previous iteration
    … …  
    121124
    122125       #### Contour plot
     126       if var2:
     127           what_I_contour, error = reducefield( getfield(nc,var2), d4=itime, d3=nvert )
     128           if not error:
     129               if typefile in ['mesoapi','meso']:    what_I_contour = dumpbdy(what_I_contour,6)
     130               zevmin, zevmax = calculate_bounds(what_I_contour)
     131               zelevels = np.linspace(zevmin,zevmax,num=20)
     132               if var2 == 'HGT':  zelevels = np.arange(-10000.,30000.,2000.)
     133               contour( x, y, what_I_contour, zelevels, colors='k', linewidths = 0.33 ) #colors='w' )# , alpha=0.5)
     134
     135       #### Shaded plot
    123136       if var:
    124137           what_I_plot, error = reducefield( getfield(nc,var), d4=itime, d3=nvert )
    … …  
    132145               if typefile in ['mesoapi','meso']:    what_I_plot = dumpbdy(what_I_plot,6)
    133146               zevmin, zevmax = calculate_bounds(what_I_plot,vmin=vmin,vmax=vmax)
    134                if colorb is True: colorb = defcolorb(fvar)
    135                palette = get_cmap(name=colorb)
     147               if colorb in ["def","nobar"]:   palette = get_cmap(name=defcolorb(fvar))
     148               else:                           palette = get_cmap(name=colorb)
    136149               if not tile:
    137150                   if not hole: what_I_plot = bounds(what_I_plot,zevmin,zevmax)
    138                    zelevels = np.linspace(zevmin,zevmax)
     151                   zelevels = np.linspace(zevmin,zevmax) #,num=20)
    139152                   contourf( x, y, what_I_plot, zelevels, cmap = palette )
    140                else:   
    141                    pcolor( x, y, what_I_plot, cmap = palette, vmin=zevmin, vmax=zevmax )
    142                #putpoints(m,fig)
    143                if var in ['HGT']: pass
    144                elif colorb:             
     153               else:
     154                   if hole:  what_I_plot = nolow(what_I_plot)
     155                   pcolor( x, y, what_I_plot, cmap = palette, \
     156                           vmin=zevmin, vmax=zevmax )
     157               if colorb != 'nobar' or var == 'HGT':             
    145158                         ndiv = 10
    146159                         colorbar(fraction=0.05,pad=0.1,format=fmtvar(fvar),\
    … …  
    148161                                           extend='neither',spacing='proportional')
    149162                                           # both min max neither
    150                   
     163           
    151164       ### Vector plot
    152165       if winds:
    … …  
    192205        pad_inches_value = 0.35
    193206        if save == 'png':
    194             makeplotres(zeplot,res=100.,pad_inches_value=pad_inches_value,erase=True)  ## a miniature
     207            makeplotres(zeplot,res=100.,pad_inches_value=pad_inches_value) #,erase=True)  ## a miniature
    195208            makeplotres(zeplot,res=200.,pad_inches_value=pad_inches_value,disp=False)
    196209        elif save in ['eps','svg','pdf']:
    … …  
    214227        outtab = ones(zelen) * -999999
    215228    else:
    216         print zelen, len(tab)
    217229        if zelen != len(tab):
    218230            print "not enough or too much values... setting same values all variables"
    … …  
    242254    parser.add_option('-t', action='store',dest='target',       type="string",  default=None,  help='destination folder')
    243255    parser.add_option('-s', action='store',dest='stride',       type="int",     default=3,     help='stride vectors (def=3)')
    244     parser.add_option('-v', action='append',dest='var',         type="string",  default=None,  help='variable contoured (append)')
     256    parser.add_option('-v', action='append',dest='var',         type="string",  default=None,  help='variable color-shaded (append)')
    245257    parser.add_option('-n', action='store',dest='numplot',      type="int",     default=2,     help='number of plots (def=1)(<0: 1 plot of LT -*numplot*)')
    246258    parser.add_option('-i', action='store',dest='interp',       type="int",     default=None,  help='interpolation method (2: press, 3: z-amr, 4:z-als)')
    247     parser.add_option('-c', action='store',dest='colorb',       type="string",  default=True,  help='change colormap')
     259    parser.add_option('-c', action='store',dest='colorb',       type="string",  default="def", help='change colormap (nobar: no colorbar)')
    248260    parser.add_option('-x', action='store_false',dest='winds',                  default=True,  help='no wind vectors')
    249261    parser.add_option('-m', action='append',dest='vmin',        type="float",   default=None,  help='bounding minimum value (append)')   
    … …  
    257269    parser.add_option('-S', action='store',dest='save',         type="string",  default="gui", help='save mode (png,eps,svg,pdf or gui)(def=gui)')
    258270    parser.add_option('-a', action='store_true',dest='anomaly',                 default=False, help='compute and plot relative anomaly in %')
     271    parser.add_option('-w', action='store',dest='var2',         type="string",  default=None,  help='variable contoured')
    259272    #parser.add_option('-V', action='store', dest='comb',        type="float",   default=None,  help='a defined combination of variables')
    260273    (opt,args) = parser.parse_args()
    … …  
    262275        print "I want to eat one file at least ! Use winds.py -f name_of_my_file. Or type winds.py -h"
    263276        exit()
     277    if opt.var is None and opt.anomaly is True:
     278        print "Cannot ask to compute anomaly if no variable is set"
     279        exit()   
    264280    print "Options:", opt
    265281
    … …  
    296312            if zefields == ''       :  zefields = listvar
    297313            else                    :  zefields = zefields + "," + listvar
     314            if opt.var2 is not None : zefields = zefields + "," + opt.var2 
    298315            print zefields
    299316            zefile = api_onelevel (  path_to_input   = '', \
    … …  
    327344                tile=opt.tile,zoom=opt.zoom,display=opt.display,\
    328345                itstep=opt.itstep,hole=opt.hole,save=opt.save,\
    329                 anomaly=opt.anomaly)
     346                anomaly=opt.anomaly,var2=opt.var2)
    330347            print 'Done: '+name
    331348            system("rm -f to_be_erased")
Note: See TracChangeset for help on using the changeset viewer.