Análise de dados Multi-paramétrica

04 de Novembro de 2013

Prof. Filipe Fernandes

Podemos classificar a maioria das propriedades oceânicas em traçadores do tipo:

  • conservativos e;
  • não-conservativos.

Análise de Massas de Água Clássica

  • Utilizamos apenas traçadores conservativos;
  • Águas tipos são caracterizadas por pares $\theta$-S;
  • Curva de mistura de duas ou mais águas tipos formando o diagrama $\theta$-S.
In [2]:
#Mistura simples de duas águas tipo.
Image('./figures/massa_dagua_mistura_01.png', retina=True)
Out[2]:
In [3]:
# Mistura de três águas tipo.
Image('./figures/massa_dagua_mistura_02.png', retina=True)
Out[3]:
In [4]:
# Detalhe da formação de um diagrama TS com uma mistura de três águas tipo.
Image('./figures/massa_dagua_mistura_03.png', retina=True)
Out[4]:

Problemas com a análise clássica:

  • Depende apenas de T e S, logo é incapaz de descrever massas d'água com características TS próximas porém de origems diferentes (Ex.: A mistura entre duas massas d'água gera uma terceira que se confunde facilmente com as massas que a originaram).
  • Podemos calcular as porcentagens de mistura apenas se conhecermos as águas tipo. Mas em geral a mistrua será subestimada por não contabilizarar massas d'água que não podemos descrever com T e S apenas.
  • Não traz informação sobre a idade da massa d'água.
In [5]:
%%latex
Modelo da Análise Clássica:
\begin{align*}
x_1 T_1 + x_2 T_2 + x_3 T_3 = T \\
x_1 S_1 + x_2 S_2 + x_3 S_3 = S \\
x_1 + x_2 + x_3 = 1 \\
\text{ou} \\
\mathbf{Ax} = \mathbf{b}
\end{align*}
Modelo da Análise Clássica: \begin{align*} x_1 T_1 + x_2 T_2 + x_3 T_3 = T \\ x_1 S_1 + x_2 S_2 + x_3 S_3 = S \\ x_1 + x_2 + x_3 = 1 \\ \text{ou} \\ \mathbf{Ax} = \mathbf{b} \end{align*}

Análise de dados Multiparamétrica

  • Atualmente, a análise de massas de água utiliza de dados que não eram disponíveis a 60 anos atrás, quando Sverdrup et al. (1942) formalizaram o estudo deste tema.
  • A combinação entre conjuntos de dados hidrográficos, que incluem os parâmetros não-conservativos da água do mar e os dados hidrográficos clássicos (T e S), tornou possível o desenvolvimento de métodos inversos capazes de extrair um maior número informações sobre as massas de água.
  • Inclui-se na categoria destes novos métodos, a OMPA (Sigla em Inglês para "Análise Multiparamétrica Ótima"), que possui uma abordagem efetivamente quantitativa.

Mas como incluir elementos não conservativos nas equações acima?

  • Nutrientes são introduzidos em sub-superfície por processos de oxidação de matéria orgânica;
  • Isso ocorre em proporções estequiométricas com o comsumo de oxigênio;
  • Essa proporção é chamada de Taxa de Redfield.
In [6]:
%%latex
\begin{align*}
NO &= 9\text{NO}_3 + \text{O}_2 \\
PO &= 135\text{PO}_4 + \text{O}_2
\end{align*}
\begin{align*} NO &= 9\text{NO}_3 + \text{O}_2 \\ PO &= 135\text{PO}_4 + \text{O}_2 \end{align*}
  • Adicionalmente, Silicato pode ser considerado quase conservativo no Atlântico. Elevando de 2 para 5 o número de parâmetros que termos para descrever massas de água.

Infelizmente não podemos aplicar a equação acima diretamente. Isso porque:

  • Existem erros nas amostras;
  • Águas tipos desconhecidas,;
  • Flutações nas regiões fonte e etc

Isso nos leva a um sistema que raramente fecha em 1 como esperado.

Por isso temos que utilizar uma solução por ajustes de mínimos quadrados do tipo:

$$D^2 = (\mathbf{Ax-b})'\mathbf{W}'\mathbf{W}(\mathbf{Ax-b})$$

Os passos para se aplicar uma OMPA?

  1. Identifique um conjunto de $M-1$ traçadores conservativos;
  1. Defina as características (ou seja os valores dos M-1 traçadores) para as suas águas tipo;
  1. Desenhe a sua matriz e os vetores de dados ($\mathbf{Ax-b}$);
  1. Padronize e crie pesos para a sua matriz $\mathbf{A}$ e para o conjunto de dados.
  1. Busque a solução inversa para $x$ para o seu conjunto de dados por Mínimos Quadrados do modelo $\mathbf{\tilde{G}x-\tilde{d}}$ sujeito a uma constrição positiva!
  • Os passos Oceanográficos mais iportantes estão nos passas 1, 2 e 4.
  • É essencial checar a distribuição espacial dos resíduos do modelo de Mínimo quadrados ($\mathbf{\tilde{G}x-\tilde{d}}$).

Exemplo de OMPA

In [7]:
# Load bottle data: http://www.nodc.noaa.gov/woce/woce_v3/wocedata_1/whp/
# These files are mMissing C14 data.

path = './WOCE/BOTTLE'
kw = dict(skiprows=5, skipfooter=1, skipinitialspace=True, header=0)

# http://www.nodc.noaa.gov/woce/woce_v3/wocedata_1/whp/data/onetime/atlantic/a16/a16s/index.htm
A16S = read_csv('%s/a16s_hy1.csv' % path, **kw)
A16S_units = A16S.iloc[0]
A16S = A16S.iloc[1:]

# http://www.nodc.noaa.gov/woce/woce_v3/wocedata_1/whp/data/onetime/atlantic/a16/a16c/index.htm
A16C = read_csv('%s/a16c_hy1.csv'  % path, **kw)
A16C_units = A16C.iloc[0]
A16C = A16C.iloc[1:]

# http://www.nodc.noaa.gov/woce/woce_v3/wocedata_1/whp/data/onetime/atlantic/a16/a16n/index.htm
A16N = read_csv('%s/a16n_hy1.csv' % path, **kw)
A16N_units = A16N.iloc[0]
A16N = A16N.iloc[1:]

# Merge all 3 section into one.
A16 = A16S.merge(A16C, how='outer').merge(A16N, how='outer')
/home/filipe/.virtualenvs/iris/lib/python2.7/site-packages/pandas/io/parsers.py:639: ParserWarning: Falling back to the 'python' engine because the 'c' engine does not support skip_footer; you can avoid this warning by specifying engine='python'.
  ParserWarning)
In [8]:
# Select a continous track.
maskN = np.logical_and(A16['STNNBR'] >= 1, A16['STNNBR'] <= 124)
maskS = np.logical_and(A16['STNNBR'] >= 237, A16['STNNBR'] <= 278)
maskC = np.logical_and(A16['STNNBR'] >= 314, A16['STNNBR'] <= 373)
mask = maskN | maskS | maskC
A16 = A16[mask]

# Use time as index.
dtime = [datetime.strptime('%s %03d' % (int(date), int(time)), '%Y%m%d %H%M') for
         date, time in zip(A16['DATE'], A16['TIME'])]

# Drop unused and non-numeric columns.
A16.drop(['DATE','TIME', 'EXPOCODE', 'SECT_ID'], axis=1, inplace=True)
A16.index = dtime

# Force all values to be numeric.
A16 = A16.convert_objects(convert_numeric=True)
In [9]:
# Select good data using flags.
# Multiple NA is used. -999.00, -999.9 and etc,
# so I'm getting all data that is marked as less than -999.
A16[A16 <= -999] = np.NaN

def apply_flag(flag):
    prop = flag.split('_')[0]
    mask = A16[flag] == 2.
    A16[prop][mask] = np.NaN

flags = ['BTLNBR_FLAG_W', 'CTDSAL_FLAG_W', 'SALNTY_FLAG_W', 'SILCAT_FLAG_W',
         'OXYGEN_FLAG_W', 'NITRAT_FLAG_W', 'NITRIT_FLAG_W', 'PHSPHT_FLAG_W',
         'CFC-11_FLAG_W', 'CFC-12_FLAG_W', 'TCARBN_FLAG_W', 'PCO2_FLAG_W']

for flag in flags:
    # Not sure if "exchange" format already applied bad flag.
    # So I'm just dropping that column.
    #apply_flag(flag)
    A16.drop([flag], axis=1, inplace=True)
In [10]:
# Compute AOU.
import seawater as sw
from oceans.sw_extras import o2sat
Osat = o2sat(A16['CTDSAL'].values, A16['THETA'].values)

A16['AOU'] = Osat - A16['OXYGEN']

# Normalize Total DIC.
A16['TCARBN_NORM'] = A16['TCARBN'] * (35 / A16['CTDSAL'])

pt = sw.ptmp(A16['CTDSAL'], A16['CTDTMP'], A16['CTDPRS'], pr=0)
A16['PDEN'] = sw.pden(A16['CTDSAL'], A16['CTDTMP'], A16['CTDPRS'], pr=0) - 1000

pandas.set_option('display.max_columns', 40)
A16.head(5)
Out[10]:
STNNBR CASTNO SAMPNO BTLNBR LATITUDE LONGITUDE DEPTH CTDPRS CTDTMP CTDSAL SALNTY OXYGEN SILCAT NITRAT NITRIT PHSPHT CFC-11 CFC-12 THETA TCARBN PCO2 PCO2TMP CTDOXY AOU TCARBN_NORM PDEN
1989-02-01 13:45:00 237 1 1 1 -32.4933 -24.985 4148 3.4 22.4242 35.6509 NaN 226.1 1.47 0.29 0.00 0.09 1.930 1.011 22.4235 NaN NaN NaN NaN -11.414709 NaN 24.592728
1989-02-01 13:45:00 237 1 2 2 -32.4933 -24.985 4148 22.0 21.8333 35.6469 NaN 229.4 1.47 0.29 0.00 0.08 NaN NaN 21.8290 NaN NaN NaN NaN -12.456727 NaN 24.757427
1989-02-01 13:45:00 237 1 3 3 -32.4933 -24.985 4148 51.1 18.6389 35.5794 NaN 251.7 1.47 0.29 0.00 0.11 2.242 1.150 18.6299 NaN NaN NaN NaN -21.700186 NaN 25.558316
1989-02-01 13:45:00 237 1 4 4 -32.4933 -24.985 4148 63.1 17.3043 35.5762 NaN 256.9 1.37 0.29 0.00 0.14 NaN NaN 17.2938 NaN NaN NaN NaN -21.029268 NaN 25.886125
1989-02-01 13:45:00 237 1 5 5 -32.4933 -24.985 4148 83.5 16.2117 35.6145 NaN 244.2 1.37 0.39 0.02 0.18 2.274 1.159 16.1983 NaN NaN NaN NaN -3.363767 NaN 26.174684
In [11]:
import os
from ctd.ctd import CTD
from glob import glob
from pandas import Panel
from collections import OrderedDict

def get_pattern(fname, pattern):
    with open(fname) as f:
        for line in f.readlines():
            if line.startswith(pattern):
                value = float(line.split('=')[1])
                continue
        return value
    
    
names = ['CTDPRS', 'CTDTMP', 'CTDSAL', 'CTDOXY']
def load_ctd_section(fname):
    kw = dict(skiprows=20, skipfooter=1, names=names, usecols=(0, 2, 4, 6), index_col='CTDPRS')
    df = read_csv(fname, **kw)
    return CTD(df)

# Using only contiguous track.
maskN = range(1, 124 + 1)
maskS = range(237, 278 + 1)
maskC = range(314, 373 + 1)
mask = maskN + maskS + maskC

path = './WOCE/CTD'
A16S = sorted(glob('%s/a16s*.csv' % path), reverse=True)
A16C = sorted(glob('%s/a16c*.csv' % path), reverse=False)
A16N = sorted(glob('%s/a16n*.csv' % path), reverse=True)
fnames = A16S + A16C + A16N

lon, lat, CTD_A16 = [], [], OrderedDict()
for fname in fnames:
    # Select a continous track.
    STNNBR = int(fname.split('_')[1])
    if STNNBR in mask:
        name = os.path.basename(fname).split('.')[0]
        lon.append(get_pattern(fname, 'LONGITUDE'))
        lat.append(get_pattern(fname, 'LATITUDE'))
        CTD_A16.update({STNNBR: load_ctd_section(fname)})
    
CTD_A16 = Panel.fromDict(CTD_A16)


CTDTMP = CTD(CTD_A16.minor_xs('CTDTMP'))
CTDTMP.lon, CTDTMP.lat = lon, lat

CTDSAL = CTD(CTD_A16.minor_xs('CTDSAL'))
CTDSAL.lon, CTDSAL.lat = lon, lat


from oceans.sw_extras import gamma_GP_from_SP_pt
p = CTDSAL.index.values.astype(float)
s, t, p, lon, lat = np.broadcast_arrays(CTDSAL.values, CTDTMP.values, p[..., None], lon, lat)
pt = sw.ptmp(s, t, p, pr=0)
#gamma = gamma_GP_from_SP_pt(s, pt, p, lon, lat)
pden = sw.pden(s, t, p, pr=0) - 1000
In [12]:
# Get bottles x, y.
from pandas import pivot_table

LON = pivot_table(A16, values='LONGITUDE', rows=['CTDPRS'], cols=['STNNBR'])
LAT = pivot_table(A16, values='LATITUDE', rows=['CTDPRS'], cols=['STNNBR'])

# Get just common stations.
common = []
for stn in CTDSAL.columns:
    if stn in LON.columns:
        common.append(stn)
common = np.array(common, np.float_)
LON = LON[common]
LAT = LAT[common]

DIST = np.r_[0, sw.dist(LAT.mean(), LON.mean())[0].cumsum()]
bottles = DataFrame(LON.notnull().values, index=LON.index, columns=DIST)

X, Y = np.meshgrid(bottles.columns.values, bottles.index.values)
X, Y = X[bottles.values], Y[bottles.values]

# Make salinity section.
lon_0 = 360 + A16['LONGITUDE'].mean()
lat_0 = A16['LATITUDE'].mean()

from mpl_toolkits.axes_grid1.inset_locator import inset_axes
def make_basemap(projection='ortho', resolution='c', ax=None):
    m = Basemap(projection=projection, resolution=resolution,
                lon_0=lon_0, lat_0=lat_0, ax=ax)
    m.drawcoastlines()
    m.fillcontinents(color='0.85')
    return fig, m

from ctd.plotting import plot_section
levels = np.arange(33.8, 37.4, 0.1)
fig, ax, cbar = plot_section(CTDSAL, cmap=cm.odv, marker=None,
                             levels=levels, figsize=(8, 4))
ax.set_ylabel('Depth [m]')
ax.set_xlabel('Distance in [km]')

axin = inset_axes(ax, width="40%", height="40%", loc=4)
fig, m = make_basemap(ax=axin)
kw = dict(marker='.', color='#FF9900', linestyle='none')
m.plot(*m(A16['LONGITUDE'].values, A16['LATITUDE'].values), **kw)
ax.plot(X, Y, 'k.', alpha=0.5, markersize=2)

levels = [22.5, 23, 24.5, 25.5, 26.5, 27, 27.5, 27.845, 27.9]
levels = [25, 26.5] + [27, 27.5] + [27.7, 27.85]
cs = ax.contour(DIST, CTDSAL.index.values.astype(float), ma.masked_invalid(pden), levels=levels, colors='k')
_ = ax.clabel(cs, fmt='%2.2f')
/home/filipe/.virtualenvs/iris/lib/python2.7/site-packages/pandas/util/decorators.py:81: FutureWarning: the 'cols' keyword is deprecated, use 'columns' instead
  warnings.warn(msg, FutureWarning)
/home/filipe/.virtualenvs/iris/lib/python2.7/site-packages/pandas/util/decorators.py:81: FutureWarning: the 'rows' keyword is deprecated, use 'index' instead
  warnings.warn(msg, FutureWarning)
/home/filipe/.virtualenvs/iris/lib/python2.7/site-packages/matplotlib/text.py:52: UnicodeWarning: Unicode equal comparison failed to convert both arguments to Unicode - interpreting them as being unequal
  if rotation in ('horizontal', None):
/home/filipe/.virtualenvs/iris/lib/python2.7/site-packages/matplotlib/text.py:54: UnicodeWarning: Unicode equal comparison failed to convert both arguments to Unicode - interpreting them as being unequal
  elif rotation == 'vertical':
In [13]:
# Seção de salinidade no Atlântico.
Image('./figures/Atlantic_salinity.png', width=2206/4., height=1122/4., retina=False)
Out[13]:
In [15]:
df = A16[['PHSPHT', 'NITRAT', 'OXYGEN']].dropna()
sckw = dict(cmap=cm.odv, alpha=0.7, vmax=350, vmin=0)

fig, ax = plt.subplots()
cs = ax.scatter(df['PHSPHT'], df['NITRAT'], c=df['OXYGEN'], **sckw)
fig.colorbar(cs, **cbarkw)
ax.set_title(r'Oxygen [$\mu$mol kg$^{-1}$]')
ax.set_xlabel(r'Phosphate [$\mu$mol kg$^{-1}$]')
ax.set_ylabel(r'Nitrate [$\mu$mol kg$^{-1}$]')
ax.axis('tight')

plot_gmfit(df['PHSPHT'], df['NITRAT'], ax)
In [16]:
Image('./figures/nvp_01.png', retina=True)
Out[16]: