<jupyter_start><jupyter_text>budget_tables_draft2017b.ipynb PurposeProduce tables of MARC aerosol budget for *"Effective radiative forcing in the aerosol-climate model CAM5.3-MARC-ARG compared to default CAM5.3"*, based on results from project [p17c-marc-comparison](https://github.com/grandey/p17c-marc-comparison). Requirements- NetCDF files produced by [budget_cdo_nco_draft2017b.ipynb](https://github.com/grandey/p17c-marc-comparison/blob/master/analysis_draft2017b/budget_cdo_nco_draft2017b.ipynb) and [analysis_cdo_nco_draft2017b.ipynb](https://github.com/grandey/p17c-marc-comparison/blob/master/analysis_draft2017b/analysis_cdo_nco_draft2017b.ipynb).- Python modules mentioned in cell 2, inlcuding [climapy](https://github.com/grandey/climapy) ([doi:10.5281/zenodo.1053020](https://doi.org/10.5281/zenodo.1053020)). Author<NAME>, 2018 AcknowledgmentsThanks are due to **<NAME>**.  The interpretation of the diagnostics are informed by *print_MARC_budget_Benjamin.pro*, an IDL script written by <NAME>.<jupyter_code>! date
import climapy
import os
import numpy as np
import pandas as pd
import scipy
from scipy import stats
import xarray as xr

# Print versions of packages
for module in [climapy, np, pd, scipy, xr]:
    try:
        print('{}.__version__ = {}'.format(module.__name__, module.__version__))
    except AttributeError:
        pass
    
# Default precision when displaying DataFrames
pd.set_option('precision', 3)<jupyter_output>climapy.__version__ = 0.1.0+b35937c
numpy.__version__ = 1.14.3
pandas.__version__ = 0.23.0
scipy.__version__ = 1.1.0
xarray.__version__ = 0.10.4<jupyter_text>Post-processed data file locationsLocation of files written by [budget_cdo_nco_draft2017b.ipynb](https://github.com/grandey/p17c-marc-comparison/blob/master/analysis_draft2017b/budget_cdo_nco_draft2017b.ipynb) and [analysis_cdo_nco_draft2017b.ipynb](https://github.com/grandey/p17c-marc-comparison/blob/master/analysis_draft2017b/analysis_cdo_nco_draft2017b.ipynb). The latter contains the column burdens.<jupyter_code>budget_dir = os.path.expandvars('$HOME/data/projects/p17c_marc_comparison/budget_cdo_nco_draft2017b/')
analysis_dir = os.path.expandvars('$HOME/data/projects/p17c_marc_comparison/analysis_cdo_nco_draft2017b/')<jupyter_output><empty_output><jupyter_text>Read global totals into dictionary<jupyter_code># Diagnostics to load data for
diagnostic_dict = {
    # Sources
    'c{}_SOURCE': 'Sources, Tg/yr',  # sum of sources
    'c{}_SF': '    Emission',
    'c{}_BNUC': '    Binary nucleation',
    'c{}_CND': '    Condensation',
    'c{}_AGEsrc': '    Aging (source)',
    'c{}_ADJsrc': '    Growth (source)',
    'c{}_CGsrc': '    Coagulation (source)',
    'c{}_EVAP': '    Hydrometeor evaporation',
    # Sinks
    'c{}_SINK': 'Sinks, Tg/yr',  # sum of sinks
    'c{}_AGEsnk': '    Aging (sink)',
    'c{}_ADJsnk': '    Growth (sink)',      
    'c{}_CGsnk': '    Coagulation (sink)',             
    'c{}_ACT': '    Nucleation scavenging by stratiform clouds',
    'c{}_SCV': '    Nucleation scavenging by convective clouds',
    'c{}_IMP': '    Impaction scavenging',
    'c{}_DRY': '    Dry deposition',
    # Burden / loading
    'c{}_LDG': 'Burden, Tg'
}
# List of variables (aerosol-diagnostic combinations) to try loading data for
variable_list = []
for aerosol in ['tSUL', 'tOC', 'tBC', 'tDST', 'tSSLT',
                'NUC', 'AIT', 'ACC', 'OC', 'MOS', 'OIM', 'SIMOS', 'BC', 'MBS', 'BIM', 'SIMBS']:
    for diagnostic in diagnostic_dict.keys():
        variable_list.append(diagnostic.format(aerosol)) 
# Read data and calculate global totals, standard errors
# Initialise dictionaries
mean_dict = {}  # mean global totals
error_dict = {}  # standard errors, based on annual means
# Loop over years and variables
for year in ['2000', '1850']:
    for variable in variable_list:
        variable_year = (variable, year)
        try:
            # Read data from input file
            if '_LDG' in variable:
                in_filename = '{}/marc_s2_{}_{}_ANN.nc'.format(analysis_dir, year, variable)
                in_ds = xr.open_dataset(in_filename, decode_times=False)[variable]
            else:
                in_filename = '{}/marc_s2_{}_{}_ANN.nc'.format(budget_dir, year, variable)
                in_ds = xr.open_dataset(in_filename, decode_times=False)[variable]
            # Exclude first two years as spin-up
            data = in_ds.isel(time=slice(2, None))
            # Calculate sum across earth: kg/m2/s -> kg/s / kg/m2 -> kg
            totals = climapy.xr_area_weighted_stat(data, stat='sum')
            # Convert units
            if '_LDG' in variable:  # kg -> Tg
                totals = totals / 1e9
            else:  # kg/s -> Tg/yr
                totals = totals / 1e9 * 60 * 60 *24 * 365
            # Calculate mean across years
            mean_dict[variable_year] = float(totals.mean())
            # Calculate standard error
            error = totals.std(dim='time', ddof=1).data.flatten() / np.sqrt(totals['time'].size)
            error_dict[variable_year] = float(error)
            # Print mean
            print('{:13}, {}: {:10.4f} ± {:7.4f}'.format(variable, year, mean_dict[variable_year],
                                                        error_dict[variable_year]))
            # Close input file
            in_ds.close()   
        except FileNotFoundError:
            pass
! date<jupyter_output>ctSUL_SOURCE , 2000:   538.7005 ±  0.7464
ctSUL_BNUC   , 2000:     0.0019 ±  0.0001
ctSUL_CND    , 2000:    11.5506 ±  0.0531
ctSUL_AGEsrc , 2000:    22.3322 ±  0.0203
ctSUL_EVAP   , 2000:   504.8158 ±  0.7477
ctSUL_SINK   , 2000:  -538.8804 ±  0.7439
ctSUL_ACT    , 2000:  -396.5978 ±  0.6930
ctSUL_SCV    , 2000:   -20.7465 ±  0.0611
ctSUL_IMP    , 2000:  -116.2898 ±  0.0614
ctSUL_DRY    , 2000:    -5.2464 ±  0.0102
ctSUL_LDG    , 2000:     1.3313 ±  0.0028
ctOC_SOURCE  , 2000:   108.8810 ±  0.0000
ctOC_SF      , 2000:   108.8810 ±  0.0000
ctOC_SINK    , 2000:  -108.8222 ±  0.0086
ctOC_ACT     , 2000:   -13.1431 ±  0.0146
ctOC_SCV     , 2000:    -0.5127 ±  0.0021
ctOC_IMP     , 2000:   -90.6522 ±  0.0175
ctOC_DRY     , 2000:    -4.5141 ±  0.0047
ctOC_LDG     , 2000:     1.4934 ±  0.0028
ctBC_SOURCE  , 2000:     7.7587 ±  0.0000
ctBC_SF      , 2000:     7.7587 ±  0.0000
ctBC_SINK    , 2000:    -7.7653 ±  0.0011
ctBC_ACT     , 2000:    -1.0153 ±  0.0019
ctBC_SCV     , 2000:    -0.0277 ± [...]<jupyter_text>Note: the values are ~0.1% lower than Alex's values, due to using slightly different values for the area of the earth.  Rearrange into DataFrames<jupyter_code># Initialise dataframes
columns = ['tSUL', 'tOC', 'tBC', 'tDST', 'tSSLT',
           'NUC', 'AIT', 'ACC', 'OC', 'MOS', 'OIM', 'SIMOS', 'BC', 'MBS', 'BIM', 'SIMBS']
mean_2000_df = pd.DataFrame(columns=columns)  # marc_s2_2000
error_2000_df = pd.DataFrame(columns=columns)
mean_1850_df = pd.DataFrame(columns=columns)  # marc_s2_1850
error_1850_df = pd.DataFrame(columns=columns)
# Populate with data for diagnostics
for k, v in diagnostic_dict.items():  # loop over diagnostics
    for col in columns:  # loop over columns
        if (k.format(col), '2000') in mean_dict:
            mean_2000_df.loc[v, col] = mean_dict[(k.format(col), '2000')]
            error_2000_df.loc[v, col] = error_dict[(k.format(col), '2000')]
            mean_1850_df.loc[v, col] = mean_dict[(k.format(col), '1850')]
            error_1850_df.loc[v, col] = error_dict[(k.format(col), '1850')]
# Lifetime
mean_2000_df.loc['Lifetime, days'] = (mean_2000_df.loc['Burden, Tg'] /
                                      mean_2000_df.loc['Sinks, Tg/yr'].abs() * 365)
error_2000_df.loc['Lifetime, days'] = (mean_2000_df.loc['Lifetime, days'] *
                                       ((error_2000_df.loc['Burden, Tg'] /
                                         mean_2000_df.loc['Burden, Tg'])**2 +
                                        (error_2000_df.loc['Sinks, Tg/yr'] /
                                         mean_2000_df.loc['Sinks, Tg/yr'].abs())**2).pow(0.5))
mean_1850_df.loc['Lifetime, days'] = (mean_1850_df.loc['Burden, Tg'] /
                                      mean_1850_df.loc['Sinks, Tg/yr'] * -365)
error_1850_df.loc['Lifetime, days'] = (mean_1850_df.loc['Lifetime, days'] *
                                       ((error_1850_df.loc['Burden, Tg'] /
                                         mean_1850_df.loc['Burden, Tg'])**2 +
                                        (error_1850_df.loc['Sinks, Tg/yr'] /
                                         mean_1850_df.loc['Sinks, Tg/yr'].abs())**2).pow(0.5))
# Print means for marc_s2_2000
mean_2000_df
# Print means for marc_s2_1850
mean_1850_df<jupyter_output><empty_output><jupyter_text>Sulfate tables<jupyter_code># Year-2000
# Sulfate columns
sul_cols = ['tSUL', 'NUC', 'AIT', 'ACC', 'SIMOS', 'SIMBS',]
# Extract sulfate data
sul_mean_2000_df = mean_2000_df[sul_cols]
sul_error_2000_df = error_2000_df[sul_cols]
# Drop rows with no data
sul_mean_2000_df = sul_mean_2000_df.dropna(how='all')
sul_error_2000_df = sul_error_2000_df.dropna(how='all')
# Convert values to strings
sul_mean_2000_df = sul_mean_2000_df.applymap(lambda x: '%+.2f' % x)
sul_error_2000_df = sul_error_2000_df.applymap(lambda x: '%.2f' % x)
# Combine mean and standard error strings
sul_2000_df = sul_mean_2000_df + ' ± ' + sul_error_2000_df
# Replace 'nan ± nan' with ''
sul_2000_df = sul_2000_df.replace('nan ± nan', '')
sul_2000_df = sul_2000_df.replace('+nan ± nan', '')
# Rename columns
sul_2000_df = sul_2000_df.rename(index=str, columns={'tSUL': 'Sulfate aerosol',
                                                     'NUC': 'in NUC',
                                                     'AIT': 'in AIT',
                                                     'ACC': 'in ACC',
                                                     'SIMOS': 'in MOS',
                                                     'SIMBS': 'in MBS'})
# Replace Tg with Tg(SO4)
for i in sul_2000_df.index:  # not an efficient method, but it works!
    if 'Tg' in i:
        sul_2000_df = sul_2000_df.rename({i: 'Tg(SO4)'.join(i.split('Tg'))})
# Save to CSV
out_filename = 'csv/t01.csv'
sul_2000_df.to_csv(out_filename)
print('Written {}'.format(out_filename))
# Save to Excel
out_filename = 'csv/t01.xlsx'
sul_2000_df.to_excel(out_filename)
print('Written {}'.format(out_filename))
# Display
sul_2000_df.style.set_caption('Year-2000')
# Year-1850
# Extract sulfate data
sul_mean_1850_df = mean_1850_df[sul_cols]
sul_error_1850_df = error_1850_df[sul_cols]
# Drop rows with no data
sul_mean_1850_df = sul_mean_1850_df.dropna(how='all')
sul_error_1850_df = sul_error_1850_df.dropna(how='all')
# Convert values to strings
sul_mean_1850_df = sul_mean_1850_df.applymap(lambda x: '%+.2f' % x)
sul_error_1850_df = sul_error_1850_df.applymap(lambda x: '%.2f' % x)
# Combine mean and standard error strings
sul_1850_df = sul_mean_1850_df + ' ± ' + sul_error_1850_df
# Replace 'nan ± nan' with ''
sul_1850_df = sul_1850_df.replace('nan ± nan', '')
sul_1850_df = sul_1850_df.replace('+nan ± nan', '')
# Rename columns
sul_1850_df = sul_1850_df.rename(index=str, columns={'tSUL': 'Sulfate aerosol',
                                                     'NUC': 'in NUC',
                                                     'AIT': 'in AIT',
                                                     'ACC': 'in ACC',
                                                     'SIMOS': 'in MOS',
                                                     'SIMBS': 'in MBS'})
# Replace Tg with Tg(SO4)
for i in sul_1850_df.index:
    if 'Tg' in i:
        sul_1850_df = sul_1850_df.rename({i: 'Tg(SO4)'.join(i.split('Tg'))})
# Save to CSV
out_filename = 'csv/tS01.csv'
sul_1850_df.to_csv(out_filename)
print('Written {}'.format(out_filename))
# Save to Excel
out_filename = 'csv/tS01.xlsx'
sul_1850_df.to_excel(out_filename)
print('Written {}'.format(out_filename))
# Display
sul_1850_df.style.set_caption('Year-1850')<jupyter_output>Written csv/tS01.csv
Written csv/tS01.xlsx<jupyter_text>Organic carbon (OC) tables<jupyter_code># Year-2000
# OC columns
oc_cols = ['tOC', 'OC', 'OIM',]
# Extract OC data
oc_mean_2000_df = mean_2000_df[oc_cols]
oc_error_2000_df = error_2000_df[oc_cols]
# Drop rows with no data
oc_mean_2000_df = oc_mean_2000_df.dropna(how='all')
oc_error_2000_df = oc_error_2000_df.dropna(how='all')
# Convert values to strings
oc_mean_2000_df = oc_mean_2000_df.applymap(lambda x: '%+.2f' % x)
oc_error_2000_df = oc_error_2000_df.applymap(lambda x: '%.2f' % x)
# Combine mean and standard error strings
oc_2000_df = oc_mean_2000_df + ' ± ' + oc_error_2000_df
# Replace 'nan ± nan' with ''
oc_2000_df = oc_2000_df.replace('nan ± nan', '')
oc_2000_df = oc_2000_df.replace('+nan ± nan', '')
# Rename columns
oc_2000_df = oc_2000_df.rename(index=str, columns={'tOC': 'Organic carbon aerosol',
                                                   'OC': 'in pure OC',
                                                   'OIM': 'in MOS'})
# Save to CSV
out_filename = 'csv/t02.csv'
oc_2000_df.to_csv(out_filename)
print('Written {}'.format(out_filename))
# Save to Excel
out_filename = 'csv/t02.xlsx'
oc_2000_df.to_excel(out_filename)
print('Written {}'.format(out_filename))
# Display
oc_2000_df.style.set_caption('Year-2000')
# Year-1850
# OC columns
oc_cols = ['tOC', 'OC', 'OIM',]
# Extract OC data
oc_mean_1850_df = mean_1850_df[oc_cols]
oc_error_1850_df = error_1850_df[oc_cols]
# Drop rows with no data
oc_mean_1850_df = oc_mean_1850_df.dropna(how='all')
oc_error_1850_df = oc_error_1850_df.dropna(how='all')
# Convert values to strings
oc_mean_1850_df = oc_mean_1850_df.applymap(lambda x: '%+.2f' % x)
oc_error_1850_df = oc_error_1850_df.applymap(lambda x: '%.2f' % x)
# Combine mean and standard error strings
oc_1850_df = oc_mean_1850_df + ' ± ' + oc_error_1850_df
# Replace 'nan ± nan' with ''
oc_1850_df = oc_1850_df.replace('nan ± nan', '')
oc_1850_df = oc_1850_df.replace('+nan ± nan', '')
# Rename columns
oc_1850_df = oc_1850_df.rename(index=str, columns={'tOC': 'Organic carbon aerosol',
                                                   'OC': 'in pure OC',
                                                   'OIM': 'in MOS'})
# Save to CSV
out_filename = 'csv/tS02.csv'
oc_1850_df.to_csv(out_filename)
print('Written {}'.format(out_filename))
# Save to Excel
out_filename = 'csv/tS02.xlsx'
oc_1850_df.to_excel(out_filename)
print('Written {}'.format(out_filename))
# Display
oc_1850_df.style.set_caption('Year-1850')<jupyter_output>Written csv/tS02.csv
Written csv/tS02.xlsx<jupyter_text>Black carbon (BC) tables<jupyter_code># Year-2000
# BC columns
bc_cols = ['tBC', 'BC', 'BIM',]
# Extract BC data
bc_mean_2000_df = mean_2000_df[bc_cols]
bc_error_2000_df = error_2000_df[bc_cols]
# Drop rows with no data
bc_mean_2000_df = bc_mean_2000_df.dropna(how='all')
bc_error_2000_df = bc_error_2000_df.dropna(how='all')
# Convert values to strings
bc_mean_2000_df = bc_mean_2000_df.applymap(lambda x: '%+.2f' % x)
bc_error_2000_df = bc_error_2000_df.applymap(lambda x: '%.2f' % x)
# Combine mean and standard error strings
bc_2000_df = bc_mean_2000_df + ' ± ' + bc_error_2000_df
# Replace 'nan ± nan' with ''
bc_2000_df = bc_2000_df.replace('nan ± nan', '')
bc_2000_df = bc_2000_df.replace('+nan ± nan', '')
# Rename columns
bc_2000_df = bc_2000_df.rename(index=str, columns={'tBC': 'Black carbon aerosol',
                                                   'BC': 'in pure BC',
                                                   'BIM': 'in MBS'})
# Save to CSV
out_filename = 'csv/t03.csv'
bc_2000_df.to_csv(out_filename)
print('Written {}'.format(out_filename))
# Save to Excel
out_filename = 'csv/t03.xlsx'
bc_2000_df.to_excel(out_filename)
print('Written {}'.format(out_filename))
# Display
bc_2000_df.style.set_caption('Year-2000')
# Year-1850
# BC columns
bc_cols = ['tBC', 'BC', 'BIM',]
# Extract BC data
bc_mean_1850_df = mean_1850_df[bc_cols]
bc_error_1850_df = error_1850_df[bc_cols]
# Drop rows with no data
bc_mean_1850_df = bc_mean_1850_df.dropna(how='all')
bc_error_1850_df = bc_error_1850_df.dropna(how='all')
# Convert values to strings
bc_mean_1850_df = bc_mean_1850_df.applymap(lambda x: '%+.2f' % x)
bc_error_1850_df = bc_error_1850_df.applymap(lambda x: '%.2f' % x)
# Combine mean and standard error strings
bc_1850_df = bc_mean_1850_df + ' ± ' + bc_error_1850_df
# Replace 'nan ± nan' with ''
bc_1850_df = bc_1850_df.replace('nan ± nan', '')
bc_1850_df = bc_1850_df.replace('+nan ± nan', '')
# Rename columns
bc_1850_df = bc_1850_df.rename(index=str, columns={'tBC': 'Black carbon aerosol',
                                                   'BC': 'in pure BC',
                                                   'BIM': 'in MBS'})
# Save to CSV
out_filename = 'csv/tS03.csv'
bc_1850_df.to_csv(out_filename)
print('Written {}'.format(out_filename))
# Save to Excel
out_filename = 'csv/tS03.xlsx'
bc_1850_df.to_excel(out_filename)
print('Written {}'.format(out_filename))
# Display
bc_1850_df.style.set_caption('Year-1850')<jupyter_output>Written csv/tS03.csv
Written csv/tS03.xlsx<jupyter_text>Dust and sea-salt<jupyter_code># Year-2000
# Columns
dust_salt_cols = ['tSSLT', 'tDST']
# Extract data
dust_salt_mean_2000_df = mean_2000_df[dust_salt_cols]
dust_salt_error_2000_df = error_2000_df[dust_salt_cols]
# Drop rows with no data
dust_salt_mean_2000_df = dust_salt_mean_2000_df.dropna(how='all')
dust_salt_error_2000_df = dust_salt_error_2000_df.dropna(how='all')
# Convert values to strings
dust_salt_mean_2000_df = dust_salt_mean_2000_df.applymap(lambda x: '%+.2f' % x)
dust_salt_error_2000_df = dust_salt_error_2000_df.applymap(lambda x: '%.2f' % x)
# Combine mean and standard error strings
dust_salt_2000_df = dust_salt_mean_2000_df + ' ± ' + dust_salt_error_2000_df
# Replace 'nan ± nan' with ''
dust_salt_2000_df = dust_salt_2000_df.replace('nan ± nan', '')
dust_salt_2000_df = dust_salt_2000_df.replace('+nan ± nan', '')
# Rename columns
dust_salt_2000_df = dust_salt_2000_df.rename(index=str, columns={'tSSLT': 'Sea salt aerosol',
                                                                 'tDST': 'Dust aerosol',})
# Save to CSV
out_filename = 'csv/t04.csv'
dust_salt_2000_df.to_csv(out_filename)
print('Written {}'.format(out_filename))
# Save to Excel
out_filename = 'csv/t04.xlsx'
dust_salt_2000_df.to_excel(out_filename)
print('Written {}'.format(out_filename))
# Display
dust_salt_2000_df.style.set_caption('Year-2000')
# Year-1850
# Columns
dust_salt_cols = ['tSSLT', 'tDST']
# Extract data
dust_salt_mean_1850_df = mean_1850_df[dust_salt_cols]
dust_salt_error_1850_df = error_1850_df[dust_salt_cols]
# Drop rows with no data
dust_salt_mean_1850_df = dust_salt_mean_1850_df.dropna(how='all')
dust_salt_error_1850_df = dust_salt_error_1850_df.dropna(how='all')
# Convert values to strings
dust_salt_mean_1850_df = dust_salt_mean_1850_df.applymap(lambda x: '%+.2f' % x)
dust_salt_error_1850_df = dust_salt_error_1850_df.applymap(lambda x: '%.2f' % x)
# Combine mean and standard error strings
dust_salt_1850_df = dust_salt_mean_1850_df + ' ± ' + dust_salt_error_1850_df
# Replace 'nan ± nan' with ''
dust_salt_1850_df = dust_salt_1850_df.replace('nan ± nan', '')
dust_salt_1850_df = dust_salt_1850_df.replace('+nan ± nan', '')
# Rename columns
dust_salt_1850_df = dust_salt_1850_df.rename(index=str, columns={'tSSLT': 'Sea salt aerosol',
                                                                 'tDST': 'Dust aerosol',})
# Save to CSV
out_filename = 'csv/tS04.csv'
dust_salt_1850_df.to_csv(out_filename)
print('Written {}'.format(out_filename))
# Save to Excel
out_filename = 'csv/tS04.xlsx'
dust_salt_1850_df.to_excel(out_filename)
print('Written {}'.format(out_filename))
# Display
dust_salt_1850_df.style.set_caption('Year-1850')
! date<jupyter_output>Fri Aug 17 17:08:04 +08 2018
{'fraction_non_alphanumeric': 0.13135400333822264, 'fraction_numerical': 0.08208992969500783, 'mean_word_length': 3.2256892498397094, 'pattern_counts': {'":': 0, '<': 30, '<?xml version=': 0, '>': 34, 'https://': 7, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 56, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '637635', 'license': 'MIT', 'n_tokens_mistral': 8555, 'n_tokens_neox': 7263, 'n_words': 1451, 'path': 'analysis_draft2017b/budget_tables_draft2017b.ipynb', 'repo_name': 'grandey/p17c-marc-comparison', 'size': '19796'}