''' Created on 21/02/2017 @author: ramato ''' import iris import calendar import sys import numpy as np import math import iris.plot as iplt import matplotlib.pyplot as plt import matplotlib.colors as mcols import cartopy.crs as ccrs # coordinate reference system iris.FUTURE.netcdf_promote = True DIR = '/home/h03/ramato/CSSP_validation' RID = 'apepb-apepc' months = calendar.month_abbr def rounddown(x): if x>10: return int(math.floor(x / 10.0)) * 10 else: return math.floor(x) def load_file(): fil = '{}/data/bigendian/{}.h9i5*.05216.be_norim.pp'.format(DIR, RID) cube = iris.load_cube(fil) return cube def mmday(cube): new_cube = cube * 86400 new_cube.units = 'mm/day' new_cube.standard_name = 'precipitation_flux' print new_cube return new_cube def get_mon(i): if i>0: mon = months[i] elif i == 0: mon = months[12]#dec return mon def plot_subplot(cube, i, mycmap, mymax, mon): ax1 = plt.subplot(3, 4, i, projection=ccrs.PlateCarree()) ax1.coastlines() #ax1.set_extent([-45, 55, 20, 80], ccrs.PlateCarree()) mesh = iplt.contourf(cube, cmap=mycmap, levels=np.arange(-mymax, mymax+1)) plt.title('{}'.format(mon)) add_gridlines(ax1) return mesh def make_colourbar(mesh, myextend, tick_levels, cbar_title): '''standard setup for plot colourbar''' cax = plt.axes([0.1, 0.1, 0.8, 0.025]) cbar = plt.colorbar(mesh, cax=cax, orientation='horizontal', extend=myextend) cbar.set_ticks(tick_levels) cbar.set_label(cbar_title) def add_gridlines(axs): '''Add gridlines, showing lat-lon values''' gls = axs.gridlines(draw_labels=True) gls.xlabels_top = None gls.ylabels_right = None def plot_log_precip(cube, i): #fig = plt.figure() #fig.text(0.5, 0.97, '{} precip ({})'.format(RID, mon), # horizontalalignment = 'center', verticalalignment = 'top', fontsize = 15, fontweight = 'bold') #for i, cube in enumerate(cubes[:], 0): # print cube.coord('time') # fig.add_subplot(2, 2, i+1) mon = get_mon(i) plt.title('{}'.format(mon)) # Define scaling levels for the logarithmic colouring. minimum_log_level = 0.1 maximum_scale_level = max(np.max(cube.data), abs(np.min(cube.data))) mymax = rounddown(maximum_scale_level) # Create a 'logarithmic' data normalization. anom_norm = mcols.SymLogNorm(linthresh=minimum_log_level, linscale=0, vmin=-maximum_scale_level, vmax=maximum_scale_level) # Setting "linthresh=minimum_log_level" makes its non-logarithmic # data range equal to our 'zero band'. # Setting "linscale=0" maps the whole zero band to the middle colour value # (i.e. 0.5), which is the neutral point of a "diverging" style colormap. ax = plt.axes(projection=ccrs.PlateCarree()) mesh = iplt.pcolormesh(cube, cmap=plt.get_cmap('Spectral'), norm=anom_norm)#bwr_r cax = plt.axes([0.1, 0.1, 0.8, 0.025]) cbar = plt.colorbar(mesh, cax=cax, orientation='horizontal', extend = 'both') cbar.set_label('logarithmic precipitation change (mm/day)') # Set some suitable fixed "logarithmic" colourbar tick positions. tick_levels = [-mymax, 0.0, mymax] #tick_levels = [-mymax, -10, -1, -0.3, 0.0, 0.3, 1, 10, mymax] cbar.set_ticks(tick_levels) # Modify the tick labels so that the centre one shows "+/-". tick_levels[1] = r'$\pm${:g}'.format(minimum_log_level) cbar.set_ticklabels(tick_levels) ax.coastlines() gl = ax.gridlines(draw_labels=True) gl.xlabels_top = None gl.ylabels_right = None #plt.show() def make_subplot(cube, plot_type): maximum_scale_level = max(np.max(cube.data), abs(np.min(cube.data))) mymax = rounddown(maximum_scale_level) if plot_type == 'difference': mycmap = plt.get_cmap('Spectral') mymin = -mymax tick_levels = np.linspace(mymin, mymax, 9) myextend = 'both' elif plot_type == 'raw data': mycmap = plt.get_cmap('YlGnBu') mymin = 0 tick_levels = np.linspace(mymin, mymax, 5) myextend = 'max' fig = plt.figure(figsize=(12, 10)) fig.text(0.5, 0.97, '{} monthly precipitation'.format(RID), horizontalalignment = 'center', verticalalignment = 'top', fontsize = 15, fontweight = 'bold') cbar_title = 'linear precipitation (mm/day)' for i in range(0, 2):#cube.shape[0]): mon = get_mon(i) print i, mon mesh = plot_subplot(cube[i], i, mycmap, mymax, mon) make_colourbar(mesh, myextend, tick_levels, cbar_title) plt.show() def main(): cube = load_file() cube_mm = mmday(cube) make_subplot(cube_mm, 'difference') #print cube_mm #[0==dec] #plot_log_precip(cube_mm) #plot_precip(cube_mm, mon) if __name__ == "__main__": main()