Tutorial (oceanografía): paquete oce
Pablo Jorgensen, basado en Kelley, D. E. (2018). Oceanographic Analysis with R. Springer
2023-05-28
Paquetes “oce” y “ocedata” para manejo de datos oceanográficos
El paquete oce es útil para quienes trabajan con datos oceanográficos. oce proporciona muchas funciones para:
leer datos oceanográficos,
procesar los datos de un instrumento de medición específico, y
visualizar los resultados siguiendo las convenciones oceanográficas (utilizando los gráficos base).
La función clave para importar datos en R es
?read.oce(), que reconoce automáticamente el tipo de
archivo. Si el reconocimiento no funciona, pruebe las funciones
específicas para el tipo de datos a leer (por ejemplo,
read.ctd() para archivos medidos con CTD,
read.odf() para archivos ODF (Ocean Data Files),
read.netcdf() para archivos NetCDF).
Para ver todas las funciones proporcionadas vean el manual de referencia
o ?oce, y luego hagan clic en el índice en la parte
inferior de la página. Les mostraré algunos gráficos que puede crear
para conjuntos de datos o propósitos específicos. La gran ventaja de
este paquete es que no es necesario transformar los datos en el típico
formato de data.frame para su visualización, y que puede manejar una
amplia variedad de formatos generados por diferentes instrumentos
oceanográficos.
El paquete también viene con algunos conjuntos de datos incluidos, por ejemplo:
ctd;
wind;
section (una sección hidrográfica);
adp (conjunto de datos del ADP (perfilador acústico-doppler));
echosounder echosounder (una submuestra atenuada de las mediciones realizadas con una ecosonda científica de Biosonics, en el marco del Experimento de Ondas Internas de Saint Lawrence (SLEIWEX));
landsat (una muestra del conjunto de datos del satélite Landsat);
argo (contiene datos de la boya argo ARGO 6900388 en el Atlántico Norte. Las boyas argo son robots que flotan y recogen información a diferentes profundidades en el mar);
lobo ( muestra de datos de lobo obtenida en el Brazo Noroeste de Halifax por Satlantic; lobo es un sistema completo de vigilancia de la calidad del agua que responde a la necesidad de realizar mediciones rutinarias, sólidas y precisas);
lisst (datos de un instrumento LISST, que es un instrumento óptico de difracción láser para medir el tamaño de las partículas y la concentración de los sedimentos en suspensión in situ).
La mayoría de los demás conjuntos de datos se han incorporado al
paquete de datos ocedata.
El paquete ocedata contiene un conjunto de datos
topográficos de 2 grados de resolución llamado topo2,
que da una matriz de elevación en metros sobre el nivel medio del mar.
La matriz puede ser contorneada con contour(topo2) que
produce un diagrama útil con intervalos de contorno auto-seleccionados y
ejes que van de 0 a 1, pero lo siguiente produce un resultado más
informativo (figura 1); noten el uso de los ejes x e y para evitar la
distracción del espacio alrededor del gráfico.
# install.packages("oce")
# install.packages("ocedata")
# install.packages("ncdf4")
library(oce)
library(ocedata)
data(topo2, package="ocedata")
lon <- seq(-179.5, 178.5, length.out=180) # see ?topo2
lat <- seq(-89.5, 88.5, length.out=90)
contour(lon, lat, topo2,
drawlabels=FALSE,
levels=c(0,-5000),
lty=c(1, 3), xaxs="i", yaxs="i", asp=1)
Figura 1. Contornos de elevación en metros sobre el nivel medio del mar
La función de R image() produce imágenes básicas, e imagep() agrega una barra de referencia de color (figura 2).
# par(mfrow=c(1, 2))
image(topo2)
imagep(topo2)
Figura 2. Mapas topográficos en metros respecto del nivel medio del mar producidos con image() (izquierda), e imagep() (derecha)
El paquete rggobi proporciona gráficos 3D interactivos, mientras que el paquete lattice proporciona gráficos estáticos. Estos gráficos son los más adecuados para funciones de suavizado, por lo que pondremos un ejemplo utilizando interpBarnes() para suavizar el conjunto de datos del viento con un diagrama de contorno, y con diagrama 3D de malla (figura 3).
library(lattice)
data(wind, package="oce")
g <- interpBarnes(wind$x, wind$y, wind$z)
# par(mfrow=c(1, 2))
contour(g$xg, g$yg, g$zg, xlab="x", ylab="y", labcex=1)
# Diagrama 3D con wireframe()
wireframe(g$zg, xlab="x", ylab="y", zlab="z", cex=5)
Figura 3. Datos tridimensionales representados con contornos con la función contour() del lattice (izquierda), y una malla reticular con la función wireframe (derecha)
La clase ts maneja datos de series temporales que se muestrean a una frecuencia constante. Por ejemplo, el conjunto de datos giss del paquete ocedata puede convertirse en un objeto ts y trazarse como en el panel de la izquierda de la figura 4. Dado que la señal mensual de giss es bastante irregular, puede querer construirse un gráfico que integre observaciones a lo largo de años o décadas con un gráfico de caja (“box-plot”). La función round() puede utilizarse para redondear a potencias de diez (figura 4, derecha).
# par(mfrow=c(1, 2))
data(giss, package="ocedata")
Ta <- ts(giss$index, start=giss$year[1], frequency=1 / mean(diff(giss$year)))
plot(Ta, ylab="Anomalía de la temperatura")
decade <- round(time(Ta), -1) # yields an indicator of decade
boxplot(Ta ~ decade, notch=TRUE, xlab="Tiempo", ylab="Anomalía de la temperatura")
Figura 4. Anomalía de la temperatura superficial mostrada con series temporales y gráficos de caja
Propiedades del agua de mar
En oceanografía se utilizan símbolos griegos para diferentes aspectos de la densidad, en particular \(\sigma_t\) para una medida de la densidad del agua de mar si se lleva a la superficie sin intercambio de calor ni cambio de salinidad, sólo descompresión: \(\sigma_t = 1000 (\rho - 1)\), donde \(\rho\) es la densidad real, normalmente un número como 1,02437 para el que \(\sigma_t = 24,37\). Sin embargo, el agua no es incompresible. A presiones oceánicas profundas se contrae sustancialmente. Entonces, a menudo se utiliza un refinamiento adicional, \(\sigma_\theta\), que tiene en cuenta el enfriamiento adiabático por expansión (expandida, tendrá la “temperatura potencial”). Para el cálculo preciso de la densidad a partir de los datos de conductividad (“C”, una medida de “S”), temperatura (“T”) y presión (“D”, porque la profundidad es proporcional a la presión, de ahí “CTD”), es necesario utilizar funciones polinómicas empíricas con muchos términos. Afortunadamente, el paquete oce proporciona varias funciones que trabajan con las propiedades del agua de mar. Sus nombres comienzan con sw, y por lo tanto su documentación se puede encontrar con help.search(“sw”, package=“oce”). Las funciones sw que pueden aceptar objetos oce como primer argumento buscan dentro del objeto para encontrar toda la información necesaria, por ejemplo, swRho() utiliza los valores de longitud y latitud almacenados en un CTD para el cálculo de la densidad del agua de mar.
library(oce)
data(ctd, package="oce")
head(swRho(ctd, eos="gsw"), 3)
#> [1] 1022.233 1022.233 1022.235El paquete oce permite calcular la densidad del agua de mar en función de la presión, la salinidad y la temperatura (\(\rho = \rho(S,T,p)\)). Por ejemplo, ¿cuál es la densidad del agua de mar a una presión de 100 dbar, una salinidad de 34 PSU y una temperatura de 10 \(^{\circ}\)C?. Considerando que no se indican la longitud y la latitud, se utiliza la ecuación de estado de la UNESCO, que es especificada mediante el argumento eos para las diversas funciones sw:
library(oce)
swRho(34, 10, 100, eos="unesco")
#> [1] 1026.624Densidad del agua de mar
- Con los datos del ejemplo anterior, ¿qué temperatura tendría la parcela de agua si se eleva adiabáticamente hasta la superficie? (AYUDA usen la función swTheta()), sabiendo que la temperatura potencial es función de S, T y p (\(\theta = \theta(S,T,p)\))
RESPUESTA
library(oce)
swTheta(34, 10, 100, eos="unesco")
#> [1] 9.988599- ¿Qué densidad tendría si se eleva adiabáticamente hasta la superficie?
RESPUESTA
library(oce)
swRho(34, swTheta(34, 10, 100, eos="unesco"), 0, eos="unesco")
#> [1] 1026.173- ¿Qué densidad tendría si se baja unos 100 m, aumentando la presión a 200 dbar?
RESPUESTA
library(oce)
theta <- swTheta(34, 10, pressure=100, referencePressure=200, eos="unesco")
swRho(34, theta, 200, eos="unesco")
#> [1] 1027.074Diagramas T-S
Los diagramas T-S (temperatura potencial vs. salinidad) se utilizan para identificar las masas de agua, considerando que la temperatura potencial y la salinidad de una parcela de agua se conservan mientras permanezca aislada de la superficie (donde puede ganar o perder calor o agua dulce) y no se mezcle con otras masas de agua. Las masas de agua profundas conservan sus características T-S durante largos períodos de tiempo y pueden identificarse fácilmente en un gráfico T-S. La temperatura y la salinidad se combinan para determinar la densidad potencial del agua de mar; los contornos de densidad potencial constante se muestran a menudo en los diagramas T-S como el de la figura 5.
library(oce)
plotTS(as.ctd(c(30,40),c(-2,20),rep(0,2)), eos="unesco", type="n")
Figura 5. Diagrama T-S
Detección de la capa de mezcla
La capa superior del océano suele ser mezclada por olas o por la convección que se genera por la pérdida de calor en la superficie, provocando turbulencias. Esta zona superficial de agua revuelta en la que la temperatura se mantiene constante con la profundidad se denomina capa de mezcla. A veces se denomina capa de mezcla superficial para distinguirla de las capas homogéneas del interior o del fondo. En la base de la capa de mezcla, un gradiente de densidad, o picnoclina, separa el agua más ligera de la capa de mezcla del agua más densa que hay debajo. Si una partícula de densidad \(\rho_2\) de debajo de la picnoclina se desplaza por encima de la misma, estará rodeada por agua de densidad \(\rho_1\). La partícula más densa será más pesada que el agua circundante y volverá a caer hacia la interfase. Del mismo modo, una partícula de agua de la capa superior desplazada por debajo de la interfaz será más ligera que el agua circundante y se desplazará hacia la interfaz. La eficacia de la barrera aumenta con el incremento de la diferencia de densidad \(\rho_2-\rho_1\). La reducción del transporte vertical significa que las capas a ambos lados de la picnoclina están bastante aisladas unas de otras. Esta separación tiene importantes consecuencias biológicas: el fitoplancton crece sobre todo en la capa superior, pero el suministro de nutrientes de la capa inferior a la superior está, hasta cierto punto, bloqueado por la picnoclina. En aguas estratificadas se observa un apilamiento vertical de la columna de agua del océano. La estabilidad del apilamiento, las profundidades de las picnoclinas prominentes (niveles de mayor cambio de densidad y, por tanto, de estratificación más estable) y las fuerzas disponibles para impulsar el flujo ascendente (afloramiento) y la mezcla vertical, varían mucho en el océano mundial, lo que afecta al potencial de producción fotosintética de las distintas regiones.
La profundiad de la capa de mezcla puede calcularse a partir de las mediciones hechas con un CTD (“Conductivity, Temperature, Depth”), un instrumento que sirve para registrar las variaciones de la Conductividad, Temperatura y Profundidad del agua en un perfil de profundidad. Una función genérica del paquete oce permite examinar los datos de un CTD, por ejemplo:
library(oce)
data(ctd)
ctd
#> ctd object, from file '/Users/kelley/git/oce/create_data/ctd/ctd.cnv', has data as follows.
#> scan[1:181]: 130, 131, ..., 309, 310
#> timeS[1:181]: 129, 130, ..., 308, 309
#> pressure[1:181]: 1.480, 1.671, ..., 43.903, 44.141
#> depth[1:181]: 1.468, 1.657, ..., 43.542, 43.778
#> temperature[1:181]: 14.225, 14.230, ..., 2.9190, 2.9194
#> salinity[1:181]: 29.921, 29.921, ..., 31.395, 31.393
#> flag[1:181]: 0, 0, ..., 0, 0
summary(ctd)
#> CTD Summary
#> -----------
#>
#> * Instrument: SBE 25
#> * Temp. serial no.: 1140
#> * Cond. serial no.: 832
#> * File: "/Users/kelley/git/oce/create_data/ctd/ctd.cnv"
#> * Original file: c:\seasoft3\basin\bed0302.hex
#> * Start time: 2003-10-15 15:38:38
#> * Sample interval: 1 s
#> * Cruise: Halifax Harbour
#> * Vessel: Divcom3
#> * Station: Stn 2
#> * Mean Location: 44.684N 63.644W
#> * Data Overview
#>
#> Min. Mean Max. Dim. NAs OriginalName
#> scan 130 220 310 181 0 scan
#> timeS [s] 129 219 309 181 0 timeS
#> pressure [dbar] 1.48 22.885 44.141 181 0 pr
#> depth [m] 1.468 22.698 43.778 181 0 depS
#> temperature [°C, IPTS-68] 2.919 7.5063 14.237 181 0 t068
#> salinity [PSS-78] 29.916 31.219 31.498 181 0 sal00
#> flag 0 0 0 181 0 flag
#>
#> * Processing Log
#>
#> - 2018-11-14 20:03:47 UTC: `create 'ctd' object`
#> - 2018-11-14 20:03:47 UTC: `read.ctd.sbe(file = file, debug = 10, processingLog = processingLog)`
#> - 2018-11-14 20:03:47 UTC: `oce.edit(x = ctd, item = "startTime", value = as.POSIXct(gsub("1903", "2003", format(ctd[["startTime"]])), tz = "UTC") + 4 * 3600, reason = "file had year=1903, instead of 2003", person = "Dan Kelley")`El objeto de R no es un data.frame, pero si se utiliza una función de graficado no es necesario convertirlo. Al cargar el paquete oce pueden utilizar la función genérica plot(), ya que reconoce automáticamente el tipo de datos y aplica el método de graficado correspondiente (figuras 6 y 7). Los argumentos adicionales permiten la configuración, con un mayor refinamiento disponible con plotTS() y plotProfile().
plot(ctd)
Figura 6. Perfil de temperatura y salinidad con un CTD con la funciión de graficado base. Noten que, por defecto, se devuelve una figura con cuatro paneles
plot(ctd, which= 1)
Figura 7. Ídem a la figura anterior, pero en esta figura indiqué que se grafique solamente el perfil de temperatura con el argumento which()
El método de graficado para los datos CTD genera automáticamente varios paneles (figuras 6 y 7). Si desea graficar parte de los paneles, recordar usar el argumento which (figura 7). En el argumento which pueden especificar un escalar o un vector de gráficos deseados. De hecho, todos los gráficos que el paquete puede producir tienen un número del 1 al 33, pero no todos pueden ser generados con el mismo tipo de datos. Si el tipo de datos de entrada es CTD, se producen los gráficos número 1,2,3 y 5 (ver también ?plot,ctd-method).
El paquete oce proporciona funciones de lectura para una amplia variedad de instrumentos oceanográficos, por lo que no será necesaria ninguna codificación especial (a menos que sea un instrumento inusual). Para la mayoría de los tipos de datos, los datos leídos por algún otro medio también pueden transformarse en un objeto de este tipo con la función as.ctd(S, T, p), dada la salinidad, la temperatura y la presión; los argumentos adicionales permiten la especificación de la longitud y la latitud.
Ejemplo de los datos de una ecosonda
La ecosonda es un tipo de sonar utilizado para determinar la profundidad de la columna de agua mediante la transmisión de ondas sonoras hacia el fondo. Pero la ecosonda también puede referirse a las “ecosondas” hidroacústicas definidas como sonido activo en el agua (sonar) utilizadas para estudiar, por ejemplo, los peces (comúnmente aplicadas en la pesca) (figura 8). Los datos de una ecosonda pueden estar en el siguiente formato.
library(oce)
data(echosounder)
echosounder
#> echosounder object, from file '/Users/kelley/Dropbox/data/archive/sleiwex/2008/fielddata/2008-07-01/Merlu/Biosonics/20080701_163942.dt4', has data as follows.
#> timeSlow[1:190]: 2008-07-01 16:39:41, 2008-07-01 16:39:42, ..., 2008-07-01 16:42:49, 2008-07-01 16:42:50
#> latitudeSlow[1:190]: 47.879, 47.879, ..., 47.882, 47.882
#> longitudeSlow[1:190]: -69.724, -69.724, ..., -69.729, -69.729
#> depth[1:54]: 39.640, 38.924, ..., 2.4412, 1.7258
#> time[1:389]: 2008-07-01 16:39:41, 2008-07-01 16:39:41, ..., 2008-07-01 16:42:50, 2008-07-01 16:42:50
#> latitude[1:389]: 47.879, 47.879, ..., 47.882, 47.882
#> longitude[1:389]: -69.724, -69.724, ..., -69.729, -69.729
#> a, a 389x54 array with value 131.375 at [1,1] position
summary(echosounder)
#> Echosounder Summary
#> -------------------
#>
#> * File source: "/Users/kelley/Dropbox/data/archive/sleiwex/2008/fielddata/2008-07-01/Merlu/Biosonics/20080701_163942.dt4"
#> * Transducer serial #: TX03022C
#> * File type: DT4 pre v2.3
#> * Beam type: single-beam
#> * Channel: 1
#> * Measurements: 2008-07-01 16:39:41 UTC to 2008-07-01 16:42:50 UTC
#> * Sound speed: 1490.32 m/s
#> * Pulse duration: 0.0004 s
#> * Frequency: 418000.000000
#> * Blanked samples: 57
#> * Pings in file: 780
#> * Samples per ping: 3399
#> * Time: 2008-07-01 16:39:41 to 2008-07-01 16:42:50 (389 samples, mean increment 0.488866 s)
#> * Data Overview
#>
#> Min. Mean Max. Dim. NAs
#> timeSlow 1214930381 1214930476 1214930570 190 0
#> latitudeSlow [°N] 47.879 47.881 47.882 190 0
#> longitudeSlow [°E] -69.729 -69.726 -69.724 190 0
#> depth [m] 1.7258 20.683 39.64 54 0
#> time 2008-07-01 16:39:41 2008-07-01 16:41:15 2008-07-01 16:42:50 389 0
#> latitude [°N] 47.879 47.881 47.882 389 0
#> longitude [°E] -69.729 -69.726 -69.724 389 0
#> a 11.325 2492.3 364234 389x54 0
#>
#> * Processing Log
#>
#> - 2019-12-18 13:18:39 UTC: `read.oce("~/Dropbox/data/archive/sleiwex/2008/fielddata/2008-07-01/Merlu/Biosonics/20080701_163942.dt4")`
#> - 2019-12-18 13:18:39 UTC: `read.echosounder("/Users/kelley/Dropbox/data/archive/sleiwex/2008/fielddata/2008-07-01/Merlu/Biosonics/20080701_163942.dt4", channel=1, soundSpeed=(missing), tz="UTC", debug=0, processingLog)`
#> - 2019-12-18 13:18:39 UTC: `subset.echosounder(x, subset=depth < 40)`
#> - 2019-12-18 13:18:41 UTC: `decimate(x = echosounder, by = c(2, 40))`
plot(echosounder, which = 2, drawTimeRange = TRUE, drawBottom = TRUE)
Figura 8. Profundidad de la columna de agua con una ecosonda
Graficado de datos Landsat
Los datos Landsat pueden descargarse del sitio web earthexplorer del USGS. El siguiente es un subconjunto de la imagen Landsat-8 designada LC80080292014065LGN00, una imagen de marzo de 2014 que cubre Nueva Escocia y partes de la Bahía de Fundy y la Plataforma de Escocia (figura 9).
library(oce)
data(landsat)
landsat
#> Landsat object, ID LC80080292014065LGN00
#> Data (bands or calculated):
#> "aerosol" has dimension c(79,80)
#> "blue" has dimension c(79,80)
#> "green" has dimension c(79,80)
#> "red" has dimension c(79,80)
#> "nir" has dimension c(79,80)
#> "swir1" has dimension c(79,80)
#> "swir2" has dimension c(79,80)
#> "panchromatic" has dimension c(158,160)
#> "cirrus" has dimension c(79,80)
#> "tirs1" has dimension c(79,80)
#> "tirs2" has dimension c(79,80)
plot(landsat, band = "temperature")
Figura 9. Ejemplo de imagen Landsat
Usando banda = “temperature” se trazará una estimación de la temperatura de brillo del satélite, calculada a partir de la banda tirs1. Otras opciones de banda para el Landsat-8 son: “aerosol”, “blue”, “green”, “red”, “nir”, “swir1”, “swir2”, “panchromatic”, “cirrus”, “tirs1”, or “tirs2”.
Graficado de líneas costeras
oce proporciona un conjunto de datos de la línea de costa descargado de Natural Earth (figura 10). Este conjunto de datos tiene una resolución gruesa con una escala de 1:110M y 10.696 puntos, lo que lo hace más adecuado para los gráficos a escala mundial trazados a un tamaño pequeño. Si quieren archivos de la línea de costa con una resolución más fina, utilicen el paquete ocedata.
library(oce)
data(coastlineWorld)
coastlineWorld
#> coastline object, from file '/Users/kelley/git/oce/create_data/coastline/ne_110m_admin_0_countries/ne_110m_admin_0_countries.shp', has data as follows.
#> longitude[1:10766]: NA, 61.211, ..., 30.323, 29.839
#> latitude[1:10766]: NA, 35.65, ..., -22.272, -22.102
plot(coastlineWorld)
Figura 10. Costas del mundo
Si quieren adicionar elementos a un mapa, como puntos, líneas, contornos o texto, pueden crear un mapa usando una proyección específica proporcionada por el paquete rgdal, con la función mapPlot() de oce. La Fig. 11 muestra, por ejemplo, la trayectoria del crucero Endeavour (proporcionada en el paquete ocedata) añadiendo la función mapPoints().
# options(timeout = 600) # para evitar que se frene la descarga.
# install.packages("sf")
library(sf)
# Carga los datos de la trayectoria.
data(coastlineWorld)
data(endeavour, package = "ocedata")
par(mar = rep(0.5, 4)) # esto elimina el margen de la figura.
mapPlot(coastlineWorld, type = 'l', col = 'gray', projection = "+proj=moll")
# Para agregar los puntos usen mapPoint().
mapPoints(endeavour$longitude, endeavour$latitude, pch = 20, col = 'red')
Figura 11. Viaje del Endeavour
Notas suplementarias
Es posible colorear los mapas topográficos (figura 12). Los datos necesarios están en el conjunto de datos topoWorld.
data(topoWorld)
topo <- decimate(topoWorld, 2) # coarsen grid: 4X faster plot
lon <- topo[["longitude"]]
lat <- topo[["latitude"]]
z <- topo[["z"]]
cm <- colormap(name = "gmt_globe")
drawPalette(colormap = cm)
par(mar = rep(0.5, 4))
mapPlot(coastlineWorld, type = 'l', col = 'gray', projection = "+proj=moll", grid = FALSE)
mapImage(lon, lat, z, colormap = cm)
Figura 12. Costas del mundo con detalle de la topografía
Cuando se hace un acercamiento a latitudes medias, se suele utilizar la proyección cónica conformacional de Lambert (figura 13).
par(mar = c(2, 2, 1, 1))
lonlim <- c(-80, 0)
latlim <- c(20, 60)
mapPlot(coastlineWorld, projection = "+proj=lcc +lat_1=30 +lat_2=50 +lon_0=-40",
col = "lightgray", longitudelim = lonlim, latitudelim = latlim)
Figura 13. Proyección cónica conformacional de Lambert
El paquete ocedata tiene un conjunto de datos llamado levitus, que almacena la temperatura de la superficie del mar (SST) y la salinidad (SSS) de la versión 2013 del World Ocean Atlas (WOA), comúnmente conocido como el atlas “Levitus”” (versión del paquete 0.1.3) (figura 14).
data(levitus, package = "ocedata")
par(mfrow = c(2,1))
imagep(levitus$longitude, levitus$latitude, levitus$SST, col = oceColorsJet, zlim = c(-2, 30))
imagep(levitus$longitude, levitus$latitude, levitus$SSS, col = oceColorsJet, zlim = c(20, 40))
Figura 14. Temperatura (arriba) y salinidad (debajo) superficial del mar en 2013