'''Plotting script for Anglo American Creates fut - his seasonal plots for Consecutive dry day data Plots country borders and a point over Santiago. Need to change the load_data() path to make the script more generic for different variables. Rosanna Amato, 09/02/2017''' import iris import iris.plot as iplt import calendar import numpy as np import matplotlib.pyplot as plt import matplotlib.ticker as mticker import cartopy.crs as ccrs # coordinate reference system import cartopy.feature as cfeature iris.FUTURE.netcdf_promote = True import sys #pylint: disable=no-member #pylint: disable=trailing-whitespace ######################################################### #MODEL = 'HadGEM2-ES' #HIS = 'mi-ao554' # #FUT = 'mi-ap163' # MODEL = 'MPI-ESM-LR' HIS = 'mi-ap448' # FUT = 'mi-ap430' # YR1 = 1985 YR2 = 2005 YR3 = 2035 YR4 = 2065 Santiago = [-70.67, -33.45] DIR = '/project/ccrs/ramato/monitoring/anglo' ######################################################### VAR = 'precipitation_flux' MYMAP = 'Spectral' STASH = '05216' #VAR = 'CDD' #LONG_VAR = 'Consecutive dry days' #MYMAP = 'Spectral_r' #STASH = '' # INFIL = '{}/{}/output/monthly/data/CDD/Santiago_CDD_monthly_max_c_mean.nc'.format(DIR, rid) #VAR = 'air_temperature' #MYMAP = 'Spectral_r' #STASH = '03236.128' ########################################################### SEASONS = {'djf' : ['Dec', 'Jan', 'Feb'], 'mam' : ['Mar', 'Apr', 'May'], 'jja' : ['Jun', 'Jul', 'Aug'], 'son' : ['Sep', 'Oct', 'Nov'], } NROW = 1 # len(SEASONS) NCOL = 3 # fut, his, diff #DATAMAX = 28.0 #DIFFMAX = 6.0 MONTHS = calendar.month_abbr abbr_to_num = {name: num for num, name in enumerate(calendar.month_abbr) if num} ############################################################ def load_data(rid): cube_data = '{}/{}/output/monthly/data/CDD/Santiago_CDD_monthly_max_c_mean.nc'.format(DIR, rid) cube = iris.load_cube(cube_data) if cube.standard_name == 'precipitation_flux': cube = mmday(cube) elif cube.units == 'K': cube.convert_units('celsius') print cube[0].coord('time') return cube def mmday(cube): '''Convert precip data from kg/m/s to mm/day''' new_cube = cube * 86400 new_cube.units = 'mm/day' new_cube.standard_name = 'precipitation_flux' return new_cube def season_average(cube, key): ''' Create a three month seasonal average cube (flattens time so 2D for plotting). This uses the dictionary 'SEASONS' which contains the names of the months for each season.''' #djf_cube = cube.extract(iris.Constraint(month=[12, 1, 2])) seas = cube.extract(iris.Constraint(month=[ abbr_to_num[SEASONS.get(key)[0]], abbr_to_num[SEASONS.get(key)[1]], abbr_to_num[SEASONS.get(key)[2]] ])) seas_av = seas.collapsed('time', iris.analysis.MEAN) return seas_av def calc_diff(his, fut, key): ''' calc future - baseline values for CSSP validation''' diff = fut - his diff.long_name = '{} {} ({} difference)'.format(MODEL, VAR, key) return diff def get_lims(cube): mymin = np.min(cube.data) mymax = np.max(cube.data) abs_min = round(abs(mymin)) abs_max = round(abs(mymax)) print mymin, mymax print abs_min, abs_max mylim = max(abs_min, abs_max) #print mylim return mylim def add_gridlines(axs): '''Add coastlines and gridlines, showing lat-lon values''' axs.add_feature(cfeature.OCEAN) axs.add_feature(cfeature.COASTLINE) axs.add_feature(cfeature.BORDERS) #gls = axs.gridlines(draw_labels=True) #gls.xlabels_top = None #gls.ylabels_right = None def plot_subplot(cube, i, mycmap, myrange, sp_title): ax1 = plt.subplot(NROW, NCOL, i, projection=ccrs.PlateCarree()) ax1.coastlines() #ax1.set_extent([-45, 55, 20, 80], ccrs.PlateCarree()) mesh = iplt.contourf(cube, cmap=mycmap, levels=myrange) ax1.plot(Santiago[0], Santiago[1], 'ro', markersize=6, transform=ccrs.PlateCarree()) plt.title('{}'.format(sp_title)) add_gridlines(ax1) return mesh def make_colourbar(cube, mesh, myextend): '''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('{} ({})'.format(VAR, cube.units)) def plot_diff(his, fut, diff, mylim, key): fig = plt.figure(figsize=(12, 6)) fig.text(0.5, 0.97, '{}: {} in {}'.format(MODEL, VAR, key.upper()), horizontalalignment = 'center', verticalalignment = 'top', fontsize = 15, fontweight = 'bold') #plot_subplot(his, 1, plt.get_cmap('YlGnBu_r'), np.arange(0, DATAMAX+1), 'Historic \n ({}:{})'.format(YR1,YR2)) plot_subplot(his, 1, plt.get_cmap('YlGnBu_r'), np.arange(0, np.max(his.data)), 'Historic \n ({}:{})'.format(YR1,YR2)) #mesh = plot_subplot(fut, 2, plt.get_cmap('YlGnBu_r'), np.arange(0, DATAMAX+1), 'RCP8.5 Future \n ({}:{})'.format(YR3,YR4)) mesh = plot_subplot(fut, 2, plt.get_cmap('YlGnBu_r'), np.arange(0, np.max(fut.data)), 'RCP8.5 Future \n ({}:{})'.format(YR3,YR4)) cax1 = plt.axes([0.02, 0.1, 0.64, 0.025]) #[ from left, from bottom, length, thickness] cbar1 = plt.colorbar(mesh, cax=cax1, orientation='horizontal', extend='max') cbar1.set_label = '{} ({})'.format(VAR, his.units) diff_range = np.arange(-mylim, mylim+1) #diff_range = np.arange(-DIFFMAX, DIFFMAX+1) mesh = plot_subplot(diff, 3, plt.get_cmap(MYMAP), diff_range, 'Difference \n (future - historic)') cax2 = plt.axes([0.68, 0.1, 0.3, 0.025]) #[ from left, from bottom, length, thickness] cbar2 = plt.colorbar(mesh, cax=cax2, orientation='horizontal', extend='both') cbar2.set_label = '{} difference ({})'.format(VAR, diff.units) plt.tight_layout() #plt.savefig('{}/plots/{}_{}_{}_diff.png'.format(DIR, key, MODEL, VAR)) # plt.show() ################################################################################ def main(): ''' plot future - baseline CDD values for Anglo American project''' all_his = load_data(HIS) all_fut = load_data(FUT) for key in SEASONS: print key his = season_average(all_his, key) fut = season_average(all_fut, key) diff = calc_diff(his, fut, key) print diff lim = get_lims(diff) plot_diff(his, fut, diff, lim, key) plt.show() sys.exit() if __name__ == '__main__': main()