Changeset 240


Ignore:
Timestamp:
Jul 22, 2011, 9:46:21 AM (13 years ago)
Author:
aslmd
Message:

MESOSCALE: corrected a bug in the previous commit in makemeso, works now fine. added anomaly calculations in python interface + better treatment for output; default is now GUI where one can zoom etc... and further explore the plotted field

Location:
trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/MESOSCALE/LMD_MM_MARS/makemeso

    r239 r240  
    538538      *) echo not yet supported ; exit ;;
    539539  esac
    540   restex=$(expr ${physx}%${divx})   
    541   restey=$(expr ${physy}%${divy})
     540  physx=$(expr ${lon} - 1)
     541  restex=$(expr ${physx} \% ${divx})   
     542  physy=$(expr ${lat} - 1)
     543  restey=$(expr ${physy} \% ${divy})
    542544  if [[ ${restex} != 0 || ${restey} != 0 ]]
    543545  then
    … …  
    546548     exit
    547549  fi
    548   physx=$(expr ${lon} - 1) ; physx=$(expr ${physx} \/ ${divx})
    549   physy=$(expr ${lat} - 1) ; physy=$(expr ${physy} \/ ${divy})
     550  physx=$(expr ${physx} \/ ${divx})
     551  physy=$(expr ${physy} \/ ${divy})
    550552  physz=$(expr ${level} - 1)
    551553

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

    r238 r240  
    100100                wlon = [-60.,60.]
    101101        elif area == "North_Pole":
    102                 wlat = [60.,90.]
     102                wlat = [50.,90.]
    103103                wlon = [-180.,180.]
    104104        elif area == "Close_North_Pole":
    105105                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.]
    106112                wlon = [-180.,180.]
    107113        return wlon,wlat
    … …  
    244250    return [wlon,wlat]
    245251
    246 def makeplotpngres (filename,res,pad_inches_value=0.25,folder='',disp=True):
     252def makeplotres (filename,res=None,pad_inches_value=0.25,folder='',disp=True,ext='png',erase=False):
    247253    import  matplotlib.pyplot as plt
    248     res = int(res)
    249     name = filename+"_"+str(res)+".png"
     254    from os import system
     255    addstr = ""
     256    if res is not None:
     257        res = int(res)
     258        addstr = "_"+str(res)
     259    name = filename+addstr+"."+ext
    250260    if folder != '':      name = folder+'/'+name
    251261    plt.savefig(name,dpi=res,bbox_inches='tight',pad_inches=pad_inches_value)
    252     if disp:              display(name)         
     262    if disp:              display(name)
     263    if ext in ['eps','ps','svg']:   system("tar czvf "+name+".tar.gz "+name+" ; rm -f "+name)
     264    if erase:   system("mv "+name+" to_be_erased")             
    253265    return
    254266
    255 def makeplotpng (filename,pad_inches_value=0.25,minres=100.,folder='',disp=True):
    256     makeplotpngres(filename,minres,     pad_inches_value=pad_inches_value,folder=folder,disp=disp)
    257     makeplotpngres(filename,minres+200.,pad_inches_value=pad_inches_value,folder=folder,disp=False)
    258     return
    259 
    260 def dumpbdy (field,stag=None):
     267def dumpbdy (field,n,stag=None):
    261268    nx = len(field[0,:])-1
    262269    ny = len(field[:,0])-1
    263270    if stag == 'U': nx = nx-1
    264271    if stag == 'V': ny = ny-1
    265     return field[5:ny-5,5:nx-5]
     272    return field[n:ny-n,n:nx-n]
    266273
    267274def getcoorddef ( nc ):   
    … …  
    270277    if typefile in ['mesoapi','meso']:
    271278        [lon2d,lat2d] = getcoord2d(nc)
    272         lon2d = dumpbdy(lon2d)
    273         lat2d = dumpbdy(lat2d)
     279        lon2d = dumpbdy(lon2d,6)
     280        lat2d = dumpbdy(lat2d,6)
    274281    elif typefile in ['gcm']:
    275282        [lon2d,lat2d] = getcoord2d(nc,nlat="latitude",nlon="longitude",is1d=True)
    … …  
    376383    import  numpy                           as np
    377384    import  matplotlib                      as mpl
     385    from mymath import max
    378386    meanlon = 0.5*(wlon[0]+wlon[1])
    379387    meanlat = 0.5*(wlat[0]+wlat[1])
    380388    if   wlat[0] >= 80.:   blat =  40.
    381     elif wlat[0] <= -80.:  blat = -40.
    382     else:                  blat = wlat[0]
     389    elif wlat[1] <= -80.:  blat = -40.
     390    else:                  blat = max([wlat[0],wlat[1]])
     391    print blat
    383392    h = 50.  ## en km
    384393    radius = 3397200.
    … …  
    482491             "ICETOT":       "%.1e",\
    483492             "TAU_ICE":      "%.2f",\
     493             "anomaly":      "%.1f",\
    484494                    }
    485495    if whichvar not in fmtvar:
    … …  
    499509             "ICETOT":       "YlGnBu",\
    500510             "TAU_ICE":      "Blues",\
     511             "anomaly":      "RdBu_r",\
    501512                     }
     513#spectral BrBG RdBu_r
    502514    if whichone not in whichcolorb:
    503515        whichone = "def"

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

    r238 r240  
    2626           display=True,\
    2727           itstep=None,\
    28            hole=False):
     28           hole=False,\
     29           save="gui",\
     30           anomaly=False):
    2931
    3032    ####################################################################################################################
    … …  
    3436    ### Load librairies and functions
    3537    from netCDF4 import Dataset
    36     from myplot import getcoord2d,define_proj,makeplotpng,makeplotpngres,simplinterv,vectorfield,ptitle,latinterv,getproj,wrfinterv,dumpbdy,\
     38    from myplot import getcoord2d,define_proj,makeplotres,simplinterv,vectorfield,ptitle,latinterv,getproj,wrfinterv,dumpbdy,\
    3739                       fmtvar,definecolorvec,defcolorb,getprefix,putpoints,calculate_bounds,errormess,definesubplot,\
    38                        zoomset,getcoorddef,getwinddef,whatkindfile,reducefield,bounds,getstralt,getfield
     40                       zoomset,getcoorddef,getwinddef,whatkindfile,reducefield,bounds,getstralt,getfield,smooth
    3941    from mymath import deg,max,min,mean
    40     from matplotlib.pyplot import contour,contourf, subplot, figure, rcParams, savefig, colorbar, pcolor
     42    from matplotlib.pyplot import contour,contourf, subplot, figure, rcParams, savefig, colorbar, pcolor, show
    4143    from matplotlib.cm import get_cmap
    4244    import numpy as np
    4345    from numpy.core.defchararray import find
    4446
     47    #rcParams['backend'] = 'PS'
     48
    4549    ######################
    4650    ### Load NETCDF object
    47     nc  = Dataset(namefile)
     51    nc  = Dataset(namefile) 
    4852   
    4953    ##################################
    … …  
    6064    ### Define plot boundaries
    6165    if proj == "npstere":             [wlon,wlat] = latinterv("North_Pole")
     66    elif proj == "spstere":           [wlon,wlat] = latinterv("South_Pole")
    6267    elif proj in ["lcc","laea"]:      [wlon,wlat] = wrfinterv(lon2d,lat2d)
    6368    else:                             [wlon,wlat] = simplinterv(lon2d,lat2d)
    … …  
    119124           what_I_plot, error = reducefield( getfield(nc,var), d4=itime, d3=nvert )
    120125           if not error:
    121                if typefile in ['mesoapi','meso']:    what_I_plot = dumpbdy(what_I_plot)
     126               fvar = var
     127               ###
     128               if anomaly:
     129                   what_I_plot = 100. * ((what_I_plot / smooth(what_I_plot,10)) - 1.)
     130                   fvar = 'anomaly'
     131               ###
     132               if typefile in ['mesoapi','meso']:    what_I_plot = dumpbdy(what_I_plot,6)
    122133               zevmin, zevmax = calculate_bounds(what_I_plot,vmin=vmin,vmax=vmax)
    123                if colorb is True: colorb = defcolorb(var)
     134               if colorb is True: colorb = defcolorb(fvar)
    124135               palette = get_cmap(name=colorb)
    125136               if not tile:
    … …  
    133144               elif colorb:             
    134145                         ndiv = 10
    135                          colorbar(fraction=0.05,pad=0.1,format=fmtvar(var),\
     146                         colorbar(fraction=0.05,pad=0.1,format=fmtvar(fvar),\
    136147                                           ticks=np.linspace(zevmin,zevmax,ndiv+1),\
    137                                            extend='max',spacing='proportional')
    138                                            # both min max neither
     148                                           extend='neither',spacing='proportional')
     149                                           # both min max neither 
    139150                 
    140151       ### Vector plot
    … …  
    144155           if not error:
    145156               if typefile in ['mesoapi','meso']:   
    146                    [vecx,vecy] = [dumpbdy(vecx,stag=uchar), dumpbdy(vecy,stag=vchar)]
     157                   [vecx,vecy] = [dumpbdy(vecx,6,stag=uchar), dumpbdy(vecy,6,stag=vchar)]
    147158                   key = True
    148159               elif typefile in ['gcm']:           
    … …  
    178189    else:               zeplot = target + "/" + zeplot 
    179190    ###
    180     if found_lct:     makeplotpng(zeplot,pad_inches_value=0.35,disp=display)   
    181     #if found_lct:     makeplotpngres(zeplot,200,disp=display)
     191    if found_lct:     
     192        pad_inches_value = 0.35
     193        if save == 'png':
     194            makeplotres(zeplot,res=100.,pad_inches_value=pad_inches_value,erase=True)  ## a miniature
     195            makeplotres(zeplot,res=200.,pad_inches_value=pad_inches_value,disp=False)
     196        elif save in ['eps','svg','pdf']:
     197            makeplotres(zeplot,         pad_inches_value=pad_inches_value,disp=False,ext=save)
     198        elif save == 'gui':
     199            show()
     200        else:
     201            print "save mode not supported. using gui instead."
     202            show()
    182203    else:             print "Local time not found"
    183204
    … …  
    210231    from netCDF4 import Dataset
    211232    from myplot import getlschar
     233    from os import system
    212234
    213235    #############################
    … …  
    233255    parser.add_option('-e', action='store',dest='itstep',       type="int",     default=None,  help='stride time (def=4)')
    234256    parser.add_option('-H', action='store_true',dest='hole',                    default=False, help='holes above max and below min')
     257    parser.add_option('-S', action='store',dest='save',         type="string",  default="gui", help='save mode (png,eps,svg,pdf or gui)(def=gui)')
     258    parser.add_option('-a', action='store_true',dest='anomaly',                 default=False, help='compute and plot relative anomaly in %')
    235259    #parser.add_option('-V', action='store', dest='comb',        type="float",   default=None,  help='a defined combination of variables')
    236260    (opt,args) = parser.parse_args()
    … …  
    275299            zefile = api_onelevel (  path_to_input   = '', \
    276300                                     input_name      = zefile, \
    277                                      #path_to_output  = opt.target, \
    278301                                     fields          = zefields, \
    279302                                     interp_method   = opt.interp, \
    … …  
    303326                addchar=lschar,interv=[zehour,zehourin],vmin=argvmin,vmax=argvmax,\
    304327                tile=opt.tile,zoom=opt.zoom,display=opt.display,\
    305                 itstep=opt.itstep,hole=opt.hole)
     328                itstep=opt.itstep,hole=opt.hole,save=opt.save,\
     329                anomaly=opt.anomaly)
    306330            print 'Done: '+name
    307    
     331            system("rm -f to_be_erased")
     332 
    308333        #########################################################
    309334        ### Generate a .sh file with the used command saved in it
Note: See TracChangeset for help on using the changeset viewer.