source: trunk/LMDZ.PLUTO/util/script_figures/mapch4cloud.py

#! /usr/bin/env python
from ppclass import pp
from    netCDF4               import    Dataset
from    numpy                 import    *
import  numpy                 as        np
import  matplotlib.pyplot     as        mpl
from matplotlib.cm import get_cmap
import pylab
from matplotlib import ticker
import matplotlib.colors as mcolors
import datetime
from mpl_toolkits.basemap import Basemap, shiftgrid

############################
filename1="diagfi2015.nc"
filename2="diagfi2015_S.nc"
var="tsurf" #variable
var1="ch4_ice_col"
var2="phisinit"
tint=["30.625","30.875","31.125","30.375"] #Time must be as written in the input file
tintstr=["03:00","09:00","15:00","21:00"] #Time must be as written in the input file

font=15

nc1=Dataset(filename1)
nc2=Dataset(filename2)

lat=nc1.variables["lat"][:]
lon=nc1.variables["lon"][:]

# altitdue file 2
alt=nc2.variables["altitude"][:]
############################

def getvar(filename,var,tint):
    myvar = pp(file=filename,var=var,t=tint,compute="nothing").getf()  # get data to be changed according to selected variable
    print((shape(myvar)))
    return myvar

def swinglon(myvar):
   # changer les longitudes pour mettre TR au centre
   vec=shape(myvar)
   myvarbis=np.zeros(vec,dtype='f')
# i lat : pas de changement
# j lon :
   for i in range(vec[0]):
       for j in range(vec[1]):
           if j < int(vec[1]/2.) :
             myvarbis[i,j]=myvar[i,j+int(vec[1]/2.)]
             # myvar2bis[i,j]=myvar2[i,j+int(vec[1]/2)]
           else:
             myvarbis[i,j]=myvar[i,j-int(vec[1]/2.)]
             # myvar2bis[i,j]=myvar2[i,j-int(vec[1]/2)]
   return myvarbis

def getwinds(lon,lat,vecx,vecy):
          svx='None'   # arrow every svx box
          svy='None'
          svx=1
          svy=1
          angle='uv'       # 'xy'
          color='black'    # arrow color
          pivot='mid'      # arrow around middle of box. Alternative : tip
          scale=2*33     # scale arrow
          width=0.002      # width arrow
          linewidths=0.5   # epaisseur contour arrow
          edgecolors='k'   # couleur contour arrow

    #  *scale*: [ *None* | float ]
    #  Data units per arrow length unit, e.g., m/s per plot width; a smaller
    #  scale parameter makes the arrow longer.  If *None*, a simple
    #  autoscaling algorithm is used, based on the average vector length
    #  and the number of vectors.  The arrow length unit is given by
    #  the *scale_units* parameter
   
    #  *scale_units*: *None*, or any of the *units* options.
    #  For example, if *scale_units* is 'inches', *scale* is 2.0, and
    #  ``(u,v) = (1,0)``, then the vector will be 0.5 inches long.
    #  If *scale_units* is 'width', then the vector will be half the width
    #  of the axes.
   
    #  If *scale_units* is 'x' then the vector will be 0.5 x-axis
    #  units.  To plot vectors in the x-y plane, with u and v having
    #  the same units as x and y, use
    #  "angles='xy', scale_units='xy', scale=1".
    
          x, y = np.meshgrid(lon,lat)
          q = mpl.quiver( x[::svy,::svx],y[::svy,::svx],\
           vecx[::svy,::svx],vecy[::svy,::svx],\
          angles=angle,color=color,pivot=pivot,\
          scale=scale,width=width,linewidths=linewidths,edgecolors=edgecolors)

          # make vector key.
          #keyh = 1.025 ; keyv = 1.05 # upper right corner over colorbar
          keyh = 0.97 ; keyv = 1.06
          keyh = 0.03 ; keyv = 1.07
          #keyh = -0.03 ; keyv = 1.08 # upper left corner
          labelpos='E'    # position label compared to arrow : N S E W
          p = mpl.quiverkey(q,keyh,keyv,\
          5.0,r'$5 m/s$',\
          fontproperties={'size': font,'weight': 'bold'},\
          color='black',labelpos=labelpos,labelsep = 0.07)

def getfigvar(nbfig,nbrow,nbcol,i):
    mpl.subplot(nbrow,nbcol,i+1)
    pal=get_cmap(name="PuBu")
    newlon=lon+180

    mymin=0.1
    mymax=50

    # log
    norm=mcolors.LogNorm()
    lvls=np.logspace(np.log10(mymin),np.log10(mymax),15)
    #titi=[1.e-14,1.e-13,1.e-12,1.e-11]
    CF=mpl.contourf(newlon, lat, cloud,levels=lvls,norm=norm,cmap=pal)
    '''
    # lin
    lev=np.linspace(mymin,mymax,6)
    CF=mpl.contourf(newlon, lat, cloud,lev,cmap=pal,extend='both')
    '''
    yticks=[-90,-60,-30,0,30,60,90]
    xticks=[0,60,120,180,240,300,360]
    cbar=mpl.colorbar(CF,shrink=1, format="%1.1f")
    cbar.ax.set_title("1E-6 kg/m2",y=1.04,fontsize=font)

    for t in cbar.ax.get_yticklabels():
        t.set_fontsize(font)

    vect=[47]
    print(('shape=',shape(newlon),shape(lat),shape(tmp)))
    CS=mpl.contour(newlon,lat,tmp,colors='k',linewidths=0.5)
    #### inline=1 : values over the line
    mpl.clabel(CS, inline=1, fontsize=20, fmt='%1.1f',inline_spacing=1)

    mpl.title('Local Time (180E) ='+str(tintstr[i]),fontsize=font)
    mpl.ylabel('Latitude (deg)',labelpad=10,fontsize=font)
    mpl.xlabel('Longitude (deg)',labelpad=10, fontsize=font)
    mpl.xticks(xticks,fontsize=font)
    mpl.yticks(yticks,fontsize=font)
    #getwinds(newlon,lat,u,v)

def getnumalt(choicealt,alt):
    numalt=np.where(abs(alt-choicealt)==min(abs(alt-choicealt)))
    numalt=numalt[0][0]
    return numalt

#######################

nbfig=size(tint)
nbrow=2
nbcol=2
mpl.figure(figsize=(30, 15))
choicealt=30 #m


for i in range(nbfig): 
   numalt=getnumalt(choicealt,alt)
   print((numalt,'alt=',alt[numalt],'m'))
   mycloud=getvar(filename1,var1,tint[i])[0,0,:,:]
   mytmp=getvar(filename1,var2,tint[i])[0,0,:,:]
   #myvar=getvar(filename2,var,tint[i])[0,0,:,:]

   cloud=swinglon(mycloud)
   cloud=cloud*1.e6
   tmp=swinglon(mytmp)
   #myvarbis=swinglon(myvar)

   getfigvar(nbfig,nbrow,nbcol,i)

left  = None  # the left side of the subplots of the figure
right = None   # the right side of the subplots of the figure
bottom = None  # the bottom of the subplots of the figure
top = None     # the top of the subplots of the figure
wspace = None  # the amount of width reserved for blank space between subplots
hspace = None  # the amount of height reserved for white space between subplots
#mpl.subplots_adjust(left, bottom, right, top, wspace, hspace)
#mpl.subplots_adjust(hspace = .1)

mpl.savefig('mapch4cloud.eps',dpi=200)
mpl.savefig('mapch4cloud.png',dpi=200)

#mpl.show()