################################################### ### Preparations (load modules, definitions...) ################################################### import matplotlib import matplotlib.pyplot as plt from matplotlib.gridspec import GridSpec import matplotlib._color_data as mcd from mpl_toolkits.axes_grid1.inset_locator import inset_axes, zoomed_inset_axes, mark_inset import numpy as np import math #as ma import scipy from scipy import stats, ndimage, interpolate import pandas as pd import netCDF4 as nc4 from netCDF4 import Dataset from netCDF4 import num2date import pickle import datetime import glob import os import time import mpl_toolkits.basemap from mpl_toolkits.basemap import Basemap, addcyclic, shiftgrid import pyferret pyferret.start(quiet=True) import importlib #my scripts: import tools_plot as jplt importlib.reload(jplt) import tools_analysis as jan importlib.reload(jan) def var_pad(var): var_p=np.pad(np.pad(var, ((1,1),(0,0)), mode='edge'), ((0,0),(1,1)), mode='wrap') return var_p ### read out biogeo_df from all_snapshots_mom_annual.nc paths=sorted(glob.glob('../c3a_model_input_output/*Ma*')) paths_1000ppm=sorted(glob.glob('../c3a_model_input_output/*Ma_1000ppm')) paths_1000ppm_S0ini=sorted(glob.glob('../c3a_model_input_output/*Ma_1000ppm_S0ini')) paths_1000ppm_S0ini_VegFix=sorted(glob.glob('../c3a_model_input_output/*Ma_1000ppm_S0ini_VegFix')) paths_1000ppm_VegHom=sorted(glob.glob('../c3a_model_input_output/*Ma_1000ppm_VegHom')) #paths_1000ppm_VegTree=sorted(glob.glob('../c3a_model_input_output/*Ma_1000ppm_VegTree')) #paths_1000ppm_VegBare=sorted(glob.glob('../c3a_model_input_output/*Ma_1000ppm_VegBare')) paths_orb=sorted(glob.glob('../c3a_model_input_output/*Ma*orb*')) paths_1000ppm_orb=sorted(glob.glob('../c3a_model_input_output/*Ma_1000ppm_orb_*')) paths_1000ppm_orb_22p0=sorted(glob.glob('../c3a_model_input_output/*Ma_1000ppm_orb_22p0-*')) paths_1000ppm_orb_24p5=sorted(glob.glob('../c3a_model_input_output/*Ma_1000ppm_orb_24p5-*')) for ii in paths_1000ppm_orb: paths_orb.remove(ii) paths_SL=sorted(glob.glob('../c3a_model_input_output/*Ma_*ppm_SL*')) paths_VegHom_SL=sorted(glob.glob('../c3a_model_input_output/*Ma_1000ppm_VegHom_SL*')) #paths_0Ma=sorted(glob.glob('../c3a_model_input_output/000Ma_*')) paths_225Ma=sorted(glob.glob('../c3a_model_input_output/225Ma_*')) paths_250ppm=[] paths_500ppm=[] paths_1500ppm=[] paths_2000ppm=[] paths_4000ppm=[] paths_1000ppm_25Ma=[] timeslices_25Ma=np.append([0,65],np.arange(75,250+25,25)) for tt in timeslices_25Ma[1:]: age_str='{:03d}'.format(tt) paths_250ppm.append(glob.glob('../c3a_model_input_output/'+age_str+'Ma_250ppm')[0]) paths_500ppm.append(glob.glob('../c3a_model_input_output/'+age_str+'Ma_500ppm')[0]) paths_1500ppm.append(glob.glob('../c3a_model_input_output/'+age_str+'Ma_1500ppm')[0]) paths_2000ppm.append(glob.glob('../c3a_model_input_output/'+age_str+'Ma_2000ppm')[0]) paths_4000ppm.append(glob.glob('../c3a_model_input_output/'+age_str+'Ma_4000ppm')[0]) paths_1000ppm_25Ma.append(glob.glob('../c3a_model_input_output/'+age_str+'Ma_1000ppm')[0]) pCO2_values_df=pd.read_pickle('./meso_BoundaryConditions_pCO2.pkl') pCO2_values_df paths_pCO2Foster=[] paths_pCO2COPSE=[] for age in pCO2_values_df.index[1:]: age_str='{:03d}'.format(age) [pCO2_Foster,pCO2_Mills]=pCO2_values_df.loc[age,['pCO2_Foster_round','pCO2_Mills_round']] paths_pCO2Foster.append('../c3a_model_input_output/'+age_str+'Ma_'+str(int(pCO2_Foster))+'ppm') paths_pCO2COPSE.append('../c3a_model_input_output/'+age_str+'Ma_'+str(int(pCO2_Mills))+'ppm') paths_problematic=[] timeslices_problematic=[] timeslices_failed=[]+timeslices_problematic paths_failed=[] # for tt in timeslices_problematic: # tt_string='{:03d}'.format(tt) # paths_problematic.append('../MesoClim/c3a_meso_'+tt_string+'Ma') # for tt in timeslices_failed: # tt_string='{:03d}'.format(tt) # paths_failed.append('../MesoClim/c3a_meso_'+tt_string+'Ma') # paths_problematic.append('../c3a_model_input_output/000Ma_280ppm_ice_runoffuni') # paths_problematic.append('../c3a_model_input_output/070Ma_1000ppm_runoff') # paths_problematic.append('../c3a_model_input_output/235Ma_1000ppm_stabilize_runoff') # paths_problematic.append('../c3a_model_input_output/235Ma_1000ppm_stabilize_topo') # paths_problematic.append('../c3a_model_input_output/000Ma_280ppm_ice_NoRunoffTuning') paths_paper_example=['../c3a_model_input_output/075Ma_1000ppm', '../c3a_model_input_output/125Ma_1000ppm', '../c3a_model_input_output/175Ma_1000ppm', '../c3a_model_input_output/225Ma_1000ppm'] # paths_paper_appendix=paths[:] # for pp in paths_SL+paths_0Ma+paths_orb+paths_1000ppm_VegTree+paths_1000ppm_VegBare+paths_problematic: # if pp in paths_paper_appendix: # paths_paper_appendix.remove(pp) # paths_paper=paths[:] # for pp in paths_0Ma+paths_orb+paths_1000ppm_VegTree+paths_1000ppm_VegBare+paths_problematic: # if pp in paths_paper: # paths_paper.remove(pp) ###Colour list color_dict=mcd.TABLEAU_COLORS #XKCD_COLORS color_list=list(color_dict.keys())[:]*30 colors_df=pd.DataFrame([], columns=['color'], index=paths) colors_df.iloc[:,0]=color_list[0:len(paths)] color_trias=np.array([129.,43.,146.])/255 color_jura=np.array([52.,178.,201.])/255 color_cret=np.array([127.,198.,78.])/255 color_trias_lower=np.array([168.,88.,158.])/255 color_trias_middle=np.array([188.,134.,186.])/255 color_trias_upper=np.array([199.,166.,207.])/255 color_jura_lower=np.array([0.,176.,233.])/255 color_jura_middle=np.array([132.,207.,237.])/255 color_jura_upper=np.array([188.,228.,250.])/255 color_cret_lower=np.array([160.,201.,109.])/255 color_cret_upper=np.array([186.,210.,95.])/255 # see Stadler: Aridity Indexes in Encyclopedia of World Climatology 2005 #https://en.wikipedia.org/wiki/Aridity_index def func_aridity_idx(prc,radbal): L_wv=2264.705 #Latent heat of water vaporization, kJ/kg, #radbal W/m^2=(J/s)/m^2; prc mm/d~kg/d annual mean aridity_budyko=100*radbal/(L_wv*prc) # (J/s)/m^2 * 1/(kJ/kg) * 1/(kg/d/m^2) = d/(1000*s) = 24*60*60/1000 = 86.4 # radbal/L_wv is potential evaporation (Arora 2002) return aridity_budyko latsop,lonsop,lonsop_2D,latsop_2D=jan.grid_gen(region='ocean', cell='p')#[0:1+1] latsot,lonsot,lonsot_2D,latsot_2D=jan.grid_gen(region='ocean', cell='T')#[0:1+1] latsap,lonsap,lonsap_2D,latsap_2D=jan.grid_gen(region='atmo', cell='p')#[0:1+1] latsat,lonsat,lonsat_2D,latsat_2D=jan.grid_gen(region='atmo', cell='T')#[0:1+1] nc_fid = Dataset('./gridarea_mom.nc', 'r') #nc_attrs, nc_dims, nc_vars = ncdump(nc_fid) gridcell_areas_mom=nc_fid.variables['cell_area'][:]/1000**2/1000**2 #in the file its in m**2, convert it to 10^6 km^2 timeslices=np.append(0,np.arange(60,255+5,5)) timeslices_25Ma=np.append(np.array([0,65]),timeslices[4::5]) timeslices_SL=[100,150,200] pi_frac=(np.pi/180)*latsot[:] timeslices_pcolor=np.append(50-2.5,np.arange(60-2.5,255+2.5+5,5)) timeslices_plot=np.append(50,np.arange(60,255+5,5)) timeslices_25Ma_plot=np.append([50,65],np.arange(75,250+25,25)) timeslices_pcolor2D,lats2D =np.meshgrid(-timeslices_pcolor,latsop) timeslices_scat2D,latsot2D =np.meshgrid(-timeslices_plot,latsot) ################################################### ### Aggregate data from model output ### (not necessary when *.pkl files already exist) ################################################### if os.path.isfile('./data_aggregated_history_p2_df.pkl')==True: data_aggregated_history_p2_df=pd.read_pickle('./data_aggregated_history_p2_df.pkl') data_aggregated_p2_arrays_df=pd.read_pickle('./data_aggregated_p2_arrays_df.pkl') data_aggregated_geography_df=pd.read_pickle('./data_aggregated_geography_df.pkl') data_aggregated_SeasonalityZonality_df=pd.read_pickle('./data_aggregated_SeasonalityZonality_df.pkl') else: data_aggregated_history_p2_df=pd.DataFrame([], columns=['run','age','pCO2','S0','years','year_last','ts_ann','ts_ann_mean','ts_ann_min','ts_ann_max','prc_ann',\ 'prc_ann_mean','snap_years','snap_radbal','snap_ts','status','ts_ann_500yr_trend'], index=paths) #, dtype=float[], columns=['dp13_ts_globmean', 'years'], index=[] paths_problematic_rev1=[] for path in paths[:]: print(path) idx_Ma=path.find('Ma') idx_ppm=path.find('ppm') age=int(path[idx_Ma-3:idx_Ma]) run=path[idx_Ma-3:] data_aggregated_history_p2_df.loc[path, ['run']]=[run] pCO2=int(path[idx_Ma+3:idx_ppm]) #np.loadtxt(path+'/pCO2.dat', dtype=float) S0=np.loadtxt(path+'/SolarConstant.dat', dtype=float) nc_f = path+'/history_p2.nc' if os.path.isfile(nc_f): nc_fid = Dataset(nc_f, 'r') dates=num2date(nc_fid.variables['Time'][:], units=nc_fid.variables['Time'].units, calendar=nc_fid.variables['Time'].calendar) years_strings=[time.strftime('%Y') for time in dates] years=np.array(list(map(int,years_strings))) year_last=years[-1] ferret_cmd='cancel data/all; use "'+nc_f+'"; let ts_ann_f=ts_ann[i=@ave,j=@ave]; ; let prc_ann_f=PRC_ANN[i=@ave,j=@ave]*360.;' pyferret.run(ferret_cmd) fldnames=['ts_ann','prc_ann'] flddata={} for ff in range(len(fldnames)): fld_dict = pyferret.getdata(fldnames[ff]+'_f', False) flddata[fldnames[ff]] = fld_dict['data'][0,0,0,:,0,0] ts_ann,prc_ann=[flddata.get(ff) for ff in fldnames] ts_ann_mean=np.mean(ts_ann[:]) prc_ann_mean=np.mean(prc_ann[:]) ts_ann_min=np.min(ts_ann[:]) ts_ann_max=np.max(ts_ann[:]) data_aggregated_history_p2_df.loc[path, ['years','year_last','ts_ann','ts_ann_mean','ts_ann_min','ts_ann_max','prc_ann','prc_ann_mean']]=\ [years,year_last,ts_ann,ts_ann_mean,ts_ann_min,ts_ann_max,prc_ann,prc_ann_mean] ts_ann_500yr=data_aggregated_history_p2_df.loc[path, ['ts_ann']][0][-500:] slope, intercept, r_value, p_value, std_err = scipy.stats.linregress(np.arange(500),ts_ann_500yr) data_aggregated_history_p2_df.loc[path,['ts_ann_500yr_trend']]=500.*slope if path in paths_problematic: status=' (not yet)' data_aggregated_history_p2_df.loc[path, ['run','age','pCO2','S0']]=[run,age, pCO2, S0] else: paths_problematic_rev1.append(path) data_aggregated_history_p2_df.to_pickle('./data_aggregated_history_p2_df.pkl') data_aggregated_p2_arrays_df=pd.DataFrame([], columns=['run','age','ts_ann_arr','ts_jja_arr','ts_djf_arr','prc_ann_arr','prc_jja_arr','prc_djf_arr',\ 'evap_ann_arr','evap_jja_arr','evap_djf_arr','slp_ann_arr','slp_jja_arr','slp_djf_arr',\ 'ts_seasonal_arr','ts_seasonal_arr_intp','ts_monthly_arr','prc_seasonal_arr','prc_seasonal_arr_intp','prc_monthly_arr',\ 'evap_seasonal_arr','evap_seasonal_arr_intp','evap_monthly_arr','slp_seasonal_arr_intp','slp_monthly_arr','slp_seasonal_arr',\ 'rbsr_ann','vegfrac_tree_globmean','vegfrac_grass_globmean','prc_monthly_arr_intp','ts_monthly_arr_intp','horo','frlnd',\ 'ts_ensmean_globmean','tsl_ensmean_globmean','tsl_ensmean_globmean_manual','merid_heattransport_atm','merid_heattransport_globint_atm',\ 'ts_zonal_ann_intp','ts_zonal_season_min_intp','ts_zonal_season_max_intp','ts_zonal_loc_min_intp','ts_zonal_loc_max_intp',\ 'ts_ann_trop','ts_ann_NorthP','ts_ann_SouthP','tsl_ann_trop','tsl_ann_NorthP','tsl_ann_SouthP','ts_gradient_slopes'], index=paths) #, dtype=float[], columns=['dp13_ts_globmean', 'years'], index=[] latsot_idx_NH=np.where((latsot>=30) & (latsot<=80))[0] latsot_idx_SH=np.where((latsot<=-30) & (latsot>=-80))[0] pred_xx_NH=latsot[latsot_idx_NH] pred_xx_SH=latsot[latsot_idx_SH] for path in paths[:]: print(path) idx_Ma=path.find('Ma') age=int(path[idx_Ma-3:idx_Ma]) run=path[idx_Ma-3:] data_aggregated_p2_arrays_df.loc[path, ['run']]=[run] if path not in paths_problematic: ferret_cmd='cancel data/all; use "'+path+'/history_p2.nc"; \ let ts_ann_arr_f=ts_ann[l=@ave]; let ts_jja_arr_f=ts_jja[l=@ave]; let ts_djf_arr_f=ts_djf[l=@ave]; \ let prc_ann_arr_f=prc_ann[l=@ave]; let prc_jja_arr_f=prc_jja[l=@ave]; let prc_djf_arr_f=prc_djf[l=@ave];\ let evap_ann_arr_f=e_ann[l=@ave]; let evap_jja_arr_f=e_jja[l=@ave]; let evap_djf_arr_f=e_djf[l=@ave]; \ let slp_ann_arr_f=slp_ann[l=@ave]; let slp_jja_arr_f=slp_jja[l=@ave]; let slp_djf_arr_f=slp_djf[l=@ave];\ let ts_ann_trop_f=ts_ann_arr_f[i=@ave,y=30S:30N@ave]; let ts_ann_NorthP_f=ts_ann_arr_f[i=@ave,y=60N:90N@ave]; let ts_ann_SouthP_f=ts_ann_arr_f[i=@ave,y=60S:90S@ave];' pyferret.run(ferret_cmd) fldnames=['ts_ann_arr','ts_jja_arr','ts_djf_arr',\ 'prc_ann_arr','prc_jja_arr','prc_djf_arr',\ 'evap_ann_arr','evap_jja_arr','evap_djf_arr',\ 'slp_ann_arr','slp_jja_arr','slp_djf_arr',\ 'ts_ann_trop','ts_ann_NorthP','ts_ann_SouthP'] fldopts=[0,0,0,0,0,0,0,0,0,0,0,0, 1,1,1] flddata={} for ff in range(len(fldnames)): fld_dict = pyferret.getdata(fldnames[ff]+'_f', False) if fldopts[ff]==0: flddata[ff] = np.ma.masked_values(np.rot90(fld_dict['data'][:,:,0,0,0,0]), -1e34); if fldopts[ff]==1: flddata[ff] = fld_dict['data'][0,0,0,0,0,0] ts_ann_arr,ts_jja_arr,ts_djf_arr,\ prc_ann_arr,prc_jja_arr,prc_djf_arr,\ evap_ann_arr,evap_jja_arr,evap_djf_arr,\ slp_ann_arr,slp_jja_arr,slp_djf_arr,\ ts_ann_trop,ts_ann_NorthP,ts_ann_SouthP=[flddata.get(ff) for ff in range(len(fldnames))] ts_zonal_ann=np.ma.average(ts_ann_arr,axis=1) ts_zonal_ann_intp=jan.interpol_atm_to_ocn_grid(ts_zonal_ann,axes='lat') slope_NH, intercept_NH, r_value, p_value, std_err = scipy.stats.linregress(pred_xx_NH,ts_zonal_ann_intp[latsot_idx_NH]) slope_SH, intercept_SH, r_value, p_value, std_err = scipy.stats.linregress(pred_xx_SH,ts_zonal_ann_intp[latsot_idx_SH]) data_aggregated_p2_arrays_df.loc[path, ['age','ts_ann_arr','ts_jja_arr','ts_djf_arr','prc_ann_arr','prc_jja_arr','prc_djf_arr',\ 'evap_ann_arr','evap_jja_arr','evap_djf_arr','slp_ann_arr','slp_jja_arr','slp_djf_arr']]=\ [age,ts_ann_arr,ts_jja_arr,ts_djf_arr,prc_ann_arr,prc_jja_arr,prc_djf_arr,\ evap_ann_arr,evap_jja_arr,evap_djf_arr,slp_ann_arr,slp_jja_arr,slp_djf_arr] #meso_v2_p2_0Ma_arrays_df.loc[path, ['age','ts_ann_arr','prc_ann_arr']]=[age,ts_ann_arr,prc_ann_arr] data_aggregated_p2_arrays_df.loc[path, ['ts_zonal_ann_intp','ts_ann_trop','ts_ann_NorthP','ts_ann_SouthP']]=\ [ts_zonal_ann_intp,ts_ann_trop,ts_ann_NorthP,ts_ann_SouthP] data_aggregated_p2_arrays_df.loc[path, ['ts_gradient_slopes']]=\ [np.array([slope_NH,slope_SH])] ferret_cmd='cancel data/all; use "'+path+'/ensmean_snapshots_potsdam2.nc"; \ let ts_monthly_f=ts; let prc_monthly_f=prc; let evap_monthly_f=e; let slp_monthly_f=slp; let rbsr_ann_f=rbsr[l=@ave] ;\ let ft_f=ft[l=1]; let fg_f=fg[l=1]; let frlnd_f=frlnd[l=1]; let horo_f=horo[l=1];\ let ts_ensmean_globmean_f=ts[i=@ave,j=@ave,l=@ave];\ let tsl_ensmean_globmean_f=tsl[i=@ave,j=@ave,l=@ave]; let tsl_ensmean_manual_f=ts[l=@ave]+gam[l=@ave]*horo[l=@ave]; let tsl_ensmean_globmean_manual_f=tsl_ensmean_manual_f[i=@ave,j=@ave];\ let tsl_ann_trop_f=tsl_ensmean_manual_f[i=@ave,y=30S:30N@ave]; let tsl_ann_NorthP_f=tsl_ensmean_manual_f[i=@ave,y=60N:90N@ave]; let tsl_ann_SouthP_f=tsl_ensmean_manual_f[i=@ave,y=60S:90S@ave];\ let mm=fta[l=@ave]+ftd[l=@ave]; let mm_abs=abs(mm); let merid_heattransport_atm_f=mm; let merid_heattransport_globint_atm_f=mm_abs[j=@ave];' pyferret.run(ferret_cmd) fldnames=['ts_monthly','prc_monthly','evap_monthly','slp_monthly',\ 'rbsr_ann','ft','fg','frlnd','horo','ts_ensmean_globmean','tsl_ensmean_globmean','tsl_ensmean_globmean_manual',\ 'merid_heattransport_atm','merid_heattransport_globint_atm',\ 'tsl_ann_trop','tsl_ann_NorthP','tsl_ann_SouthP'] fldopts=[0,0,0,0,\ 1,1,1,1,1,1,1,1,\ 2,2,\ 3,3,3] flddata={} for ff in range(len(fldnames)): fld_dict = pyferret.getdata(fldnames[ff]+'_f', False) if fldopts[ff]==0: flddata[ff] = np.ma.masked_values(np.rot90(fld_dict['data'][:,:,0,:,0,0]), -1e34); if fldopts[ff]==1: flddata[ff] = np.ma.masked_values(np.rot90(fld_dict['data'][:,:,0,0,0,0]), -1e34); if fldopts[ff]==2: flddata[ff] = np.ma.masked_values(fld_dict['data'][0,::-1,0,0,0,0], -1e34); if fldopts[ff]==3: flddata[ff] = fld_dict['data'][0,0,0,0,0,0]; ts_monthly,prc_monthly,evap_monthly,slp_monthly,\ rbsr_ann,ft,fg,frlnd,horo,ts_ensmean_globmean,tsl_ensmean_globmean,tsl_ensmean_globmean_manual,\ merid_heattransport_atm,merid_heattransport_globint_atm,\ tsl_ann_trop,tsl_ann_NorthP,tsl_ann_SouthP=[flddata.get(ff) for ff in range(len(fldnames))] vegfrac_tree_globmean=jan.global_mean(np.ma.masked_array(ft,mask=1.-1.*(frlnd>0.))) vegfrac_grass_globmean=jan.global_mean(np.ma.masked_array(fg,mask=1.-1.*(frlnd>0.))) prc_seasonal=np.stack((np.mean(np.stack((prc_monthly[:,:,11],prc_monthly[:,:,0],prc_monthly[:,:,1]),axis=2),axis=2),\ np.mean(prc_monthly[:,:,2:4+1],axis=2),np.mean(prc_monthly[:,:,5:7+1],axis=2),np.mean(prc_monthly[:,:,8:10+1],axis=2)),axis=2) evap_seasonal=np.stack((np.mean(np.stack((evap_monthly[:,:,11],evap_monthly[:,:,0],evap_monthly[:,:,1]),axis=2),axis=2),\ np.mean(evap_monthly[:,:,2:4+1],axis=2),np.mean(evap_monthly[:,:,5:7+1],axis=2),np.mean(evap_monthly[:,:,8:10+1],axis=2)),axis=2) ts_seasonal=np.stack((np.mean(np.stack((ts_monthly[:,:,11],ts_monthly[:,:,0],ts_monthly[:,:,1]),axis=2),axis=2),\ np.mean(ts_monthly[:,:,2:4+1],axis=2),np.mean(ts_monthly[:,:,5:7+1],axis=2),np.mean(ts_monthly[:,:,8:10+1],axis=2)),axis=2) slp_seasonal=np.stack((np.mean(np.stack((slp_monthly[:,:,11],slp_monthly[:,:,0],slp_monthly[:,:,1]),axis=2),axis=2),\ np.mean(slp_monthly[:,:,2:4+1],axis=2),np.mean(slp_monthly[:,:,5:7+1],axis=2),np.mean(slp_monthly[:,:,8:10+1],axis=2)),axis=2) ts_seasonal_intp=np.ma.empty((48,96,4)) prc_seasonal_intp=np.ma.empty((48,96,4)) evap_seasonal_intp=np.ma.empty((48,96,4)) slp_seasonal_intp=np.ma.empty((48,96,4)) for ss in range(4): ts_seasonal_intp[:,:,ss]=jan.interpol_atm_to_ocn_grid(ts_seasonal[:,:,ss]) prc_seasonal_intp[:,:,ss]=jan.interpol_atm_to_ocn_grid(prc_seasonal[:,:,ss]) evap_seasonal_intp[:,:,ss]=jan.interpol_atm_to_ocn_grid(evap_seasonal[:,:,ss]) slp_seasonal_intp[:,:,ss]=jan.interpol_atm_to_ocn_grid(slp_seasonal[:,:,ss]) prc_monthly_intp=np.ma.empty((48,96,12)) ts_monthly_intp=np.ma.empty((48,96,12)) for ss in range(12): ts_monthly_intp[:,:,ss]=jan.interpol_atm_to_ocn_grid(ts_monthly[:,:,ss]) prc_monthly_intp[:,:,ss]=jan.interpol_atm_to_ocn_grid(prc_monthly[:,:,ss]) ts_zonal_season_min=np.ma.average(np.ma.min(ts_seasonal,axis=2),axis=1) ts_zonal_season_max=np.ma.average(np.ma.max(ts_seasonal,axis=2),axis=1) ts_zonal_loc_min=np.min(np.min(ts_monthly,axis=2),axis=1) ts_zonal_loc_max=np.max(np.max(ts_monthly,axis=2),axis=1) ts_zonal_season_min_intp=jan.interpol_atm_to_ocn_grid(ts_zonal_season_min,axes='lat') ts_zonal_season_max_intp=jan.interpol_atm_to_ocn_grid(ts_zonal_season_max,axes='lat') ts_zonal_loc_min_intp=jan.interpol_atm_to_ocn_grid(ts_zonal_loc_min,axes='lat') ts_zonal_loc_max_intp=jan.interpol_atm_to_ocn_grid(ts_zonal_loc_max,axes='lat') data_aggregated_p2_arrays_df.loc[path, ['ts_monthly_arr','prc_monthly_arr','evap_monthly_arr','slp_monthly_arr','ts_monthly_arr_intp','prc_monthly_arr_intp']]=\ [ts_monthly,prc_monthly,evap_monthly,slp_monthly,ts_monthly_intp,prc_monthly_intp] data_aggregated_p2_arrays_df.loc[path, ['ts_seasonal_arr','prc_seasonal_arr','evap_seasonal_arr','slp_seasonal_arr']]=\ [ts_seasonal,prc_seasonal,evap_seasonal,slp_seasonal] data_aggregated_p2_arrays_df.loc[path, ['ts_seasonal_arr_intp','prc_seasonal_arr_intp','evap_seasonal_arr_intp','slp_seasonal_arr_intp','rbsr_ann']]=\ [ts_seasonal_intp,prc_seasonal_intp,evap_seasonal_intp,slp_seasonal_intp,rbsr_ann] data_aggregated_p2_arrays_df.loc[path, ['vegfrac_tree_globmean','vegfrac_grass_globmean']]=\ [vegfrac_tree_globmean,vegfrac_grass_globmean] data_aggregated_p2_arrays_df.loc[path, ['horo','frlnd']]=\ [horo,frlnd] data_aggregated_p2_arrays_df.loc[path, ['ts_ensmean_globmean','tsl_ensmean_globmean','tsl_ensmean_globmean_manual']]=\ [ts_ensmean_globmean,tsl_ensmean_globmean,tsl_ensmean_globmean_manual] data_aggregated_p2_arrays_df.loc[path, ['tsl_ann_trop','tsl_ann_NorthP','tsl_ann_SouthP']]=\ [tsl_ann_trop,tsl_ann_NorthP,tsl_ann_SouthP] data_aggregated_p2_arrays_df.loc[path, ['merid_heattransport_atm','merid_heattransport_globint_atm']]=\ [merid_heattransport_atm,merid_heattransport_globint_atm] #meso_v2_p2_0Ma_arrays_df.loc[path, ['ts_monthly_arr','prc_monthly_arr']]=[ts_monthly,prc_monthly] data_aggregated_p2_arrays_df.loc[path, ['ts_zonal_season_min_intp','ts_zonal_season_max_intp','ts_zonal_loc_min_intp','ts_zonal_loc_max_intp']]=\ [ts_zonal_season_min_intp,ts_zonal_season_max_intp,ts_zonal_loc_min_intp,ts_zonal_loc_max_intp] data_aggregated_p2_arrays_df.to_pickle('./data_aggregated_p2_arrays_df.pkl') data_aggregated_geography_df=pd.DataFrame([], columns=['run','age','continental_mask','shallowcoastal_mask','landfrac_glob','landarea_zonal','ocn_depth','horo_globmean'], index=paths) data_aggregated_geography_df.loc[:,['run','age']]=data_aggregated_history_p2_df.loc[:,['run','age']] for pp in paths[:]: #print(pp) if pp not in paths_problematic: [run,age]=data_aggregated_geography_df.loc[pp,['run','age']] nc_fid=Dataset(pp+'/topog.dta.nc') ht=1./100.*nc_fid.variables['ht'][::-1,:] #from cm to m kmt=nc_fid.variables['kmt'][::-1,:] data_aggregated_geography_df.loc[pp,['ocn_depth']]=[ht] #depth=np.loadtxt(pp+'/depth_adj.dat', dtype=int)[:,:] continental_mask=1.-1.*ht.mask #ht.np.zeros(ht.shape,dtype=int) #continental_mask[ht>0]=1 data_aggregated_geography_df.loc[pp,['continental_mask']]=[continental_mask] landarea_zonal=np.sum((1.-continental_mask)*gridcell_areas_mom,axis=1) data_aggregated_geography_df.loc[pp,['landarea_zonal']]=[landarea_zonal] dummymask=np.ones(continental_mask.shape) dummymask[kmt>12]=0 dummymask=ndimage.generic_filter(var_pad(dummymask.astype(int)), np.any, footprint=np.ones((3, 3))) dummymask=dummymask[1:-1,1:-1] dummymask[(dummymask==0) | (continental_mask==0)]=0 data_aggregated_geography_df.loc[pp,['shallowcoastal_mask']]=[1-dummymask] data_aggregated_geography_df.loc[pp,['landfrac_glob']]=jan.global_mean(np.ma.MaskedArray(1.-continental_mask)) #plt.imshow(shallowcoastal_mask,cmap='Reds') # nc_fid=Dataset(pp+'/topog.dta.nc') # ht=1./100.*nc_fid.variables['ht'][::-1,:] #from cm to m # data_aggregated_geography_df.loc[pp,['ocn_depth']]=[ht] nc_fid=Dataset(pp+'/ensmean_snapshots_potsdam2.nc') horo=nc_fid.variables['horo'][0,::-1,:] #from cm to m #horo=np.ma.MaskedArray(horo, mask=(horo==0)) horo_globmean=jan.global_mean(horo) data_aggregated_geography_df.loc[pp,['horo_globmean']]=[horo_globmean] data_aggregated_geography_df.to_pickle('./data_aggregated_geography_df.pkl') paths_indxs=['paths_1000ppm','paths_pCO2Foster','paths_pCO2COPSE','paths_250ppm','paths_500ppm','paths_2000ppm',\ 'paths_orb','paths_1000ppm_orb','paths_1000ppm_S0ini','paths_1000ppm_S0ini_VegFix','paths_SL','paths_1000ppm_VegHom','paths_VegHom_SL'] paths_list=[[paths_1000ppm],[paths_pCO2Foster],[paths_pCO2COPSE],[paths_250ppm],[paths_500ppm],[paths_2000ppm],\ [paths_orb],[paths_1000ppm_orb],[paths_1000ppm_S0ini],[paths_1000ppm_S0ini_VegFix],[paths_SL],[paths_1000ppm_VegHom],[paths_VegHom_SL]] data_aggregated_SeasonalityZonality_df=pd.DataFrame([],columns=['paths_list','prc_lonmon_maxs','prc_lonmon_maxs_lats','prc_lonmon_maxs_onlyChina','prc_lonmon_maxs_lats_onlyChina',\ 'monsoon_area_zonalsum','monsoon_area_zonalsum_onlyChina','mask_rectangle_onlyChina','monsoon_area_low_zonalsum','monsoon_area_low_zonalsum_onlyChina',\ 'prc_anns_zonalmean','prc_anns_globalmean','ts_anns_zonalmean','prc_variations_globalmean','ts_variations_globalmean','prc_season_variation_zonalmean',\ 'ts_season_variation_zonalmean','ts_anns_zonality','prc_anns_zonality','ts_zonality_globalmean','prc_zonality_globalmean',\ 'evap_anns_zonalmean','evap_variations_globalmean','evap_season_variation_zonalmean','evap_anns_zonality','evap_zonality_globalmean',\ 'aridity_continents_zonalmean','aridity_continents_globalmean','aridity_continents_zonality','aridity_continents_zonality_globalmean',\ 'rbsr_anns_zonalmean','rbsr_anns_globalmean','aridity_1p5_globarea','aridity_1p5_zonalarea','landareas_zonal',\ 'ts_anns_air_zonalmean','ts_grads',\ 'sst_anns_zonality','sst_zonality_globalmean','hmxl_anns_zonalmean',\ 'lats_land_ave',\ 'merid_heattransport_atms','merid_heattransport_ocnsurfs',\ 'ts_mon_variations_globalmean','ts_mon_variation_zonalmean'],\ index=paths_indxs) monsoon_maskOnlyChina_df=pd.DataFrame([],columns=['mask_rectangle_onlyChina'],index=timeslices) data_aggregated_SeasonalityZonality_df['paths_list']=paths_list for paths_ll in data_aggregated_SeasonalityZonality_df.index[:]: paths_pp=data_aggregated_SeasonalityZonality_df.loc[paths_ll,['paths_list']][0][0] prc_anns_zonalmean=np.ma.empty((48,len(paths_pp))) ts_anns_zonalmean=np.ma.empty((48,len(paths_pp))) evap_anns_zonalmean=np.ma.empty((48,len(paths_pp))) rbsr_anns_zonalmean=np.ma.empty((48,len(paths_pp))) sst_anns_zonalmean=np.ma.empty((48,len(paths_pp))) hmxl_anns_zonalmean=np.ma.empty((48,len(paths_pp))) prc_anns_globalmean=np.ma.empty(len(paths_pp)) rbsr_anns_globalmean=np.ma.empty(len(paths_pp)) prc_mons_variation_zonalmean=np.ma.empty((48,len(paths_pp))) ts_mons_variation_zonalmean=np.ma.empty((48,len(paths_pp))) evap_mons_variation_zonalmean=np.ma.empty((48,len(paths_pp))) prc_season_variation_zonalmean=np.ma.empty((48,len(paths_pp))) ts_season_variation_zonalmean=np.ma.empty((48,len(paths_pp))) ts_mon_variation_zonalmean=np.ma.empty((48,len(paths_pp))) evap_season_variation_zonalmean=np.ma.empty((48,len(paths_pp))) prc_variations_globalmean=np.ma.empty(len(paths_pp)) ts_variations_globalmean=np.ma.empty(len(paths_pp)) ts_mon_variations_globalmean=np.ma.empty(len(paths_pp)) evap_variations_globalmean=np.ma.empty(len(paths_pp)) prc_anns_zonality=np.ma.empty((48,len(paths_pp))) ts_anns_zonality=np.ma.empty((48,len(paths_pp))) evap_anns_zonality=np.ma.empty((48,len(paths_pp))) rbsr_anns_zonality=np.ma.empty((48,len(paths_pp))) aridity_continents_zonality=np.ma.empty((48,len(paths_pp))) sst_anns_zonality=np.ma.empty((48,len(paths_pp))) prc_zonality_globalmean=np.ma.empty(len(paths_pp)) prc_zonality_globalmean_test=np.ma.empty(len(paths_pp)) ts_zonality_globalmean=np.ma.empty(len(paths_pp)) evap_zonality_globalmean=np.ma.empty(len(paths_pp)) aridity_continents_zonality_globalmean=np.ma.empty(len(paths_pp)) sst_zonality_globalmean=np.ma.empty(len(paths_pp)) aridity_continents_zonalmean=np.ma.empty((48,len(paths_pp))) aridity_1p5_globarea=np.ma.empty(len(paths_pp)) aridity_1p5_zonalarea=np.ma.empty((48,len(paths_pp))) landareas_zonal=np.ma.empty((48,len(paths_pp))) ts_grads=np.ma.empty((2,len(paths_pp))) ts_anns_air_zonalmean=np.ma.empty((48,len(paths_pp))) lats_land_ave=np.ma.empty(len(paths_pp)) merid_heattransport_atms=np.ma.empty((48,len(paths_pp))) merid_heattransport_ocnsurfs=np.ma.empty((48,len(paths_pp))) for pp in range(len(paths_pp)): path=paths_pp[pp] if path not in paths_problematic: prc_ann_arr=data_aggregated_p2_arrays_df.loc[path, ['prc_ann_arr']].values[0] ts_ann_arr=data_aggregated_p2_arrays_df.loc[path, ['ts_ann_arr']].values[0] evap_ann_arr=prc_ann_arr-data_aggregated_p2_arrays_df.loc[path, ['evap_ann_arr']].values[0] # prc_mon_arr=data_aggregated_p2_arrays_df.loc[path, ['prc_monthly_arr']].values[0] ts_mon_arr=data_aggregated_p2_arrays_df.loc[path, ['ts_monthly_arr']].values[0] # evap_mon_arr=prc_mon_arr-data_aggregated_p2_arrays_df.loc[path, ['evap_monthly_arr']].values[0] prc_season_arr=data_aggregated_p2_arrays_df.loc[path, ['prc_seasonal_arr']].values[0] ts_season_arr=data_aggregated_p2_arrays_df.loc[path, ['ts_seasonal_arr']].values[0] evap_season_arr=prc_season_arr-data_aggregated_p2_arrays_df.loc[path, ['evap_seasonal_arr']].values[0] #sst_ann_arr=meso_rev1_mom_arrays_df.loc[path, ['sst_ann_arr']][0] #hmxl_ann_arr=meso_rev1_mom_arrays_df.loc[path, ['mixedlayerdepth_ann_arr']][0] rbsr_ann_arr=data_aggregated_p2_arrays_df.loc[path, ['rbsr_ann']][0] continental_mask=data_aggregated_geography_df.loc[path,['continental_mask']][0] landarea_per_lat=np.sum((1.-continental_mask)*gridcell_areas_mom,axis=1) landareas_zonal[:,pp]=landarea_per_lat merid_heattransport_atm=data_aggregated_p2_arrays_df.loc[path, ['merid_heattransport_atm']].values[0] #merid_heattransport_ocnsurf=meso_rev1_overturn_df.loc[path, ['merid_heattransport_ocnsurf']][0] ### calculate the average latitude of land area, se also Torsvik 2020 AGU Fig.9 latsot_2D_masked=(1.-continental_mask)*np.abs(latsot_2D) latsot_2D_masked=np.ma.MaskedArray(latsot_2D_masked,mask=continental_mask) lat_land_ave=np.ma.average(latsot_2D_masked,weights=gridcell_areas_mom) lats_land_ave[pp]=lat_land_ave # prc_mon_variation=(np.max(prc_mon_arr,axis=2)-np.min(prc_mon_arr,axis=2))/prc_ann_arr prc_season_variation=(np.max(prc_season_arr,axis=2)-np.min(prc_season_arr,axis=2))/prc_ann_arr # evap_mon_variation=np.max(evap_mon_arr,axis=2)-np.min(evap_mon_arr,axis=2)/evap_ann_arr evap_season_variation=(np.max(evap_season_arr,axis=2)-np.min(evap_season_arr,axis=2))#/evap_ann_arr #/evap_ann_arr ts_mon_variation=np.max(ts_mon_arr,axis=2)-np.min(ts_mon_arr,axis=2) ts_season_variation=np.max(ts_season_arr,axis=2)-np.min(ts_season_arr,axis=2) prc_ann_intp=jan.interpol_atm_to_ocn_grid(prc_ann_arr) ts_ann_intp=jan.interpol_atm_to_ocn_grid(ts_ann_arr) evap_ann_intp=jan.interpol_atm_to_ocn_grid(evap_ann_arr) rbsr_ann_intp=jan.interpol_atm_to_ocn_grid(rbsr_ann_arr) prc_season_variation_intp=jan.interpol_atm_to_ocn_grid(prc_season_variation) # evap_mon_variation_intp=jan.interpol_atm_to_ocn_grid(evap_mon_variation) evap_season_variation_intp=jan.interpol_atm_to_ocn_grid(evap_season_variation) ts_mon_variation_intp=jan.interpol_atm_to_ocn_grid(ts_mon_variation) ts_season_variation_intp=jan.interpol_atm_to_ocn_grid(ts_season_variation) merid_heattransport_atm_intp=jan.interpol_atm_to_ocn_grid(merid_heattransport_atm,axes='lat') prc_ann_intp_masked=np.ma.MaskedArray(prc_ann_intp,mask=continental_mask) ts_ann_intp_masked=np.ma.MaskedArray(ts_ann_intp,mask=continental_mask) evap_ann_intp_masked=np.ma.MaskedArray(evap_ann_intp,mask=continental_mask) rbsr_ann_intp_masked=np.ma.MaskedArray(rbsr_ann_intp,mask=continental_mask) aridity_continents=func_aridity_idx(prc_ann_intp_masked,rbsr_ann_intp_masked) aridity_1p5_zonalarea[:,pp]=np.ma.sum(1.*(aridity_continents>=1.5)*gridcell_areas_mom,axis=1) aridity_1p5_globarea[pp]=np.ma.sum(1.*(aridity_continents>=1.5)*gridcell_areas_mom)/np.ma.sum(1.*(1-continental_mask)*gridcell_areas_mom) # prc_mon_variation_intp_masked=np.ma.MaskedArray(prc_mon_variation_intp,mask=continental_mask) prc_season_variation_intp_masked=np.ma.MaskedArray(prc_season_variation_intp,mask=continental_mask) # evap_mon_variation_intp_masked=np.ma.MaskedArray(evap_mon_variation_intp,mask=continental_mask) evap_season_variation_intp_masked=np.ma.MaskedArray(evap_season_variation_intp,mask=continental_mask) ts_mon_variation_intp_masked=np.ma.MaskedArray(ts_mon_variation_intp,mask=continental_mask) ts_season_variation_intp_masked=np.ma.MaskedArray(ts_season_variation_intp,mask=continental_mask) prc_variations_globalmean[pp]=jan.global_mean(prc_season_variation_intp_masked) ts_variations_globalmean[pp]=jan.global_mean(ts_season_variation_intp_masked) ts_mon_variations_globalmean[pp]=jan.global_mean(ts_mon_variation_intp_masked) evap_variations_globalmean[pp]=jan.global_mean(evap_season_variation_intp_masked) prc_anns_zonalmean[:,pp]=np.ma.mean(prc_ann_intp_masked,axis=1) ts_anns_zonalmean[:,pp]=np.ma.mean(ts_ann_intp_masked,axis=1) ts_anns_air_zonalmean[:,pp]=np.ma.mean(ts_ann_intp,axis=1) evap_anns_zonalmean[:,pp]=np.ma.mean(evap_ann_intp_masked,axis=1) rbsr_anns_zonalmean[:,pp]=np.ma.mean(rbsr_ann_intp_masked,axis=1) #aridity_continents_zonalmean[:,pp]=np.ma.mean(rbsr_ann_intp_masked,axis=1) aridity_continents_zonalmean=func_aridity_idx(prc_anns_zonalmean,rbsr_anns_zonalmean) ts_ann_equator=np.ma.average(ts_ann_intp[np.where(latsot<=30)[0][0]:np.where(latsot>=-30)[0][-1]+1,:],weights=gridcell_areas_mom[np.where(latsot<=30)[0][0]:np.where(latsot>=-30)[0][-1]+1,:]) ts_ann_NP=np.ma.average(ts_ann_intp[:np.where(latsot>=60)[0][-1]+1,:],weights=gridcell_areas_mom[:np.where(latsot>=60)[0][-1]+1,:]) ts_ann_SP=np.ma.average(ts_ann_intp[np.where(latsot<=-60)[0][0]:,:],weights=gridcell_areas_mom[np.where(latsot<=-60)[0][0]:,:]) ts_grads[:,pp]=[ts_ann_equator-ts_ann_NP,ts_ann_equator-ts_ann_SP] # prc_mons_variation_zonalmean[:,pp]=np.ma.mean(prc_mon_variation_intp_masked,axis=1) prc_season_variation_zonalmean[:,pp]=np.ma.mean(prc_season_variation_intp_masked,axis=1) # evap_mons_variation_zonalmean[:,pp]=np.ma.mean(evap_mon_variation_intp_masked,axis=1) evap_season_variation_zonalmean[:,pp]=np.ma.mean(evap_season_variation_intp_masked,axis=1) ts_mon_variation_zonalmean[:,pp]=np.ma.mean(ts_mon_variation_intp_masked,axis=1) ts_season_variation_zonalmean[:,pp]=np.ma.mean(ts_season_variation_intp_masked,axis=1) ts_anns_zonality[:,pp]=np.ma.mean(np.ma.abs(np.transpose(np.tile(ts_anns_zonalmean[:,pp],(96,1)))-ts_ann_intp_masked),axis=1) #ts_anns_zonality[:,pp]=np.ma.max(ts_ann_intp_masked,axis=1)-np.ma.min(ts_ann_intp_masked,axis=1) prc_anns_zonality[:,pp]=np.ma.mean(np.ma.abs(np.transpose(np.tile(prc_anns_zonalmean[:,pp],(96,1)))-prc_ann_intp_masked),axis=1)/prc_anns_zonalmean[:,pp] #prc_anns_zonality[:,pp]=(np.ma.max(prc_ann_intp_masked,axis=1)-np.ma.min(prc_ann_intp_masked,axis=1))/prc_anns_zonalmean[:,pp] evap_anns_zonality[:,pp]=np.ma.max(evap_ann_intp_masked,axis=1)-np.ma.min(evap_ann_intp_masked,axis=1)#/np.ma.mean(evap_ann_intp_masked,axis=1) #(np.ma.max(evap_ann_intp_masked,axis=1)-np.ma.min(evap_ann_intp_masked,axis=1))/np.ma.mean(evap_ann_intp_masked,axis=1) rbsr_anns_zonality[:,pp]=np.ma.mean(np.ma.abs(np.transpose(np.tile(rbsr_anns_zonalmean[:,pp],(96,1)))-rbsr_ann_intp_masked),axis=1) aridity_continents_zonality[:,pp]=func_aridity_idx(prc_anns_zonality[:,pp],rbsr_anns_zonality[:,pp]) aridity_continents_zonality[:,pp]=np.ma.mean(np.ma.abs(np.transpose(np.tile(aridity_continents_zonalmean[:,pp],(96,1)))-aridity_continents),axis=1) #/np.ma.mean(aridity_continents,axis=1) #sst_anns_zonality[:,pp]=np.ma.mean(np.ma.abs(np.transpose(np.tile(sst_anns_zonalmean[:,pp],(96,1)))-sst_ann_arr),axis=1) # Should be weighted by area per latitude! Or not? prc_zonality_globalmean[pp]=np.ma.average(prc_anns_zonality[:,pp],weights=np.cos(pi_frac)) evap_zonality_globalmean[pp]=np.ma.average(evap_anns_zonality[:,pp],weights=np.cos(pi_frac))#,weights=landarea_per_lat) #ts_zonality_globalmean[pp]=np.ma.average(ts_anns_zonality[:,pp],weights=np.cos(pi_frac)) #,weights=landarea_per_lat) #prc_zonality_globalmean_test[pp]=np.ma.average(prc_anns_zonality[:,pp],weights=np.cos(pi_frac)) ts_zonality_globalmean[pp]=jan.global_mean(np.ma.abs(np.transpose(np.tile(ts_anns_zonalmean[:,pp],(96,1)))-ts_ann_intp_masked)) #sst_zonality_globalmean[pp]=jan.global_mean(np.ma.abs(np.transpose(np.tile(sst_anns_zonalmean[:,pp],(96,1)))-sst_ann_arr)) #rbsr_anns_globalmean[pp]=np.ma.average(rbsr_anns_zonalmean[:,pp],weights=landarea_per_lat) rbsr_anns_globalmean[pp]=jan.global_mean(rbsr_ann_intp_masked) prc_anns_globalmean[pp]=jan.global_mean(prc_ann_intp_masked) aridity_continents_zonality_globalmean[pp]=np.ma.average(aridity_continents_zonality[:,pp],weights=landarea_per_lat) merid_heattransport_atms[:,pp]=merid_heattransport_atm_intp #merid_heattransport_ocnsurfs[:,pp]=np.ma.MaskedArray(merid_heattransport_ocnsurf,mask=np.all(1.-continental_mask,axis=1)) aridity_continents_globalmean=func_aridity_idx(prc_anns_globalmean,rbsr_anns_globalmean) data_aggregated_SeasonalityZonality_df.loc[paths_ll,['prc_anns_zonalmean','ts_anns_zonalmean','prc_variations_globalmean','ts_variations_globalmean','prc_season_variation_zonalmean',\ 'ts_season_variation_zonalmean','ts_anns_zonality','prc_anns_zonality','ts_zonality_globalmean','prc_zonality_globalmean',\ 'evap_anns_zonalmean','evap_variations_globalmean','evap_season_variation_zonalmean','evap_anns_zonality','evap_zonality_globalmean',\ 'aridity_continents_zonalmean','aridity_continents_globalmean','aridity_continents_zonality','aridity_continents_zonality_globalmean']]\ =[prc_anns_zonalmean,ts_anns_zonalmean,prc_variations_globalmean,ts_variations_globalmean,prc_season_variation_zonalmean,ts_season_variation_zonalmean,\ ts_anns_zonality,prc_anns_zonality,ts_zonality_globalmean,prc_zonality_globalmean,\ evap_anns_zonalmean,evap_variations_globalmean,evap_season_variation_zonalmean,evap_anns_zonality,evap_zonality_globalmean,\ aridity_continents_zonalmean,aridity_continents_globalmean,aridity_continents_zonality,aridity_continents_zonality_globalmean] data_aggregated_SeasonalityZonality_df.loc[paths_ll,['ts_mon_variations_globalmean','ts_mon_variation_zonalmean']]=[ts_mon_variations_globalmean,ts_mon_variation_zonalmean] data_aggregated_SeasonalityZonality_df.loc[paths_ll,['rbsr_anns_zonalmean','rbsr_anns_globalmean','prc_anns_globalmean','aridity_1p5_globarea','aridity_1p5_zonalarea','landareas_zonal']]=\ [rbsr_anns_zonalmean,rbsr_anns_globalmean,prc_anns_globalmean,aridity_1p5_globarea,aridity_1p5_zonalarea,landareas_zonal] data_aggregated_SeasonalityZonality_df.loc[paths_ll,['ts_anns_air_zonalmean','ts_grads']]=\ [ts_anns_air_zonalmean,ts_grads] data_aggregated_SeasonalityZonality_df.loc[paths_ll,['lats_land_ave']]=[lats_land_ave] #data_aggregated_SeasonalityZonality_df.loc[paths_ll,['merid_heattransport_atms','merid_heattransport_ocnsurfs']]=[merid_heattransport_atms,merid_heattransport_ocnsurfs] data_aggregated_SeasonalityZonality_df.loc[paths_ll,['merid_heattransport_atms']]=[merid_heattransport_atms] data_aggregated_SeasonalityZonality_df.to_pickle('./data_aggregated_SeasonalityZonality_df.pkl')