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*}
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*}
- 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?
- Identifique um conjunto de \(M-1\) traçadores conservativos;
- Defina as características (ou seja os valores dos M-1 traçadores) para as suas águas tipo;
- Desenhe a sua matriz e os vetores de dados (\(\mathbf{Ax-b}\));
- Padronize e crie pesos para a sua matriz \(\mathbf{A}\) e para o conjunto de dados.
- 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')
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]:
In [11]:
import os
from 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')
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]:
In [17]:
df = A16[['PHSPHT', 'TCARBN_NORM', 'AOU']].dropna()
sckw = dict(cmap=cm.odv, alpha=0.7, vmax=300, vmin=-25)
fig, ax = plt.subplots()
cs = ax.scatter(df['PHSPHT'], df['TCARBN_NORM'], c=df['AOU'], **sckw)
fig.colorbar(cs, **cbarkw)
ax.grid(True)
ax.set_title(r'AOU [$\mu$mol kg$^{-1}$]')
ax.set_xlabel(r'Phosphate [$\mu$mol kg$^{-1}$]')
ax.set_ylabel(r'nTDIC (Sal normalized) [$\mu$mol kq$^{-1}$]')
ax.axis('tight')
plot_gmfit(df['PHSPHT'], df['TCARBN_NORM'], ax)
In [18]:
# Fosfato vs Carbono Orgânico Dissolvido (oxigênio em cores).
Image('./figures/ntdicvp_01.png', retina=True)
Out[18]:
In [19]:
df = A16[['PHSPHT', 'AOU']].dropna()
fig, ax = plt.subplots()
ax.scatter(df['PHSPHT'], df['AOU'], c='k', alpha=0.7)
ax.grid(True)
ax.set_xlabel(r'Phosphate [$\mu$mol kg$^{-1}$]')
ax.set_ylabel(r'AOU [$\mu$mol kg$^{-1}$]')
ax.axis('tight')
plot_gmfit(df['PHSPHT'], df['AOU'], ax)
In [20]:
Image('./figures/aouvsp.png', retina=True)
Out[20]:
In [21]:
df = A16[['NITRAT', 'AOU']].dropna()
fig, ax = plt.subplots()
ax.scatter(df['NITRAT'], df['AOU'], c='k', alpha=0.7)
ax.grid(True)
ax.set_xlabel(r'Nitrate [$\mu$mol kg$^{-1}$]')
ax.set_ylabel(r'AOU [$\mu$mol kg$^{-1}$]')
ax.axis('tight')
plot_gmfit(df['NITRAT'], df['AOU'], ax)
In [22]:
Image('./figures/aouvn.png', retina=True)
Out[22]:
In [23]:
df = A16[['TCARBN_NORM', 'AOU']].dropna()
fig, ax = plt.subplots()
ax.scatter(df['TCARBN_NORM'], df['AOU'], c='k', alpha=0.7)
ax.grid(True)
ax.set_xlabel(r'nTDIC (Sal normalized) [$\mu$mol kq$^{-1}$]')
ax.set_ylabel(r'AOU [$\mu$mol kg$^{-1}$]')
ax.axis('tight')
plot_gmfit(df['TCARBN_NORM'], df['AOU'], ax)
In [24]:
Image('./figures/auovdic.png', retina=True)
Out[24]:
A relação estequiométrica C:AOU:N:P é 110.47:66.48:15.3:1. Mas apenas Nitrato vs Fosfato e DIC vs Fosfato mostraram uma relação linear.
AOU:P não mostra a mesma tendência linear devido a alta concentração de fosfato pré-formado na superfície.
Nutrientes pré-formados = Total de nutrientes - Nutrientes Regenerados
Nutrientes pré-formados são característicos da região de onde se original, podemos ser usados como traçadores de massas d'água!
Como se espera que a concentração de nutrientesaumente com o aumento do AOU todos os nutrientes que estão presentes em uma concentração maior que zero em AOU zero podem ser chamados de pré-formados.
Utilização Aparente de Oxigênio (UAO) [Apparent Oxygen Utilisation (AOU)]
É a diferença entre a concentração medida de Oxigênio dissolvido e a concentração do seu Estado de Equilíbrio de Saturação
Tais diferenção ocorrem por atividade biológica que muda o oxigênio disponível na água.
Consequentemte AOU de uma amostra de água representa a some da atividade biológica que essa amostra sofreu desde seu equilíbrio com a atmosfera.
Para melhor enxergar o que são os nutrientes pré-formados vamos utilizar apenas as camadas superficiais.
In [25]:
df = A16[['PHSPHT', 'AOU', 'PDEN']].dropna()
mask = np.logical_and(df['PDEN'] >=25, df['PDEN'] <= 26.7)
df = df[mask]
fig, ax = plt.subplots()
ax.scatter(df['PHSPHT'], df['AOU'], c='k', alpha=0.7)
ax.grid(True)
ax.set_xlabel(r'Phosphate [$\mu$mol kg$^{-1}$]')
ax.set_ylabel(r'AOU [$\mu$mol kg$^{-1}$]')
ax.axis('tight')
plot_gmfit(df['PHSPHT'], df['AOU'], ax)
In [26]:
# AOU vs Fosfato apenas das camadas superficiais.
Image('./figures/aouvp.png', retina=True)
Out[26]:
In [27]:
df = A16[['NITRAT', 'AOU', 'PDEN']].dropna()
mask = np.logical_and(df['PDEN'] >=25, df['PDEN'] <= 26.7)
df = df[mask]
fig, ax = plt.subplots()
ax.scatter(df['NITRAT'], df['AOU'], c='k', alpha=0.7)
ax.grid(True)
ax.set_xlabel(r'Nitrate [$\mu$mol kg$^{-1}$]')
ax.set_ylabel(r'AOU [$\mu$mol kg$^{-1}$]')
ax.axis('tight')
plot_gmfit(df['NITRAT'], df['AOU'], ax)
In [28]:
# AOU vs Carbônico orgânico Dissolvido apenas das camadas superficiais.
Image('./figures/aouvdic.png', retina=True)
Out[28]:
A Taxa de PO\(_4^{-3}\) é estimado usando um ajuste linear entre Fosfato e idade.
A idade é derivada de dados de radio carbono. Isso porque o método do CFC é limitado a idades entre 5-35 anos, logo idades por 14C são recomendadas para massas d'água mais velhas (profundas).
A inclinação da reta é de 0.00139, o que significa que 1,39 \(\times\) 10\(^{-3} \mu\)mol de P é formado por ano.
In [29]:
Image('./figures/pvage.png', retina=True)
Out[29]:
Para maiores informações chequem o paper original no link:
http://dx.doi.org/10.1016/S0967-0637(99)00025-4
In [30]:
columns = [r'Temperature ($^\circ$C)', 'Salinity (PSU)', 'Silicate', r'NO ($\mu$mol kg$^{-1}$)']
index = ['AAIW', 'UCDW', 'UNADW', 'TDDW', 'AABW', 'Uncertainties']
data = np.c_[
[ 3.65, 2.53, 4, 2.3, 0.14, 0.16],
[34.18, 34.65, 34.97, 34.91, 34.68, 0.016],
[ 15, 67, 18, 36, 113, 1.9],
[ 506, 468, 422, 442, 526, 8.8]
]
df = DataFrame(data, index=index, columns=columns)
df.index.name = 'Water type'
df
Out[30]:
- AAIW: Água Antártica Intermediária (AIA) [600--800 metros]
- UCDW: Água Circumploar Superior (ACS) [1000-1200 metros]
- UNADW: Água Profunda do Atlântico Norte "Superior" [1000-2000 metros]
- TDDW: Água Profunda de Dois Graus (Mistura entre APAN e AFA no Atlântico Oeste)
- AABW: Água Profuda Antártica (AFA)
In [31]:
%%latex
\begin{align*}
x_1 T_1 + x_2 T_2 + x_3 T_3 + x_4 T_4 &= T \\
x_1 S_1 + x_2 S_2 + x_3 S_3 + x_4 S_4 &= S \\
x_1 O_1 + x_2 O_2 + x_3 O_3 + x_4 O_4 + -\alpha J_o &= O \\
x_1 N_1 + x_2 N_2 + x_3 N_3 + x_4 N_4 + R_N\alpha J_o &= N \\
x_1 P_1 + x_2 P_2 + x_3 P_3 + x_4 P_4 + R_P\alpha J_o &= P \\
x_1 + x_2 + x_3 + x_4 &= 1
\end{align*}
- \(O\), \(N\), \(P\) representam Oxigênio, Nitrato e fosfato inorgânicos;
- \(\alpha\) e \(J_o\) representam respectivamente uma idade efetiva da massa de água e consumo de oxigênio in-situ;
- \(R_N\) e \(R_P\) São as taxa estequiométrica de Redfiled para nitrato e fosfato.