<jupyter_start><jupyter_text>Flux estimates for each aircraft campaign<jupyter_code>%load_ext autoreload
%autoreload 2
import cftime
import numpy as np
import xarray as xr

import emergent_constraint as ec
import util
clobber = False
clobber_deep = False
air_parms = ec.get_parameters('default')

obj_air = {};
acs = {}; scs = {};
for constraint_type in ['ocean_constraint',]:
    model_input_lists = ec.get_model_tracer_lists(constraint_type)
    obj_air[constraint_type] = ec.whole_enchilada(clobber=clobber_deep, **model_input_lists)    
    acs[constraint_type] = obj_air[constraint_type].get_ac(**air_parms, clobber=clobber)
ac = acs['ocean_constraint']
dfo = ac.gradient_summary[['vg_obs', 'vg_obs_std', 'vg_ext', 'vg_ext_std',]]
dfo = dfo.rename(
    columns=dict(
        vg_obs='Observed Vertical Gradient [ppm]',
        vg_obs_std='Observed Gradient Uncertainty [ppm]',
        vg_ext='Simulated Land+Fossil Gradient [ppm]',
        vg_ext_std='Simulated Land+Fossil Gradient [ppm]',        
    )
)

# add columns for fluxes, flux errors, and start/end of 90-day windows?
dfo['Flux Date Start'] = [str(d[0]) for d in ac.campaign_flux.loc[dfo.index].flux_window_date_bound]
dfo['Flux Date End'] = [str(d[1]) for d in ac.campaign_flux.loc[dfo.index].flux_window_date_bound]
dfo['Flux Estimate [Pg C/yr]'] = ac.campaign_flux.loc[dfo.index].flux
dfo['Flux Uncertainty [Pg C/yr]'] = ac.campaign_flux.loc[dfo.index].flux_error

dfo
dfo2 = ac.campaign_flux[['doy_mid', 'flux', 'flux_error']].rename(
    columns={
    'doy_mid': 'Day of Year',
    'flux': 'Flux Estimate [Pg C/yr]',
    'flux_error': 'Flux Uncertainty [Pg C/yr]',
})

# add columns for fluxes, flux errors, and start/end of 90-day windows?
dfo2['Flux Date Start'] = [str(d[0]) for d in ac.campaign_flux.loc[dfo2.index].flux_window_date_bound]
dfo2['Flux Date End'] = [str(d[1]) for d in ac.campaign_flux.loc[dfo2.index].flux_window_date_bound]
dfo2['Date'] = [d[0] + (d[1] - d[0])/2 for d in ac.campaign_flux.loc[dfo2.index].flux_window_date_bound]
dfo2 = dfo2[['Date', 'Flux Date Start', 'Flux Date End', 'Flux Estimate [Pg C/yr]', 'Flux Uncertainty [Pg C/yr]']]

comment_block = """# Southern Ocean flux estimates from aircraft observations described in Long et al. [2021]:
#   <NAME>, <NAME>, <NAME>, <NAME>, <NAME>, <NAME>, <NAME>, 
#     <NAME>, <NAME>, <NAME>, <NAME>, <NAME>, <NAME>, <NAME>, 
#     <NAME>, <NAME>, <NAME>, <NAME>, <NAME>, <NAME>, <NAME>, <NAME>, 
#     <NAME>, Strong Southern Ocean carbon uptake evident in airborne observations, 
#     Science (2021), https://doi.org/10.1126/science.abi4355.
#
# Flux estimates are valid for the region south of 45°S and over a period of 90 days
# prior to the end of each aircraft campaign.
#
# Column descriptions:
# Date: Mid-point of the 90 day period over which flux estimate applies
# Flux Date Start: start of 90 day period
# Flux Date End: end of 90 day period
# Flux Estimate [Pg C/yr]: the flux estimate for each campaign
# Flux Uncertainty [Pg C/yr]: the uncertainty in the flux estimate
#
"""
with open('flux-results/Long-etal-2021-aircraft-campaign-fluxes.csv', 'w') as fid:
    fid.write(comment_block)
    dfo2.to_csv(fid)
    
dfo2    
abcd, pcov = ac.harmonic_fit

x = np.arange(1., 366., 1.) - 0.5 
y = ec.harmonic(x/365.25, *abcd) 

time = np.concatenate((
    cftime.date2num(xr.cftime_range('2009-07-01', '2009-12-31', freq='1D').values, units='days since 2009-01-01'),
    cftime.date2num(xr.cftime_range('2009-01-01', '2009-06-30', freq='1D').values, units='days since 2009-01-01'),
)) + 0.5

ndx = np.argsort(time)
assert (x == time[ndx]).all()
y = y[ndx]

climatology_bounds = xr.DataArray(
    np.vstack((
        time[ndx] - 0.5, 
        time[ndx] + 0.5 + cftime.date2num(cftime.datetime(2018, 1, 1), units='days since 2009-01-01'))).T,
    dims=('time', 'd2'),
    attrs={'long_name': 'nominal bounds for climatology'}
)        

dso = xr.Dataset(
    {'time': xr.DataArray(
        time[ndx],
        dims=('time'), 
        name='time',
    ),
     'flux': xr.DataArray(
        y,
        dims=('time'),
        coords={'time': x},
        name='fgco2',
        attrs={'long_name': 'Air-sea CO$_2$ flux (south of 45°S)', 
               'units': 'Pg C yr$^{-1}$',
              },
     ),
    }
)

dso["climatology_bounds"] = climatology_bounds
dso.time.attrs = {'units': 'days since 2009-01-01', 'long_name': 'time', 'climatology': 'climatology_bounds'}
dso.time.encoding['_FillValue'] = None

dso.attrs['source'] = 'flux climatology generated as a harmonic fit to Southern Ocean aircraft campaign CO2 data'
dso.attrs['reference'] = 'Long et al., Science (2021) https://doi.org/10.1126/science.abi4355'

dso
dso.flux.plot()
dso.to_netcdf('flux-results/Long-etal-2021-seasonal-flux-fit.nc')
dso.info()
with xr.open_dataset('flux-results/Long-etal-2021-seasonal-flux-fit.nc') as ds:
    ds.flux.plot()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1636697997548018, 'fraction_numerical': 0.03637106661217818, 'mean_word_length': 4.056818181818182, 'pattern_counts': {'":': 0, '<': 28, '<?xml version=': 0, '>': 28, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '190256', 'license': 'MIT', 'n_tokens_mistral': 2019, 'n_tokens_neox': 1842, 'n_words': 434, 'path': 'so-co2-airborne-obs/fluxes-aircraft-campaign-table.ipynb', 'repo_name': 'NCAR/so-co2-airborne-obs', 'size': '4994'}