#! /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=getvar(nc1,"latitude") lon=getvar(nc1,"longitude") # 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()