[3823] | 1 | #! /usr/bin/env python |
---|
| 2 | from netCDF4 import Dataset |
---|
| 3 | from numpy import * |
---|
| 4 | import numpy as np |
---|
| 5 | import matplotlib.pyplot as mpl |
---|
| 6 | from matplotlib.cm import get_cmap |
---|
| 7 | import pylab |
---|
| 8 | from matplotlib import ticker |
---|
| 9 | import matplotlib.colors as mcolors |
---|
| 10 | from mpl_toolkits.basemap import Basemap |
---|
[3833] | 11 | from FV3_utils import * |
---|
[3823] | 12 | |
---|
| 13 | ############################ |
---|
[3833] | 14 | filename=name+".nc" |
---|
[3823] | 15 | nb_dataset=1 #change vect as well |
---|
[3833] | 16 | tint=[0,0.001] #Time must be as written in the input file |
---|
[3823] | 17 | #zint=["0"] #alt in km |
---|
| 18 | #xarea="-180,179" |
---|
| 19 | #yarea="-90,90" |
---|
| 20 | #step=5 #step between each bin |
---|
[3833] | 21 | var="ALB" #variable |
---|
[3823] | 22 | var2="phisinit" |
---|
| 23 | ############################ |
---|
| 24 | |
---|
| 25 | nc=Dataset(filename) |
---|
[3833] | 26 | lat=getvar(nc,"latitude") |
---|
| 27 | lon=getvar(nc,"longitude") |
---|
| 28 | tim=getvar(nc,"Time") |
---|
[3823] | 29 | |
---|
[3833] | 30 | myvar = getvar(nc,var,tint,t_mean=True) |
---|
| 31 | myvar2 = getvar(nc,var2) |
---|
[3823] | 32 | print(('shape var1 =',shape(myvar))) |
---|
| 33 | print(('shape lat = ',shape(lat))) |
---|
| 34 | print(('lat = ',lat)) |
---|
| 35 | |
---|
| 36 | |
---|
| 37 | myvar=myvar |
---|
| 38 | myvar2=myvar2/0.6192/1000. # geopotential ==> altitude |
---|
| 39 | |
---|
| 40 | ################ |
---|
| 41 | # colorbar |
---|
| 42 | def make_colormap(seq): |
---|
| 43 | """Return a LinearSegmentedColormap |
---|
| 44 | seq: a sequence of floats and RGB-tuples. The floats should be increasing |
---|
| 45 | and in the interval (0,1). |
---|
| 46 | """ |
---|
| 47 | |
---|
| 48 | seq = [(None,) * 3, 0.0] + list(seq) + [1.0, (None,) * 3] |
---|
| 49 | print(seq) |
---|
| 50 | cdict = {'red': [], 'green': [], 'blue': []} |
---|
| 51 | for i, item in enumerate(seq): |
---|
| 52 | if isinstance(item, float): |
---|
| 53 | r1, g1, b1 = seq[i - 1] |
---|
| 54 | r2, g2, b2 = seq[i + 1] |
---|
| 55 | cdict['red'].append([item, r1, r2]) |
---|
| 56 | cdict['green'].append([item, g1, g2]) |
---|
| 57 | cdict['blue'].append([item, b1, b2]) |
---|
| 58 | print(cdict) |
---|
| 59 | return mcolors.LinearSegmentedColormap('CustomMap', cdict) |
---|
| 60 | ###### |
---|
| 61 | c = mcolors.ColorConverter().to_rgb |
---|
| 62 | #rvb1 = make_colormap([c('red'), c('violet'), 0.33, c('violet'), c('blue'), 0.66, c('blue')]) |
---|
| 63 | rvb1 = make_colormap([c('SaddleBrown'), c('violet'), 0.33, c('violet'), c('white'), 0.66, c('white')]) |
---|
| 64 | |
---|
| 65 | pp=(0.5450980392156862-0.15, 0.27058823529411763-0.05*3, 0.07450980392156863-0.01) |
---|
| 66 | rvb1 = make_colormap([pp,0.33, c('indianred'),0.71, c('lightcyan'),0.99,c('lightsalmon')]) |
---|
| 67 | #rvb1 = make_colormap([c('saddlebrown'),0.33, c('indianred'),0.66, c('white')]) |
---|
| 68 | ###### |
---|
| 69 | |
---|
| 70 | div=3. |
---|
| 71 | nbticksbar=3. |
---|
| 72 | #pal=get_cmap(name="Blues") |
---|
| 73 | pal=get_cmap(name="gist_earth_r") |
---|
| 74 | yticks=[-90,-60,-30,0,30,60,90] |
---|
| 75 | #xticks=[-180,-120,-60,0,60,120,180] |
---|
| 76 | xticks=[0,60,120,180,240,300,360] |
---|
| 77 | |
---|
| 78 | # changer les longitudes pour mettre TR au centre |
---|
| 79 | vec=shape(myvar) |
---|
| 80 | myvarbis=np.zeros(vec,dtype='f') |
---|
| 81 | myvar2bis=np.zeros(vec,dtype='f') |
---|
| 82 | # i lat : pas de changement |
---|
| 83 | # j lon : |
---|
| 84 | for i in range(vec[0]): |
---|
| 85 | for j in range(vec[1]): |
---|
| 86 | if j < int(vec[1]/2.) : |
---|
| 87 | myvarbis[i,j]=myvar[i,j+int(vec[1]/2)] |
---|
| 88 | myvar2bis[i,j]=myvar2[i,j+int(vec[1]/2)] |
---|
| 89 | else: |
---|
| 90 | myvarbis[i,j]=myvar[i,j-int(vec[1]/2)] |
---|
| 91 | myvar2bis[i,j]=myvar2[i,j-int(vec[1]/2)] |
---|
| 92 | |
---|
[3833] | 93 | if lon[0]<0: |
---|
| 94 | lon=lon+180. |
---|
[3823] | 95 | |
---|
| 96 | |
---|
| 97 | mpl.figure(figsize=(17, 10)) |
---|
| 98 | |
---|
| 99 | c=mpl.contour(lon, lat, myvar2bis, 10,levels=[-3,-2,-1,-0.5,-0.2], colors = 'k', linewidths = 0.5) |
---|
| 100 | mpl.clabel(c, fmt='%2.1f', colors='k', fontsize=14) |
---|
| 101 | |
---|
| 102 | mpl.imshow(myvarbis,extent=[0,360,-90,90],cmap=rvb1) #,extend='both') |
---|
| 103 | |
---|
| 104 | #mpl.colorbar() |
---|
| 105 | mpl.ylabel('Latitude (deg)',labelpad=20,fontsize=29) |
---|
| 106 | mpl.xlabel('Longitude (deg)',labelpad=20, fontsize=29) |
---|
| 107 | mpl.xticks(xticks,fontsize=26) |
---|
| 108 | mpl.yticks(yticks,fontsize=26) |
---|
| 109 | mpl.grid() |
---|
| 110 | mpl.savefig('mapalbini.eps',dpi=200) |
---|
| 111 | mpl.savefig('mapalbini.png',dpi=200) |
---|
| 112 | mpl.show() |
---|
| 113 | |
---|
| 114 | |
---|