Este html es parte del material de apoyo que preparé en febrero de 2019 para dar unas clases en el Máster de Oceanografía. Se trata de una descripción simplificada de cómo realizar el análisis del porcentaje de masas de agua que componen una determinada muestra. En este link puedes encontrar una presentación en pdf y hay más materiales relacionados, incluído el archivo R-markdown que generó este mismo html en mi repositorio de github.

¿Qué necesitamos?

Datos de inicio:

Observaciones
Tabla de características tipo (Source Water Type, SWT)
Posibles figuras de mezcla

Pasos a seguir:

Carga los datos del ejemplo:

Todos los datos necesarios están en la carpeta /rawdata del repositorio CO2AO.

# Load data ---------------------------------------------------------------

data <- read.delim("C:/Users/MFontela/Nextcloud/ULPGC/OMPGT10/INPUT_botellas.txt", na.strings="-999")
SWT <- read_excel("C:/Users/MFontela/Nextcloud/ULPGC/OMPGT10/INPUT_SWT.xlsx")
MixingFigures <- read.delim("C:/Users/MFontela/Nextcloud/ULPGC/OMPGT10/INPUT_MixingFigures.txt", header=FALSE)
names(MixingFigures) <- c("wm1", "wm2", "wm3", "wm4") #rename table
MixingFigures <- mutate(MixingFigures, wm1=as.character(wm1),wm2=as.character(wm2),wm3=as.character(wm3),wm4=as.character(wm4))

# WEIGHTS FOR W_MATRIX (MATRIX OF WEIGHTS)
W<-data.frame(Mass=100, Tpot=20, Sal=10, SiO4=2,NO3=3, NO=3,NO30=3,PO4=2, PO=2,PO40=2,O2=2, O20=2)

# % REDFIELD RATIOS FOR EXTENDED OMP
Redfield<-data.frame(NO3=9.3, PO4=145, Si=37, O2=-1)

Ya tienes cargados los datos iniciales (puedes ver la tabla entera haciendo scroll hacia la derecha ——>):

Y también las propiedades tipo de las masas de agua (SWT) que, apoyándote en bibliografía has considerado:

También cargaste las posibles figuras de mezcla (“Mixing Figures”):

Diagrama T-S

Siempre es bonito ver los datos de botella con las SWT consideradas, ¿no crees? Además te ayuda a ver si escogiste bien las masas de agua tipo y las figuras de mezcla.

Añade más variables

El siguiente paso es calcular aquellos trazadores que no estén entre los datos iniciales: NO, PO, los nutrientes preformados…

# Add variables NO, PO, AOU, preformed...--------------------------------------

data<-data %>%
  mutate(NO=O2+(Redfield$NO3*NO3),
        PO=O2+(Redfield$PO4*PO4),
         AOU=(exp(-135.29996 + 157228.8/(Tpot+273.15) - 66371490/(Tpot+273.15)^2 + 12436780000/(Tpot+273.15)^3 - 862106100000/(Tpot+273.15)^4 -(Sal*(0.020573 - 12.142/(Tpot+273.15) + 2363.1/(Tpot+273.15)^2))))-O2, #la primera parte es la sw_satO2 del mítico SEAWATER toolbox del matlab. reference: "Weiss, R. F. 1970: "The solubility of nitrogen, oxygen and argon in water and seawater."Deep-Sea Research., 1970, Vol 17, pp721-735.
                 O20=O2+AOU,
        NO30=NO3-(AOU/Redfield$NO3),
        PO40=PO4-(AOU/Redfield$PO4))

SWT<-SWT %>%
  mutate(NO=O2+(Redfield$NO3*NO3),
         PO=O2+(Redfield$PO4*PO4),
         AOU=O20-O2,
         NO30=NO3-(AOU/Redfield$NO3),
         PO40=PO4-(AOU/Redfield$PO4))

Si vuelves a ver ahora las dos tablas, verás que hay nuevas columnas (cada vez más a la derecha, sí…):

Y también en la tabla de las SWT:

Escogemos una observación

Asegúrate que eliminaste los datos de superficie a 100 metros, que no entran en el análisis. En nuestro caso la menor profundidad que existe en los datos es ya 100.3m, por lo tanto podemos seguir.

Lo siguiente fue definir una muestra concreta, en este ejemplo vamos a caracterizar la composición de masas de agua de la muestra 200. La muestra seleccionada tiene estas propiedades:

	Tpot	Sal	SiO4	NO	PO
200	3.129	34.96	28.8	429.172	440.881

Vamos a comprobar que proporción de masas de agua componen esa muestra concreta. Empezamos encontrando cuál es la figura de mezcla que mejor explica la composición de esta muestra concreta a través de un OMP clasico (de ahí la ‘c’ de ‘cOMP’ y del nombre del resto de las variables:). Recuerda: clásico hace referencia a que únicamente considera variables conservativas: Temperatura potencial, salinidad, el macronutriente silicato (asumido conservativo) y los trazadores NO & PO.

Empezamos por la posibilidad de mezcla 1:

sample<-data[observation,] 
sample<- sample %>%
  select("Tpot", "Sal", "SiO4", "NO", "PO")
# preparas matriz de pesos
W_matrixSel<-W%>%select("Mass", "Tpot", "Sal", "SiO4", "NO", "PO")
W_matrixSel<-bind_rows(W_matrixSel,W_matrixSel,W_matrixSel,W_matrixSel)
#¿Será la mejor figura de mezcla la número 1?
FigurademezclaUNO<-MixingFigures[1,] #La primera fila
  selectedSWTc<-SWT %>%
  select("wm","Tpot", "Sal", "SiO4", "NO", "PO")%>%
    filter(wm %in% FigurademezclaUNO)%>%
    arrange(factor(wm, levels = unlist(FigurademezclaUNO)))%>% #Para que la tabla conserve el mismo orden que en el archivo MixingFigures original
    select(-"wm") #Deseleccionas las masas de agua por comodidad
  #Calcula el valor medio de esas masas de agua:
  meanselectedSWTc<-selectedSWTc%>% mutate_all(funs(mean)) 
  #y calculas la desviación estándar
  stdselectedSWTc<-selectedSWTc%>%  mutate_all(funs(sd)) 
  
  #Normalizas! #normalizas la Mixing figure y le añades una columna de unos (NOTA: normalizar un valor X es restarle la media y dividir el resultado por la std)
  normSWTc<-as.data.frame(c(1,(selectedSWTc-meanselectedSWTc)/stdselectedSWTc)) #y añades un uno (la masa, NO TOCAR)
  # y ponderas por peso:
  wnormSWTc<-normSWTc* W_matrixSel
  
  #Repites el proceso para la sample (normalizar y ponderar por peso)
  normsample<-as.data.frame(c(1,(sample-meanselectedSWTc[1,])/stdselectedSWTc[1,])) #y le pones un 1 delante (la masa)
  wnormsample<-normsample*W_matrixSel[1,] 
  
  #la función lsqnonneg admite como entrada una matriz y un vector, los preparo:
  #¿Cómo? convirtiéndolos desde lo que son ahora (data.frame) a lo que pide la función:
  wnormSWTc<-t(as.matrix.data.frame(wnormSWTc)) #preparo matriz
  wnormsample<-as.vector(as.matrix.data.frame(wnormsample)) #preparo vector
  lsqnonneg(wnormSWTc,wnormsample)

## $x
## [1] 0.000000 0.000000 0.000000 1.099624
## 
## $resid.norm
## [1] 1738.029

# Interpreta la salida

¡Tremendo residuo.normalizado! Parece que esta figura de mezcla no es la que mejor le va, ¿no crees? Lo cierto es que la comprobación de la figura 1 podías habértela ahorrado, ya que la temperatura potencial de la muestra era de unos 3 grados centígrados, lo que quiere decir que mejor hubieras empezado por una figura de mezcla de masas de agua profundas…

Que te toca ahora… pues volver a plantear este mismo fragmento del script con la siguiente figura de mezcla… o meterlo en un loop for que lo haga por ti… ¡ya verás como no ye tan difícil!

classic OMP (cOMP): selección de figura de mezcla

# SWT classic OMP ---------------------------------------------------------
# {'Tpot'  'Sal'  'SiO4'  'NO'  'PO'}
SWTc<-SWT %>%
  select("wm","Tpot", "Sal", "SiO4", "NO", "PO")

sample<-data[observation,] 
sample<- sample %>%
  select("Tpot", "Sal", "SiO4", "NO", "PO")
#-------
# preparas matriz de pesos
W_matrixSel<-W%>%select("Mass", "Tpot", "Sal", "SiO4", "NO", "PO")
W_matrixSel<-bind_rows(W_matrixSel,W_matrixSel,W_matrixSel,W_matrixSel)

cOMP<-list() #Generas una lista para guardar los datos
   for (nMF in 1:nrow(MixingFigures)) { #Aquí es dónde planteas el for loop para todas las mixing figures posibles
      selectedSWTc<-SWTc %>%
        filter(wm %in% MixingFigures[nMF,])%>%
        arrange(factor(wm, levels = unlist(MixingFigures[nMF,])))%>% #Para que la tabla conserve el mismo orden que en el archivo MixingFigures original
        select(-"wm") #Selecciona las masas de agua
            meanselectedSWTc<-selectedSWTc%>%
              mutate_all(funs(mean)) #Calcula el valor medio de esas masas de agua
            stdselectedSWTc<-selectedSWTc%>%
              mutate_all(funs(sd)) #y la desviación estándar
                  normSWTc<-as.data.frame(c(1,(selectedSWTc-meanselectedSWTc)/stdselectedSWTc))
                  #solución provisional
                  if (nMF==9) { #Porque hay 3 watermasses en vez de 4
                    wnormSWTc<-normSWTc* W_matrixSel[1:3,]
                  } else {
                    wnormSWTc<-normSWTc* W_matrixSel}#normalizas la Mixing figure y le añades una columna de unos (NOTA: normalizar un valor X es restarle la media y dividir el resultado por la std)
  #Después, normalizas la sample
  normsample<-as.data.frame(c(1,(sample-meanselectedSWTc[1,])/stdselectedSWTc[1,])) #y le pones un 1 delante (la masa)
  wnormsample<-normsample*W_matrixSel[1,] #normalizas la muestra según el peso
  #la función lsqnonneg necesita una matriz y un vector, los preparo
  wnormSWTc<-t(as.matrix.data.frame(wnormSWTc)) #preparo matriz
  wnormsample<-as.vector(as.matrix.data.frame(wnormsample)) #preparo vector
   cOMP[[nMF]]<-(lsqnonneg(wnormSWTc,wnormsample)) 
   } 
#en realidad sólo quieres encontrar el menor residuo, que será la mejor figura posible:
resid<-c()
for (i in 1:nrow(MixingFigures)) {
    resid[i]<-(cOMP[[i]]$resid.norm)
}
bestMF=which(resid==min(resid)) #bestMF será la mejor Mixing Figure posible

extended OMP (eOMP)

Una vez que con el classicOMP has decidido qué figura de mezcla es la mejor, que en este caso es la figura número 4 (la formada por las masas de agua AAIW, LSW, NADWu, MW) pasas a aplicar a esa figura en concreto el extended OMP. Muy similar pero con una columna/incógnita más:

sample<-data[observation,] 
sample<- sample %>%
  select("Tpot", "Sal", "SiO4", "NO3", "PO4","O2")
#-------
# preparas matriz de pesos
W_matrixSel<-W%>%select("Mass", "Tpot", "Sal", "SiO4", "NO3", "PO4", "O2")
W_matrixSel<-bind_rows(W_matrixSel,W_matrixSel,W_matrixSel,W_matrixSel,W_matrixSel)

eOMP<-list()
      selectedSWTe<-SWT%>%
        select("wm","Tpot", "Sal", "SiO4", "NO30", "PO40","O20") %>% 
        filter(wm %in% MixingFigures[bestMF,])%>% #te quedas solo con las masas de agua que están en la mejor mixing figure
        arrange(factor(wm, levels = unlist(MixingFigures[nMF,])))%>% #Para que la tabla conserve el mismo orden que en el archivo MixingFigures original
        select(-"wm") #Deselecciona las masas de agua (por comodidad de cara los cálculos)
            meanselectedSWTe<-selectedSWTe%>%
              mutate_all(funs(mean)) #Calcula el valor medio de esas masas de agua
            stdselectedSWTe<-selectedSWTe%>%
              mutate_all(funs(sd)) #y la desviación estándar
                  normSWTe<-as.data.frame(c(1,(selectedSWTe-meanselectedSWTe)/stdselectedSWTe))
   # -------- extended: with Redfield ratios!
  Redfieldratios<-c(0, 0,Redfield$Si, Redfield$NO3, Redfield$PO4, Redfield$O2)
  columntoadd<-data.frame(X1=0, 1/(Redfieldratios*stdselectedSWTe[1,])) #el primer cero es de la masa (columna de nombre X1 por el momento)
  columntoadd[2:3]=0 #Para que te quite los Inf que salen de hacer la inversa de 0.
      normSWTe<-bind_rows(normSWTe, columntoadd) #Le añades la columna del extended
                                     wnormSWTe<-normSWTe* W_matrixSel#normalizas la Mixing figure y le añades una columna de unos
  #Después, normalizas la sample
  normsample<-as.data.frame(c(1,(sample-meanselectedSWTe[1,])/stdselectedSWTe[1,])) #y le pones un 1 delante (la masa)
  wnormsample<-normsample*W_matrixSel[1,] #normalizas la muestra según el peso


  #la función lsqnonneg necesita una matriz y un vector, los preparo
  wnormSWTe<-t(as.matrix.data.frame(wnormSWTe)) #preparo matriz
  wnormsample<-as.vector(as.matrix.data.frame(wnormsample)) #preparo vector
   eOMP<-(lsqnonneg(wnormSWTe,wnormsample)) 
   #Ponerle los nombres
     if (bestMF==9) { #Porque hay 3 watermasses en vez de 4 en la mixingFigure 9
    colnames<-c(paste(MixingFigures[bestMF,1:3]), "bio")
  } else {
    colnames<-c(paste(MixingFigures[bestMF,]), "bio")}
   names(eOMP$x)<-colnames
   eOMP #print the result

## $x
##        AAIW         LSW       NADWu          MW         bio 
##  0.01697792  0.02773752  0.84056464  0.11471992 29.57026700 
## 
## $resid.norm
## [1] 0.001115317

Resultado final

La composición de la muestra 200, explicada por la figura de mezcla AAIW, LSW, NADWu, MW, es:

1.7% de AAIW
2.77% de LSW
84.06% de NADWu
11.47% de MW

Ya sabes hacerlo para una muestra, así que ahora que ya eres un experto en R, que será mejor, ¿hacerlo muestra a muestra o montar un for loop para todas, cómo hiciste para elegir la figura de mezcla? ;) Dale caña!!

Distribución de masas de agua: análisis OMP extendido - Caso práctico

Marcos Fontela, febrero 2019