#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ Created on Mon Oct 18 13:58:07 2021 @author: jene """ import glob # import os # import datetime as dt import xarray as xr # import netCDF4 as nc # from netCDF4 import MFDataset import numpy as np import pickle from matplotlib import pyplot as plt import matplotlib as mpl import matplotlib.path as mpath import cartopy.crs as ccrs import cartopy.feature as cfeature # import warnings # import scipy.io def cm_to_inch(value): return value/2.54 def calc_spatial_mean(xr_da, weights, lon_name="lon", lat_name="lat"): """ Calculate spatial mean of xarray.DataArray with grid cell weighting. Parameters ---------- xr_da: xarray.DataArray Data to average lon_name: str, optional Name of x-coordinate lat_name: str, optional Name of y-coordinate radius: float Radius of the planet [metres], currently assumed spherical (not important anyway) Returns ------- Spatially averaged xarray.DataArray. """ return (xr_da*weights).sum(dim=[lon_name, lat_name],skipna=True)/weights.sum(skipna=True) def calc_spatial_integral(xr_da, weights, lon_name="lon", lat_name="lat"): """ Calculate spatial integral of xarray.DataArray with grid cell weighting. Parameters ---------- xr_da: xarray.DataArray Data to average lon_name: str, optional Name of x-coordinate lat_name: str, optional Name of y-coordinate radius: float Radius of the planet [metres], currently assumed spherical (not important anyway) Returns ------- Spatially averaged xarray.DataArray. """ return (xr_da * weights).sum(dim=[lon_name, lat_name], skipna=True) def load_data(var, s): fname = glob.glob(dataPath + '' + var + '*' + s + '*.dat')[0] data = pickle.load(open(fname, 'rb')) grid = pickle.load(open(dataPath + 'grid.dat', 'rb')) mask = pickle.load(open(dataPath + 'mask.dat', 'rb')) grid = grid.where(((grid.lat>60) & (mask>0))) # grid.values[np.isnan(grid)] = 0 # attrs = data.attrs # data = xr.open_mfdataset('netCDF/' + var + '_' + s + '_*.nc') # if var == 'ph': # data = data['phlvl'] # else: # data = data[var] # data.attrs = attrs data = data.where(((data.lat>60) & (mask>0))) data = data.groupby('time.year').mean('time') # data = data.where(data.year<=2100, drop=True) if var == 'fgco2': ts = calc_spatial_integral(data, grid, 'x', 'y') ts = ts/1e9*(60*60*24*365.25) else: ts = calc_spatial_mean(data, grid, 'x', 'y') return data, ts def plot_map(ax, data, T): cmap = mpl.cm.Spectral_r if data.attrs['var'] == 'sst': bounds = np.arange(-1.5,1.51,.25) if data.attrs['var'] == 'sss': bounds = np.arange(-1,1.1,.2) if data.attrs['var'] == 'ph': bounds = np.arange(-.1,.11,.02) if data.attrs['var'] == 'fgco2': bounds = np.arange(-.5,.51,.1) # bounds = np.arange(-1.5,1.51,.25) data = data/(0.012)*(60*60*24*365.25) norm = mpl.colors.BoundaryNorm(bounds, cmap.N, extend='both') timeDiff = (data.where(data['year']>=2070 & (data['year']<=2100)).mean('year') - data.where((data['year']>=1985) & (data['year']<=2014)).mean('year')) im = ax.pcolormesh(data.lon, data.lat,timeDiff, transform=ccrs.PlateCarree(), cmap=cmap, norm=norm) ax.set_extent([-180, 180, 60, 90], ccrs.PlateCarree()) ax.add_feature(cfeature.LAND, color='lightgray', edgecolor='k') ax.coastlines() ax.gridlines(draw_labels=False, linewidth=.5, color='k', alpha=.5, linestyle=':') theta = np.linspace(0, 2*np.pi, 100) center, radius = [0.5, 0.5], 0.5 verts = np.vstack([np.sin(theta), np.cos(theta)]).T circle = mpath.Path(verts * radius + center) ax.set_boundary(circle, transform=ax.transAxes) ax.set_title(T, fontsize=fontsize['tit']) return im def plot_20yMap(ax, data, T): cmap = mpl.cm.Spectral_r if data.attrs['var'] == 'sst': bounds = np.arange(-2,8.1,1) if data.attrs['var'] == 'sss': bounds = np.arange(25,35.1,1) if data.attrs['var'] == 'ph': bounds = np.arange(7.5,8.01,.05) if data.attrs['var'] == 'fgco2': bounds = np.arange(-6,0.01,.6) data = data/(0.012)*(60*60*24*365.25) norm = mpl.colors.BoundaryNorm(bounds, cmap.N, extend='both') im = ax.pcolormesh(data.lon, data.lat, data.where(data['year']>=2080).mean('year'), transform=ccrs.PlateCarree(), cmap=cmap, norm=norm ) ax.set_extent([-180, 180, 60, 90], ccrs.PlateCarree()) ax.add_feature(cfeature.LAND, color='lightgray', edgecolor='k') ax.coastlines() ax.gridlines(draw_labels=False, linewidth=.5, color='k', alpha=.5, linestyle=':') theta = np.linspace(0, 2*np.pi, 100) center, radius = [0.5, 0.5], 0.5 verts = np.vstack([np.sin(theta), np.cos(theta)]).T circle = mpath.Path(verts * radius + center) ax.set_boundary(circle, transform=ax.transAxes) ax.set_title(T, fontsize=fontsize['tit']) return im fontsize = {'supt':18, 'tit':16, 'axtit':13, 'axtick':11} figNum =['a) ','b) ','c) ','d) ','e) ','f) ','g) ','h) ','i) ','j) ','k) ', 'l) ','m) ','n) ','o) ','p) ','q) ','r) ','s) ','t) ','u) ','v) ', 'w) ','x) ','y) ','z) ', r'$\alpha$)', r'$\beta)$', r'$\gamma)$', r'$\delta)$'] figExt = '.png' basePath = '/media/ipcc/' dataPath = basePath + 'NS9252K/KeyCLIM_CMOR/' dataPath = basePath + 'NS9252K/cmorout/' dataPath = 'Data/' refDataPath = basePath + 'NS9034K/CMIP6/' variable = [ 'sst', 'fgco2', 'sss', # 'ph' ] scenario = [ 'cloud2', 'eddy', 'iceSheet3', 'snow' ] grid = pickle.load(open(dataPath + 'grid.dat', 'rb')) ############################################################################### ### Timeseries ############################################################################### idx = 1 fig = plt.figure(figsize=(cm_to_inch(29.7),cm_to_inch(21)), dpi=300) for var in variable: dataRef, tsRef = load_data(var, 'ref') ax = fig.add_subplot(2,2,idx) ax.plot(tsRef.year, tsRef, label='Reference') for s in scenario: print('Processing ' + var + ', ' + s) data, ts = load_data(var, s) if s == 'ref': ls = 'Reference' if s == 'cloud2': ls = 'Cloud' if s == 'iceSheet3': ls = 'IceSheet' if s == 'snow': ls = 'Snow' if s == 'eddy': ls = 'Eddy' # if data.attrs['var'] == 'fgco2': # ts = ts/1e9*(60*60*24*365.25) # tsRef = tsRef/1e9*(60*60*24*365.25) ax.plot(ts.year, ts, label=ls) if data.attrs['var'] == 'fgco2': ax.set_ylabel('CO$_2$ flux [Tg yr${-1}$]', fontsize=fontsize['axtit']) T = 'CO$_2$ flux' if data.attrs['var'] == 'sst': ax.set_ylabel('SST [degC]', fontsize=fontsize['axtit']) T = 'SST' if data.attrs['var'] == 'sss': ax.set_ylabel('SSS [psu]', fontsize=fontsize['axtit']) T = 'SSS' if data.attrs['var'] == 'ph': ax.set_ylabel('pH', fontsize=fontsize['axtit']) T = 'pH' ax.legend() ax.set_xlim([1850,2100]) ax.grid(linestyle=':') ax.set_title(figNum[idx-1] + ' ' + T, fontsize=fontsize['tit']) if idx> 2: ax.set_xlabel('Year',fontsize=fontsize['axtit']) idx+=1 plt.subplots_adjust(wspace=.2,hspace=.3) plt.savefig('Plots/timeseries' + figExt) plt.show() plt.close() # ############################################################################### # ### Maps # ############################################################################### for var in variable: fig = plt.figure(figsize=(cm_to_inch(29.7),cm_to_inch(21)), dpi=300) gs = fig.add_gridspec(2,6) dataRef, tsRef = load_data(var, 'ref') # ax1 = fig.add_subplot(3,2,1) im = plot_map(fig.add_subplot(gs[0, 0:2],projection=ccrs.NorthPolarStereo()), dataRef, figNum[0] + ' Reference') idx = 1 for s in scenario: print('Processing ' + var + ', ' + s) data, ts = load_data(var, s) if 'ref' in s : ls = 'Reference' if s == 'cloud2': ls = 'Cloud' if s == 'iceSheet3': ls = 'IceSheet' if s == 'snow': ls = 'Snow' if s == 'eddy': ls = 'Eddy' if idx>2: plot_map(fig.add_subplot(gs[1,(idx-3)*3:(idx-2)*3], projection=ccrs.NorthPolarStereo()), data, figNum[idx] + ' ' + ls) else: plot_map(fig.add_subplot(gs[0,(idx)*2:(idx+1)*2], projection=ccrs.NorthPolarStereo()), data, figNum[idx] + ' ' + ls) # cbar.set_label('mean(' + variable + ') [' + unit + ']', fontsize=fontsize['axtit'])#, rotation=270, labelpad=+10) # cbar.ax.tick_params(labelsize=fontsize['axtick']) # plt.savefig('Plot/mean_' + Z[3:] + '.png') idx+=1 if data.attrs['var'] == 'fgco2': # fig.suptitle(data.attrs['var'], # fontsize=fontsize['supt'],fontweight="bold") fig.suptitle('$\Delta$CO$_2$ flux', fontsize=fontsize['supt'],fontweight="bold", y=0.95) cbarT = 'CO$_2$ flux [mol C m${-2}$ yr${-1}$]' if data.attrs['var'] == 'sst': fig.suptitle('$\Delta$SST', fontsize=fontsize['supt'],fontweight="bold", y=0.95) cbarT = 'SST [degC]' if data.attrs['var'] == 'sss': fig.suptitle('$\Delta$SSS', fontsize=fontsize['supt'],fontweight="bold", y=0.95) cbarT = 'SSS [psu]' if data.attrs['var'] == 'ph': fig.suptitle('$\Delta$pH', fontsize=fontsize['supt'],fontweight="bold", y=0.95) cbarT = 'pH' cbar_ax = fig.add_axes([0.25, 0.05, 0.5, 0.025]) cbar = fig.colorbar(im, orientation="horizontal", cax=cbar_ax) cbar.set_label('$\Delta$' + cbarT, fontsize=fontsize['axtit'])#, rotation=270, labelpad=+10) cbar.ax.tick_params(labelsize=fontsize['axtick']) plt.subplots_adjust(wspace=.1, hspace=.1) plt.savefig('Plots/Maps_' + var + figExt, bbox_inches='tight') plt.show() plt.close() for var in variable: fig = plt.figure(figsize=(cm_to_inch(29.7),cm_to_inch(21)), dpi=300) gs = fig.add_gridspec(2,6) dataRef, tsRef = load_data(var, 'ref') # ax1 = fig.add_subplot(3,2,1) im = plot_20yMap(fig.add_subplot(gs[0, 0:2],projection=ccrs.NorthPolarStereo()), dataRef, figNum[0] + ' Reference') idx = 1 for s in scenario: print('Processing ' + var + ', ' + s) data, ts = load_data(var, s) # attrs = data.attrs # data = data+dataRef # data.attrs = attrs if 'ref' in s : ls = 'Reference' if s == 'cloud2': ls = 'Cloud' if s == 'iceSheet3': ls = 'IceSheet' if s == 'snow': ls = 'Snow' if s == 'eddy': ls = 'Eddy' if idx>2: plot_20yMap(fig.add_subplot(gs[1,(idx-3)*3:(idx-2)*3], projection=ccrs.NorthPolarStereo()), data, figNum[idx] + ' ' + ls) else: plot_20yMap(fig.add_subplot(gs[0,(idx)*2:(idx+1)*2], projection=ccrs.NorthPolarStereo()), data, figNum[idx] + ' ' + ls) # cbar.set_label('mean(' + variable + ') [' + unit + ']', fontsize=fontsize['axtit'])#, rotation=270, labelpad=+10) # cbar.ax.tick_params(labelsize=fontsize['axtick']) # plt.savefig('Plot/mean_' + Z[3:] + '.png') idx+=1 if data.attrs['var'] == 'fgco2': # fig.suptitle(data.attrs['var'], # fontsize=fontsize['supt'],fontweight="bold") fig.suptitle('CO$_2$ flux (2080-2100)', fontsize=fontsize['supt'],fontweight="bold", y=0.95) cbarT = 'CO$_2$ flux [mol C m-2 yr-1]' if data.attrs['var'] == 'sst': fig.suptitle('SST (2080-2100)', fontsize=fontsize['supt'],fontweight="bold", y=0.95) cbarT = 'SST [degC]' if data.attrs['var'] == 'sss': fig.suptitle('SSS (2080-2100)', fontsize=fontsize['supt'],fontweight="bold", y=0.95) cbarT = 'SSS [psu]' if data.attrs['var'] == 'ph': fig.suptitle('pH (2080-2100)', fontsize=fontsize['supt'],fontweight="bold", y=0.95) cbarT = 'pH' cbar_ax = fig.add_axes([0.25, 0.05, 0.5, 0.025]) cbar = fig.colorbar(im, orientation="horizontal", cax=cbar_ax) cbar.set_label('' + cbarT, fontsize=fontsize['axtit'])#, rotation=270, labelpad=+10) cbar.ax.tick_params(labelsize=fontsize['axtick']) plt.subplots_adjust(wspace=.1, hspace=.1) plt.savefig('Plots/Maps_2080-2100_' + var + figExt, bbox_inches='tight') plt.show() plt.close()