Gráficos Avanzados y series de tiempo [GEO1012]
MSc. José A. Lastra, Elizabeth Daille & Montsserrat Truyol
2021-06-10
Introducción
En este laboratorio práctico se ahondará a un nivel básico-intermedio la manipulación de datos junto con la preparación de material gráfico. Adicionalmente se introducirán conceptos elementales del trabajo con series de tiempo.
Para este laboratorio práctico utilizaremos nuevas librerías como son: ggExtra, oce, corrplot, xts y scales. Estas librerías nos permitirán extender las capacidades de análisis de datos base de R.
1.- Lectura de datos
Abra el programa RStudio, cree y configure un nuevo script siguiendo el modelo básico que aparece en el siguiente código:
# Correlación y análisis de series de tiempo
# José A. Lastra
# 13.02.2021
#-----------------------------------------------------
#librerias
library(tidyverse)
library(oce)
library(ggExtra)
library(corrplot)
library(xlsx)
library(npphen)
library(xts)
library(scales)
#-----------------------------------------------------
setwd('C:/lab_2/')En este laboratorio se utilizará la tabla de nombre Araucania2018_CTD_C1.xls que se encuentra en el directorio de trabajo en la carpeta tables. Los datos de esta tabla corresponden a un crucero oceanográfico realizado por el Instituto de Fomento Pesquero (IFOP) durante el año 2018 enmarcado en el “Programa de manejo de las Floraciones de algas nocivas y toxinas marinas en el Oceáno Pacífico desde Biobío a Aysén (I Etapa) 2018”. La información se compone de 3 transectos y un total de 8 estaciones de monitoreo, ubicadas a 2, 5 y 10 millas naúticas de la costa.
Esta información es de libre acceso y usted puede disponer de la totalidad de esta en el sitio web de IFOP, en la sección Buscador de Informes.
1.1.- Gráficos básicos
Siguiendo con la idea del laboratorio número 1, lo primero será generar gráficos oceanográficos clásicos pero empleando las funciones de la librería oce. Esta librería está completamente dedicada al análisis de datos oceanográficos, permitiendo el trabajo con información de CTD, información acústica, mareógrafos, etc.
Para poder emplear adecuadamente la información, primero debemos convertir la información a un tipo de dato más familiar para la librería.
#filtrando solo una estación
datos.subset <- CTD.df %>% filter(Estaciones == 'Nehuentue' & Transecto == '10 millas')Con nuestra estación seleccionada procederemos a convertir esa tabla a un ojeto de tipo CTD con el código que aparece a continuación.
#creando objeto ctd
estacion <- as.ctd(salinity = datos.subset$PSU,
temperature = datos.subset$Temp,
pressure = datos.subset$Prof)Ahora podemos graficar de forma simple, empleando la función plot(). Podemos ver argumentos adicionales para un objeto ctd usando help(“plot,ctd-method”).
¿Qué podemos ver en este gráfico?
Si disponemos de más información en nuestra tabla como por ejemplo el oxígeno disuelto, podemos incorporarla a nuestro objeto mediante el siguiente código.
#creando variable con datos
o2 <- datos.subset$oxi_diss
#agregando datos a objeto
estacion <- oceSetData(estacion, name = "dissolved oxygen", value = o2 ,
unit = list(unit = expression(ml/L), scale = ""))
plot(estacion, which = c('temperature','salinity','density','dissolved oxygen'), type='l')
1.2.- Gráficos avanzados
1.2.1.- Frecuencias
Una forma clásica de trabajar con datos es mostrar gráficamente la distribución de los datos, para ver agrupamientos en torno a ciertos valores o conocer algunas características adicionales a nivel de resumen.
Ahora podemos crear gráficos que proporcionan mucha más información y pueden ser de mayor aporte para realizar inferencias sobre los datos.
Un ejemplo son los gráficos de dispersión, que muestran el nivel de asociación (o no) entre variables. El siguiente ejemplo muestra como realizar un diagrama de dispersión empleando ggplot2.
#diagrama de dispersión base
ggplot(data = CTD.df, aes(x = Temp, y = sigma))+
geom_point(color = "blue") +
geom_smooth(method = "lm", formula = y~x, se = T, color = "red") +
ylab("Sigma-t") + xlab("Temperatura °C")
#relacione temperatura/salinidad y contenido de oxígeno
ggplot(data = CTD.df, aes(x = Temp, y = PSU, color = oxi_diss)) +
geom_point() +
ylab("PSU") + xlab("Temperatura °C") +
scale_color_viridis_c()
Finalmente un gráfico que permite combinar información consiste en un histograma de dos dimensiones. Un ejemplo concreto de esta aplicación es el uso de la temperatura y la salinidad del mar para la caracterización de masas de agua.
# "histograma" 2D
ggplot(CTD.df, aes(x=Temp, y=PSU) ) +
geom_density_2d_filled(contour_var = "count")
Para agregar los histogramas de las variables usaremos la función ggMarginal().
#histograma 2D con distribución por variable
g <- ggplot(CTD.df, aes(x=Temp, y=PSU)) +
geom_point()+
geom_density_2d_filled(alpha = 0.7,contour_var = "count")+
theme(legend.position="bottom")
#agregar histograma de variables
ggMarginal(g, type="histogram")
1.2.2.- Gráficos de transecto
La información oceanográfica puede ser también representada de a nivel de transectos como veremos a continuación.
#filtro de transecto
datos.seccion <- CTD.df %>% filter(Estaciones == 'Nehuentue')
#as section data
seccion <- as.section(salinity = datos.seccion$PSU,
temperature = datos.seccion$Temp,
pressure = datos.seccion$Prof,
longitude = datos.seccion$Longitu,
latitude = datos.seccion$Latitud,
station = datos.seccion$Codigo,
sectionId = 'Nehuentué')
#plot simple
plot(seccion)
De igual manera, podemos representar los datos a nivel puntual para cada estación, siguiendo el modelo de transecta. Este tipo de gráficos es útil considerando varias variables
plot(sectionGrid(seccion,trim = T,p = standardDepths(5)),
ylim = c(0,35),ztype = "image",xtype = 'longitude',
ytpye = 'depth',showStations = T,showStart = T,
zcol = oce.colorsPalette(100))
1.2.3.- Correlación
El cálculo de correlaciones o nivel de asociación entre variables numéricas es un ejercicio común en el análisis de información. De forma elemental, nos permite entender ciertas relaciones de dependencia y asignar o entender el nivel de significancia de dos procesos que podrían estar relacionados. La correlación simple entre dos variables puede ser obtenida de forma sencilla empleando la función cor().
#ejemplo simple
#filtrar estación de interés
est.2mn <- CTD.df %>% filter(Estaciones == 'Nehuentue' & Transecto == '2 millas')
#calculo de correlación
cor(est.2mn$Temp,
est.2mn$sigma,
use ='pairwise.complete.obs')## [1] -0.9900476Ahora, podemos realizar gráficos que permitan evidenciar de forma visual y cuantitativa los niveles de asociación entre variables. Para esto usaremos las funciones de la librería corrplot y tomaremos el objeto est.2mn.
# seleccion de columnas numericas de interes
m.est <- est.2mn[,10:14]
# calculo de matriz de correlacion
c.est <- cor(m.est)
# plot
corrplot(c.est)
Podemos revisar opciones adicionales empleando help(‘corrplot’); por ejemplo:
#aparecen coeficientes de correlacion y solo panel superior
corrplot(c.est,method = 'number',type = 'upper',number.cex=1.2)
#incorporando significancia
#calculo de valor p
p <- cor.mtest(c.est, conf.level = .99)
#graficando correlaciones significativas
corrplot(c.est,method="number", type="upper",
sig.level = 0.05,p.mat = p$p,insig = "blank")
2.- Series de tiempo
A nivel general, muchos de los procesos físicos, químicos o biológicos tienen una dimensión espacial y a su vez una dimensión temporal que puede ser cuantificada y analizada de diferentes maneras. El monitoreo a diferentes frecuencias terporales permite comprender de mejor manera los sistemas en estudio y las conexiones existentes por ejemplo entre océano y atmósfera.
Para esto emplearemos algunas series de tiempo obtenidas desde: (i) el Goddard Earth Science Research correspondientes a la extensión de hielo sobre el mar Antártico y (ii) serie promedio de temperatura superficial del mar (Sea Surface Temperature, en inglés) sobre el Mar de Bellingshausen y Amundsen obtenida desde NOAA.
Para esto trabajaremos con la librería xts, la que nos permitirá manipular, remuestrear y ajustar nuestras series de tiempo dentro del entorno de R. Para cargar las series emplee el siguiente código:
#leyendo tablas de datos
#Temperatura superficial
sst <- read.xlsx('tables/NOAA_Bellingshausen_Amundsen_Seas.xls',1)
#Extensión de hielo
sea_ice <- read.xlsx('tables/SH_IceExt_Monthly_1978-2017.xls',1)Podemos visualizar las series empleando las funciones de gráficos de base o empleando ggplot2.
#SST
ggplot(data = sst, aes(x = fechas, y = sst))+
geom_line()+
theme_bw()+
xlab('')+ ylab('SST °C')
## sea ice extent
ggplot(data = sea_ice, aes(x = fechas, y = BA_sea/1000000))+
geom_line()+
theme_bw()+
xlab('')+ ylab('Sea ice extent square km./10^6')
2.1.- Análisis convencional
Lo primero que haremos para trabajar de forma conjunta con ambas series es considerar un periodo de tiempo común.
#revisar inicios de ambas series de tiempo
#sea ice
min(sea_ice$fechas)
max(sea_ice$fechas)
# sst
min(sst$fechas)
max(sst$fechas)#filtrar periodo de tiempo comun
sst_filter <- sst %>% filter(fechas >= '1982-01-01' & fechas <= '2017-12-31')
seaIce_filter <- sea_ice %>% filter(fechas >= '1982-01-01' & fechas <= '2017-12-31')Ahora que nuestra información está filtrada, procederemos a crear objetos de tipo serie de tiempo empleando la librería xts.
2.1.1.- Técnicas de suavizado
Una forma de obtener una serie de tiempo acotada a una ventana temporal (anual, semi-anual, semanal, etc.), ya sea por un fenómeno de interés o para ajustar frecuencias de muestreo entre series de tiempo; es posible mediante el uso de funciones como apply.monthly, apply.weekly, apply.yearly, etc.
Tomemos como ejemplo la serie de SST que tiene una frecuencia de muestreo diaria y queremos que tenga la misma frecuencia que nuestra serie de Sea Ice extent que es mensual.
NOTA: Es importante destacar que estas series de tiempo empleadas para el laboratorio no tienen celdas sin información, por lo cual podemos emplear la función mean sin mayores preocupaciones. En el caso de trabajar con series que no están completamente rellenas puede crear una función propia y aplicarla como se muestra en el siguiente ejemplo:
##función propia que borra NA's antes del cálculo
promedio <- function(x){
mean(x, na.rm=T)
}
#aplicando nuestra función
sst.m1 <- apply.monthly(sst.ts,FUN = promedio)Si queremos guardar nuestra nueva serie de tipo mensual debemos crear un objeto y escribirlo en el formato de salida que sea más conveniente para nuestra línea de procesamiento.
#guardando serie mensual
fechas <- index(sst.m)
datos <- as.vector(sst.m)
#creación de data frame
tabla.save <- data.frame(fechas = fechas, sst = datos)
#guardar en disco
write_csv(x = tabla.save,file = 'sst_monthly_Bellingshausen_Amundsen_1982_2017.csv')
write.xlsx(x = tabla.save,file = 'sst_monthly_Bellingshausen_Amundsen_1982_2017.xlsx')Igualmente, existen otras librerías como forecast, tidyquant o TTR que disponen de más funciones para el análisis de series de tiempo.
2.1.2.- Descomposición
Por lo general, las series de tiempo en el océano tienden a ser la sumatoria de diversos procesos de diferentes escalas de variabilidad. Por esto por lo que se pueden ver series con grandes fluctuaciones “aleatorias”. Para esto se recurre a la descomposición de los componentes de la serie o un “suavizado de la señal” más sofisticado.
Sigamos empleando nuestras series de SST y Sea Ice Extent mensuales, creando un objeto que sea amigable para la función de descomposición.
# descomposición
#extraer valores sst
sst.data <- as.vector(sst.m)
#crear objeto ts
serie.sst <- ts(data = sst.data, start = c(1982,1),
end = c(2017,12),frequency = 12)
#extraer valores sea ice
ice.data <- as.vector(ice.ts)
#crear objeto ts
serie.ice <- ts(data = ice.data, start = c(1982,1),
end = c(2017,12),frequency = 12)Con los objetos empleados pasaremos a generar la descomposición de las series.
#en type= "additive" para estacionaria y "mult" para no estacionaria.
sst.d <- decompose(serie.sst,type = 'additive')
plot(sst.d)
#extrayendo componentes para graficar con ggplot
fechas <- index(sst.m)
#tabla
components.df <- data.frame(fechas = fechas,
sst = sst.d$x,
ice.area = ice.d$x,
sst.trend = sst.d$trend,
ice.trend = ice.d$trend,
sst.seas = sst.d$seasonal,
ice.seas = ice.d$seasonal)Creación del gráfico combinado:
#trend plot
#base del gráfico
g <- ggplot(components.df, aes(x = fechas))
#agregando tendencia de la temperatura
g <- g + geom_line(aes(y = sst.trend, color = 'SST Trend'))
#agregando tendencia del sea ice
g <- g + geom_line(aes(y = ice.trend/1000000, color = 'Sea Ice Extent Trend'))
#agregando eje secundario
g <- g + scale_y_continuous(sec.axis = sec_axis(~., name = "Sea Ice Extent (km2/10^6)"))
#configurando el eje del tiempo con scales
g <- g + scale_x_date(labels=date_format("%Y"), breaks=date_breaks(width="2 year"),
limits=c(as.Date('1982-1-1'), as.Date('2017-12-31')))
#dando unos toques finales
g <- g + scale_colour_manual(values = c("blue", "red"))+
labs(y = "Trend SST [°C]",x = "Time", colour = "") +
theme_bw()+
theme(legend.position="bottom")
#revisando el gráfico
g
2.2.- Non parametric phenology npphen
Una aproximación más reciente, desarrollada por el LabGRS de la PUCV implica el análisis de series de tiempo empleando métodos no paramétricos. Esto implica que se pueden trabajar series con comportamientos anómalas considerando las probabilidades asociadas a un comportamiento estándar anual, dado por la serie histórica de datos y que permite la obtención de umbrales de confidencia para establecer anomalías positivas o negativas en base a un período de referencia.
Ahora utilizaremos npphen para 1 vector numérico de nuestra serie diaria de SST original, para eso use la función Phen().
## [1] -1.579158317 -1.579158317 -1.579158317 -1.579158317 -1.579158317
## [6] -1.579158317 -1.579158317 -1.579158317 -1.579158317 -1.579158317
## [11] -1.579158317 -1.579158317 -1.579158317 -1.579158317 -1.579158317
## [16] -1.579158317 -1.579158317 -1.579158317 -1.579158317 -1.579158317
## [21] -1.579158317 -1.579158317 -1.579158317 -1.579158317 -1.579158317
## [26] -1.579158317 -1.569138277 -1.569138277 -1.569138277 -1.569138277
## [31] -1.569138277 -1.559118236 -1.559118236 -1.559118236 -1.559118236
## [36] -1.549098196 -1.549098196 -1.549098196 -1.549098196 -1.549098196
## [41] -1.539078156 -1.539078156 -1.539078156 -1.539078156 -1.539078156
## [46] -1.539078156 -1.539078156 -1.539078156 -1.549098196 -1.549098196
## [51] -1.559118236 -1.569138277 -1.569138277 -1.579158317 -1.579158317
## [56] -1.579158317 -1.579158317 -1.579158317 -1.579158317 -1.579158317
## [61] -1.579158317 -1.579158317 -1.579158317 -1.579158317 -1.579158317
## [66] -1.579158317 -1.579158317 -1.579158317 -1.579158317 -1.579158317
## [71] -1.579158317 -1.569138277 -1.569138277 -1.569138277 -1.569138277
## [76] -1.569138277 -1.569138277 -1.569138277 -1.569138277 -1.559118236
## [81] -1.559118236 -1.559118236 -1.559118236 -1.559118236 -1.559118236
## [86] -1.559118236 -1.549098196 -1.549098196 -1.559118236 -1.559118236
## [91] -1.559118236 -1.559118236 -1.569138277 -1.569138277 -1.569138277
## [96] -1.569138277 -1.569138277 -1.569138277 -1.569138277 -1.569138277
## [101] -1.569138277 -1.569138277 -1.569138277 -1.559118236 -1.559118236
## [106] -1.559118236 -1.549098196 -1.549098196 -1.549098196 -1.539078156
## [111] -1.539078156 -1.529058116 -1.529058116 -1.519038076 -1.509018036
## [116] -1.509018036 -1.498997996 -1.498997996 -1.488977956 -1.488977956
## [121] -1.478957916 -1.478957916 -1.478957916 -1.468937876 -1.468937876
## [126] -1.458917836 -1.458917836 -1.448897796 -1.448897796 -1.438877756
## [131] -1.438877756 -1.428857715 -1.318637275 -1.308617234 -1.308617234
## [136] -1.298597194 -1.288577154 -1.288577154 -1.278557114 -1.268537074
## [141] -1.268537074 -1.258517034 -1.248496994 -1.238476954 -1.238476954
## [146] -1.228456914 -1.218436874 -1.208416834 -1.198396794 -1.188376754
## [151] -1.188376754 -1.178356713 -1.168336673 -1.158316633 -1.158316633
## [156] -1.148296593 -1.138276553 -1.138276553 -1.128256513 -1.128256513
## [161] -1.118236473 -1.108216433 -1.108216433 -0.777555110 -0.777555110
## [166] -0.767535070 -0.767535070 -0.767535070 -0.757515030 -0.757515030
## [171] -0.747494990 -0.747494990 -0.747494990 -0.366733467 -0.356713427
## [176] -0.356713427 -0.346693387 -0.346693387 -0.336673347 -0.336673347
## [181] -0.326653307 -0.316633267 -0.306613226 -0.306613226 -0.286573146
## [186] -0.276553106 -0.266533066 -0.256513026 -0.246492986 -0.236472946
## [191] -0.226452906 -0.216432866 0.334669339 0.354709419 0.364729459
## [196] 0.364729459 0.374749499 0.384769539 0.394789579 0.404809619
## [201] 0.414829659 0.424849699 0.575150301 0.575150301 0.585170341
## [206] 0.595190381 0.595190381 0.605210421 0.605210421 0.615230461
## [211] 0.625250501 0.625250501 0.635270541 0.645290581 0.655310621
## [216] 0.665330661 0.675350701 0.685370741 0.695390782 0.705410822
## [221] 0.705410822 0.705410822 0.715430862 0.715430862 0.715430862
## [226] 0.715430862 0.715430862 0.715430862 0.705410822 0.705410822
## [231] 0.695390782 0.695390782 0.685370741 0.675350701 0.675350701
## [236] 0.665330661 0.655310621 0.655310621 0.655310621 0.655310621
## [241] 0.655310621 0.645290581 0.645290581 0.645290581 0.635270541
## [246] 0.615230461 0.214428858 0.214428858 0.204408818 0.204408818
## [251] 0.204408818 0.204408818 0.204408818 0.204408818 0.204408818
## [256] 0.494989980 0.494989980 -0.106212425 -0.106212425 -0.106212425
## [261] -0.106212425 -0.106212425 -0.106212425 -0.106212425 -0.096192385
## [266] -0.006012024 0.004008016 0.014028056 0.024048096 0.024048096
## [271] 0.024048096 0.034068136 0.034068136 0.034068136 0.024048096
## [276] 0.024048096 -0.346693387 -0.346693387 -0.607214429 -0.607214429
## [281] -0.597194389 -0.597194389 -0.587174349 -0.587174349 -0.587174349
## [286] -0.577154309 -0.577154309 -0.577154309 -0.577154309 -0.577154309
## [291] -0.567134269 -0.567134269 -0.567134269 -0.567134269 -0.567134269
## [296] -0.567134269 -0.577154309 -0.577154309 -0.587174349 -0.607214429
## [301] -0.617234469 -0.637274549 -0.637274549 -0.647294589 -0.657314629
## [306] -0.657314629 -0.657314629 -0.667334669 -0.667334669 -0.667334669
## [311] -0.667334669 -0.677354709 -0.677354709 -1.228456914 -1.228456914
## [316] -1.228456914 -1.238476954 -1.238476954 -1.238476954 -1.238476954
## [321] -1.238476954 -1.238476954 -0.937875752 -0.937875752 -0.937875752
## [326] -0.937875752 -0.937875752 -0.937875752 -0.927855711 -1.248496994
## [331] -1.248496994 -1.258517034 -1.418837675 -1.428857715 -1.438877756
## [336] -1.438877756 -1.448897796 -1.448897796 -1.458917836 -1.458917836
## [341] -1.468937876 -1.468937876 -1.468937876 -1.468937876 -1.478957916
## [346] -1.478957916 -1.478957916 -1.478957916 -1.478957916 -1.488977956
## [351] -1.488977956 -1.488977956 -1.488977956 -1.488977956 -1.498997996
## [356] -1.498997996 -1.498997996 -1.498997996 -1.509018036 -1.509018036
## [361] -1.519038076 -1.519038076 -1.529058116 -1.529058116 -1.529058116Por otro lado, podemos obtener también los intervalos de confianza para medir la variabilidad de nuestra información y su comportamiento de largo plazo.
PhenKplot(x = sst$sst, dates = sst$fechas,h = 2,nGS = 365,
rge = c(-2,3),xlab = 'HS Day',ylab = 'SST°C')
3.- Ejercicios
A continuación, dispone de una sección con algunos ejercicios prácticos para fortalecer lo visto en este laboratorio e incorporar análisis y aplicación de conceptos.
Fecha de entrega: 23 de Junio, 2021
- Seleccione una de las 6 series de tiempo disponibles en la carpeta Series_de_tiempo:
- Realice un gráfico de la serie de tiempo empleando ggplot y genere un pequeño reporte de lo observado
- Realice un suavizado promedio mensual y un suavizado promedio anual en su serie, grafique la información, genere un reporte y guarde sus nuevas series en una tabla.
- Utilizando el archivo CIMAR_FIORDO.xlsx
- Elija tres estaciones y calcule correlaciones para temperatura, salinidad y densidad
- Repita el ejercicio, pero ahora filtre la información de sus estaciones, de acuerdo con los siguientes rangos de profundidad: 1) 0-10 metros, 2) 10-20 metros y 3) > 20 metros. ¿Hay diferencias? ¿Que las provoca?
- Grafique las tres estaciones empleando la librería oce.