Contents

Seasonal amplitude of vertical gradient

Seasonal amplitude of vertical gradient

Examine aircraft observations relative to deseasonalized record at South Pole Observatory (SPO).

%load_ext autoreload
%autoreload 2

import os
import numpy as np
import xarray as xr
xr.set_options(display_style='text')

import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec

import datasets
import emergent_constraint as ec
import figure_panels
import obs_aircraft
import obs_surface
import regression_models
import util

Aircraft profiles

dsets_prof = datasets.aircraft_profiles('obs')[['co2_med']]
dsets_prof

<xarray.Dataset>
Dimensions:    (profile: 361, theta: 27)
Coordinates:
    campaign   (profile) <U32 'HIPPO-1' 'HIPPO-1' ... 'ORCAS-F' 'ORCAS-F'
    doy        (profile) float64 20.0 20.0 20.0 20.0 20.0 ... 56.0 56.0 0.0 0.0
    flight_id  (profile) <U32 'HIPPO-001-007' ... 'ORCAS-001-019'
    lat        (profile) float64 -44.73 -46.5 -49.75 ... -54.85 -51.65 -45.44
    lon        (profile) float64 170.4 169.6 170.1 ... -68.29 -72.46 -76.02
    month      (profile) float64 1.0 1.0 1.0 1.0 1.0 1.0 ... 2.0 2.0 2.0 2.0 2.0
  * profile    (profile) <U17 'HIPPO-001-007-074' ... 'ORCAS-001-019-207'
  * theta      (theta) float64 270.0 275.0 280.0 285.0 ... 390.0 395.0 400.0
    time       (profile) datetime64[ns] 2009-01-20 2009-01-20 ... 2016-02-29
    year       (profile) float64 2.009e+03 2.009e+03 ... 2.016e+03 2.016e+03
Data variables:
    co2_med    (profile, theta) float64 nan nan nan nan nan ... nan nan nan nan
spo_trend = xr.DataArray(
    np.interp(dsets_prof.time, spo_fit.time, spo_fit.trend),
    dims=('profile'), 
)
dco2 = dsets_prof.co2_med - spo_trend
dco2

<xarray.DataArray (profile: 361, theta: 27)>
array([[nan, nan, nan, ..., nan, nan, nan],
       [nan, nan, nan, ..., nan, nan, nan],
       [nan, nan, nan, ..., nan, nan, nan],
       ...,
       [nan, nan, nan, ..., nan, nan, nan],
       [nan, nan, nan, ..., nan, nan, nan],
       [nan, nan, nan, ..., nan, nan, nan]])
Coordinates:
  * profile    (profile) object 'HIPPO-001-007-074' ... 'ORCAS-001-019-207'
    campaign   (profile) <U32 'HIPPO-1' 'HIPPO-1' ... 'ORCAS-F' 'ORCAS-F'
    doy        (profile) float64 20.0 20.0 20.0 20.0 20.0 ... 56.0 56.0 0.0 0.0
    flight_id  (profile) <U32 'HIPPO-001-007' ... 'ORCAS-001-019'
    lat        (profile) float64 -44.73 -46.5 -49.75 ... -54.85 -51.65 -45.44
    lon        (profile) float64 170.4 169.6 170.1 ... -68.29 -72.46 -76.02
    month      (profile) float64 1.0 1.0 1.0 1.0 1.0 1.0 ... 2.0 2.0 2.0 2.0 2.0
  * theta      (theta) float64 270.0 275.0 280.0 285.0 ... 390.0 395.0 400.0
    time       (profile) datetime64[ns] 2009-01-20 2009-01-20 ... 2016-02-29
    year       (profile) float64 2.009e+03 2.009e+03 ... 2.016e+03 2.016e+03

Seasonal amplitude in the column

Use a harmonic function to estimate seasonal amplitude.

xhat = np.arange(-5, 365+5, 1)    
yhat = xr.DataArray(np.ones((len(dco2.theta), len(xhat)))*np.nan, 
                    dims=('theta', 'doy'),
                    coords=dict(
                        theta=dco2.theta,
                        doy=xhat,
                    ),
                   )

doy = dco2.doy.values
for k in range(len(dco2.theta)):
    x, y = util.antyear_daily(doy, dco2.isel(theta=k).values)
    missing = np.isnan(x) | np.isnan(y)
    if np.sum(~missing) < 5:
        continue
    p, pcov = ec.curve_fit(
        ec.harmonic, 
        xdata=x[~missing]/365.25, 
        ydata=y[~missing], 
    )

    yhat.data[k, :] = ec.harmonic(xhat/365.25, *p)

    
seasonal_amplitude = (yhat.max('doy') - yhat.min('doy'))
yhat.plot()

/glade/work/mclong/miniconda3/envs/so-co2/lib/python3.7/site-packages/scipy/optimize/minpack.py:834: OptimizeWarning: Covariance of the parameters could not be estimated
  category=OptimizeWarning)

<matplotlib.collections.QuadMesh at 0x2b0a91513c90>
_images/gradients-seasonal-amplitude_14_2.png

Hovmöller visualization

campaign_info = obs_aircraft.get_campaign_info(verbose=False, lump_orcas=True)

fig = plt.figure(figsize=(8, 4)) #dpi=300)

c_spec = figure_panels.marker_spec_campaigns(lump_orcas=True)
txt_box_props = dict(boxstyle='square,pad=0', facecolor='none', edgecolor='none')

# set up plot grid
gs = gridspec.GridSpec(
    nrows=1, ncols=4, 
    width_ratios=[0.65, 0.015, 0.075, 0.375],
    left=0., right=1.,
    bottom=0.05, top=0.95,
    hspace=0.25, wspace=0.05,
)
axs = np.empty((1, 2)).astype(object)
axs[0, 0] = plt.subplot(gs[0, 0])
cax = plt.subplot(gs[0, 1])
axs[0, 1] = plt.subplot(gs[0, 3])

ax = axs[0, 0] 
pc = ax.contourf(
              yhat.doy,
              yhat.theta.sel(theta=slice(None, 320)), 
              yhat.sel(theta=slice(None, 320)),
              levels=figure_panels.levels,
              norm=figure_panels.divnorm,
              cmap=figure_panels.cmap,
             )

for c, info in campaign_info.items():
    tb = util.day_of_year(info['time_bound'])
    x, _ = util.antyear_daily(tb, np.ones(2))
    ax.axvspan(
        x[0], x[1], 
        ymin=0.84 if c == 'ORCAS' else 0.86, 
        ymax=0.98,
        color=c_spec[c]['color'], 
        alpha=1 if c in ['HIPPO-1', 'ATOM-2'] else 0.75, 
    )
    xtxt = x.mean() + np.diff(x)*0.18
    if c == 'ORCAS':
        xtxt += 10
    ax.text(
        xtxt, 320., c,
        rotation=90, 
        ha='center', 
        verticalalignment='top', 
        color='k', 
        fontsize=7,
        bbox=txt_box_props,
    )
    
    
ax.set_ylim((269., 321.))
ax.set_title('Interpolated obs: CO$_2$ minus SPO trend')


cb = plt.colorbar(pc, cax=cax)
cb.ax.set_title('$\Delta$CO$_2$ [ppm]      ', loc='center')

ax.set_ylabel('$\\theta$ [K]')

ax.set_xlim((-10, 375))
ax.set_xticks(figure_panels.bomday)
ax.set_xticklabels([f'        {m}' for m in figure_panels.monlabs_ant]+['']);

ax = axs[0, 1]
ax.plot(
    seasonal_amplitude.sel(theta=slice(None, 320)), 
    yhat.theta.sel(theta=slice(None, 320)), 
    '.-', label='Seasonal amplitude')

ax.set_yticklabels([]);
ax.set_xlabel('Seasonal amplitude [ppm]')
ax.axvline(seasonal_amplitude.sel(theta=300.), lw=0.5, c='dimgray', linestyle='--')
ax.set_ylim((269., 321.))
ax.set_title('Amplitude of seasonal cycle')
util.label_plots(fig, [ax for ax in axs.ravel()], xoff=-0.02)

util.savefig('seasonal-amplitude')
_images/gradients-seasonal-amplitude_16_0.png

previous

Gradient metrics

next

Trend in horizontal gradient metric