Segundo parcial - Instrumentación geofísica
Miguel Cabrera
2024-10-15
Segundo parcial - Instrumentación Geofísica
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors#install.packages("zoo")
library(TSstudio)
library(moments)
library(signal)##
## Adjuntando el paquete: 'signal'
##
## The following object is masked from 'package:dplyr':
##
## filter
##
## The following objects are masked from 'package:stats':
##
## filter, polyEl correcto procesamiento de los datos implica que estén en un formato tabular. Este procesamiento se realizó por medio de Python al acceder a la información y convertirla a un data.frame. Posteriormente, esta se convirtió a csv para ser leída en el entorno de R.
Según el Copernicus Marine Service (2024), el código para realizar este procedimiento se presenta a continuación como línea de código comentada.
#import xarray as xr
#DS = xr.open_dataset("input_filename.nc")
#DS.to_dataframe().to_csv("output_filename.csv")
#https://help.marine.copernicus.eu/en/articles/4839794-how-to-convert-netcdf-to-csvA pesar de que logar esto en R también es posible, el proceso es mucho más largo y enrevesado. Aunque no es recomendable, se optó por cambiar de entorno para el cargue de datos dadas las facilidades anteriormente expuesta. Al haber realizado esto se cargan los datos.
Manejo de datos y uso de ERA5
Gráficas de SST de promedios mensuales
Para todos los años presentes en la serie de tiempo de la temperatura superficial se busca obtener los promedios mensuales y gráficarlo. Para realizar este procedimiento, en primer lugar, se cargan los datos en un formato tabular, los cuales se extrajeron como se explicó anteriormente.
Temperatura =
read_csv("Temperatura.csv", show_col_types = FALSE) %>%
select(valid_time, sst)
head(Temperatura) ## # A tibble: 6 × 2
## valid_time sst
## <dttm> <dbl>
## 1 1990-01-01 18:00:00 300.
## 2 1990-01-02 18:00:00 300.
## 3 1990-01-03 18:00:00 300.
## 4 1990-01-04 18:00:00 300.
## 5 1990-01-05 18:00:00 300.
## 6 1990-01-06 18:00:00 300.Como se puede apreciar, la variable Temperatura, correspondiente a los datos del SST presenta una frecuencia diaría (365.25). Ya que se desea obtener los promedios mensuales, se sobre escribe la variable y a partir de las fechas en un formato POSIXct se extrae el año, mes y día usando el paquete lubridate.
Temperatura = Temperatura %>%
mutate(
año = year(valid_time), # Extraer el año
meses = month(valid_time, label = TRUE, abbr = TRUE), # Extraer el mes
día = day(valid_time) # Extraer el día
) %>%
select(sst, año, meses, día)
head(Temperatura)## # A tibble: 6 × 4
## sst año meses día
## <dbl> <dbl> <ord> <int>
## 1 300. 1990 ene 1
## 2 300. 1990 ene 2
## 3 300. 1990 ene 3
## 4 300. 1990 ene 4
## 5 300. 1990 ene 5
## 6 300. 1990 ene 6A pesar de que deliveradamente se optó por eleminar la información de la hora en la que se tomó el dato de sst, ahora se tiene una separación útil de las fechas. Por lo tanto, se agruparán los años y meses para sacar el promedio a partir de estos criteríos aplicando una función mena() a todo el data.frame. Posteriormente, estos datos se convierten en un ts object (objeto de serie de tiempo) para emplear el paquete TSstudio con el fin de graficarlo.
Temperatura_promedio = Temperatura %>%
group_by(año, meses) %>%
summarise(sst = mean(sst))## `summarise()` has grouped output by 'año'. You can override using the `.groups`
## argument.Temperatura_promedio_ts = ts(data= Temperatura_promedio$sst,
end = c(2009, 12),
frequency = 12
)
ts_plot(Temperatura_promedio_ts)Gráficas de los promedios climatológcios
Un promedio climatológico es una medida estadística que representa las condiciones de, en este caso, una variable oceanográfica a lo largo de una serie de tiempo multianual. En particular, se pretende observar la variación mensual de la temperatura superficial del mar a partir de los promedios de todos los años de la serie de tiempo.
Para cumplir este cometido, se optó por reializar una box plot debido a que además de mostrar la media también indica los cuartiles y los valores atípicos de manera que se podrá realizar un análisis más profundo.
Temperatura %>%
ggplot(aes(x = factor(meses), y = sst)) +
geom_boxplot()+
theme_classic()
Distribución empírica de probabilidad y probabilidad acumulada
Antes de desarrollar la distribución empírica de probabilidad y probabilidad acumulada, se considera importante hacer un análisis de la forma de la distribución de la temperatura superficial de manera cuantitativa. Para ello se emplearán las medidas de curtosis y asimetría. Así, se conocerá la homogenidad y distribución de los datos de temperatura superficial entre los 20 años comprendidos entre 1989 y 2009.
La asimetría es la medida que indica la simetría de la distribución de una variable respecto a la media aritmética. Por otro lado, la curtosis indica la cantidad de datos que hay alrededor de la media. Eventualmente esto se verá gráficamente al realizar una gráfica de densidad de los datos.
#calcular la asimetría
skewness(Temperatura$sst)## [1] -0.01995451#calcular la curtosis
kurtosis(Temperatura$sst)## [1] 2.380878La asimetría es dec -0.0199; muy cerca de 0. Esto quiere decir que la distribución de los datos es prácticamente simétrica. Lo cual, implica que los valores de temperaturas están distribuidos de manera bastante equilibrada en torno a la media. No obstante, la curtosis es de 2.380, lo cual indica que la distribución es platicúrtica. Es decir, la distribución de las temperaturas tiene una forma más plana en el centro comparada con una distribución normal. Esto sugiere que los valores anómalos o extremos son menos frecuentes de lo que se esperaría.
Todas estas apreciaciones se pueden notar también en la distribución empiríca de probabilidad que se presenta acontinuación.
Temperatura %>%
ggplot(aes(x = sst)) +
geom_histogram(fill = "cornsilk", colour = "grey60", aes(y = ..density..)) +
geom_density()+
theme_classic()## Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(density)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
A patir de la gráfica de probabilidad acumulada se puede conocer la probabilidad que existe de que se presente un valor de temperatura superficial.
Temperatura %>%
ggplot(aes(x = sst)) +
stat_ecdf(geom = "step", color = "blue", size = 1) +
theme_classic()## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
En este gráfico se puede observar que la probabildidad es análoga a la información de los cuartiles. Es decir, el 50% de probalidad corresponde al 301.1, el 25% es de 300.5, y el 75% de 301.7.
summary(Temperatura$sst)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 299.2 300.5 301.1 301.1 301.7 303.4Calculo de indice Degree Heating Week
El Degree Heating Week (DHW), se utiliza para el monitoreo de los eventos de blaqueo de coral en los arrecifes de coral. El DHW es una medida que se utiliza para cuantificar el estrés térmico que experimentan los arrefices de coral debido a temperaturas del agua elevadas y sostenidas.
El DHW se calcula sumando las cantidades de grados celsius por encima de la temperatura máxima promedio del mes más cádio histórico, durante un periodo de tres meses. Este valor se utiliza para evaluar el riesgo de blaqueo de colar, ya que los colare sson organismos sensibles a los cambios en la temperatura del agua.
En primera instancia se halla la temperatura máxima promedio del mes más cádio histórico.
Temperatura %>%
group_by(meses) %>%
summarise(media_sst = mean(sst, na.rm = TRUE)) %>%
summarise(MMM = max(media_sst)) %>% # máxima media mensual (MMM)
pull(MMM) ## [1] 302.0658Aunque se podría filtrar para escoger dos años únicamente, se decidió hacer el ejercicio con todos los datos de la serie de tiempo. Después de conocer cuál era la temperatura máxima promedio del mes más cálido histórico, se calcula el exceso de calor, etendido como la diferencia en grados si la temperatura supera el máximo de temperatura definido anteriormente, que para este caso es de 302. Después, se definen los trimestres del año. Finalmente, se suman las cantidades de grados celsius por encima de la temperatura máxima promedio del mes más cádio histórico, esto durante un periodo de tres meses.
DHW = Temperatura %>%
#filter(año == 1991) %>% # Filtrar el año que te interesa
mutate(Exceso_calor = ifelse(sst > 302, sst - 302, 0)) %>% # Calcular el exceso de calor
mutate(trimestre = case_when(
meses %in% c("ene", "feb", "mar") ~ 1,
meses %in% c("abr", "may", "jun") ~ 2,
meses %in% c("jul", "ago", "sep") ~ 3,
meses %in% c("oct", "nov", "dic") ~ 4
)) %>%
group_by(año, trimestre) %>% # Agrupar por año y trimestre
summarise(DHW = sum(Exceso_calor, na.rm = TRUE)) # Calcular la suma trimestral del exceso de calor## `summarise()` has grouped output by 'año'. You can override using the `.groups`
## argument.head(DHW)## # A tibble: 6 × 3
## # Groups: año [2]
## año trimestre DHW
## <dbl> <dbl> <dbl>
## 1 1990 1 0
## 2 1990 2 0
## 3 1990 3 0.702
## 4 1990 4 1.20
## 5 1991 1 0
## 6 1991 2 0Los datos se convierte en objeto de serie de tiempo con una frecuencia trimestral y posteriormente se gráfican.
DHW_ts = ts(data= DHW$DHW,
end = c(2009),
frequency = 4) #Frecuencia trimestral
ts_plot(DHW_ts)Es impresionante notar que el último trimestre de 1994 se tienen valores tan altos del indice en su último trimestre.
Análisis de tendencia, estacionalidad y anomalías
A partir de los datos de oleaje de ERA5 para la zona asignada a cada uno gnererar para una coordenada aleatoria:
Se cargan los datos de oleaje
Altura = read_csv("Altura.csv")## Rows: 7305 Columns: 6
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): expver
## dbl (4): latitude, longitude, number, swh
## dttm (1): valid_time
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.Altura = Altura %>%
mutate(
año = year(valid_time), # Extraer el año
meses = month(valid_time, label = TRUE, abbr = TRUE), # Extraer el mes
día = day(valid_time) # Extraer el día
) %>%
select(swh, año, meses, día)
head(Altura)## # A tibble: 6 × 4
## swh año meses día
## <dbl> <dbl> <ord> <int>
## 1 1.15 1995 ene 1
## 2 1.65 1995 ene 2
## 3 1.75 1995 ene 3
## 4 1.94 1995 ene 4
## 5 1.61 1995 ene 5
## 6 1.58 1995 ene 6Gráfias de Hs de promedios diarios
Para hacer una gráfica de serie de tiempo de los promedios diaríos, en primer lugar, se cargan los datos en un formato tabular, los cuales se extrajeron del formato .nc a partir del código de python que se presentó en un inicio.
Altura =
read_csv("Altura.csv", show_col_types = FALSE) %>%
select(valid_time, swh)Para hacer la gráfica de serie de tiempo se convierte a la variable Altura en un objeto ts, teniendo en cuenta que se trata de una frecuencia diaría.
Altura_ts = ts(data= Altura$swh,
end = c(2014, 12),
frequency = 365.25)
ts_plot(Altura_ts)Gráficas de Hs de promedios mensuales
Se realiza un proceso análogo al que se empleó en el punto 1.1. Primero se convierten los datos, después se promedian los meses, se convierte en objeto ts y se gráfica.
Altura = Altura %>%
mutate(
año = year(valid_time), # Extraer el año
meses = month(valid_time, label = TRUE, abbr = TRUE), # Extraer el mes
día = day(valid_time) # Extraer el día
) %>%
select(swh, año, meses, día)
Altura_promedio = Altura %>%
group_by(año, meses) %>%
summarise(swh = mean(swh))## `summarise()` has grouped output by 'año'. You can override using the `.groups`
## argument.Altura_promedio_ts = ts(data= Altura_promedio$swh,
end = c(2014, 12),
frequency = 12
)
ts_plot(Altura_promedio_ts)Aumento anual de HS a partir de las gráficas de promedios diaros y mensuales
Para estudiar el incremento anual se empleará un modelo lienal.
Primero se desarrolla un modelo lineal a partir de los datos, mensuales, anuales y diaríos.
modelo_altura = lm(swh ~ año + meses + día, data = Altura)
summary(modelo_altura)##
## Call:
## lm(formula = swh ~ año + meses + día, data = Altura)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.16261 -0.30252 -0.03494 0.26535 2.82108
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.9726178 1.7478739 1.129 0.25911
## año -0.0002814 0.0008720 -0.323 0.74694
## meses.L -0.3924543 0.0174439 -22.498 < 2e-16 ***
## meses.Q 0.4595565 0.0173601 26.472 < 2e-16 ***
## meses.C 0.2194386 0.0173329 12.660 < 2e-16 ***
## meses^4 0.4022120 0.0174037 23.111 < 2e-16 ***
## meses^5 0.1195221 0.0175252 6.820 9.83e-12 ***
## meses^6 -0.3876288 0.0176061 -22.017 < 2e-16 ***
## meses^7 -0.1765738 0.0174398 -10.125 < 2e-16 ***
## meses^8 0.1294398 0.0173843 7.446 1.07e-13 ***
## meses^9 0.0849917 0.0174352 4.875 1.11e-06 ***
## meses^10 -0.0498545 0.0173758 -2.869 0.00413 **
## meses^11 -0.0345361 0.0173683 -1.988 0.04680 *
## día -0.0011836 0.0005719 -2.070 0.03852 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4297 on 7291 degrees of freedom
## Multiple R-squared: 0.2655, Adjusted R-squared: 0.2642
## F-statistic: 202.7 on 13 and 7291 DF, p-value: < 2.2e-16Ahora se extrae el coeficiente corespondiente al año del modelo de Altura.
coef_anual_altura = coef(modelo_altura)["año"]
cat("Aumento anual estimado (modelo diario):", coef_anual_altura, "\n")## Aumento anual estimado (modelo diario): -0.000281369En segundo lugar, se desarrolla un modelo lineal a partir de los datos de promedios mensuales anuales.
modelo_altura_promedio = lm(swh ~ año + meses, data = Altura_promedio)
summary(modelo_altura_promedio)##
## Call:
## lm(formula = swh ~ año + meses, data = Altura_promedio)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.41650 -0.12286 -0.00765 0.11153 0.47683
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.0102841 4.1387186 0.486 0.62763
## año -0.0003094 0.0020647 -0.150 0.88100
## meses.L -0.3929672 0.0412425 -9.528 < 2e-16 ***
## meses.Q 0.4597498 0.0412425 11.147 < 2e-16 ***
## meses.C 0.2195506 0.0412425 5.323 2.44e-07 ***
## meses^4 0.4014767 0.0412425 9.735 < 2e-16 ***
## meses^5 0.1200145 0.0412425 2.910 0.00397 **
## meses^6 -0.3885118 0.0412425 -9.420 < 2e-16 ***
## meses^7 -0.1759573 0.0412425 -4.266 2.92e-05 ***
## meses^8 0.1287081 0.0412425 3.121 0.00204 **
## meses^9 0.0849576 0.0412425 2.060 0.04054 *
## meses^10 -0.0505963 0.0412425 -1.227 0.22117
## meses^11 -0.0341935 0.0412425 -0.829 0.40793
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1844 on 227 degrees of freedom
## Multiple R-squared: 0.6742, Adjusted R-squared: 0.657
## F-statistic: 39.15 on 12 and 227 DF, p-value: < 2.2e-16Ahora se extrae el coeficiente corespondiente al año del modelo de Altura.
coef_anual_altura_promedio<- coef(modelo_altura_promedio)["año"]
cat("Aumento anual estimado (modelo mensual):", coef_anual_altura_promedio, "\n")## Aumento anual estimado (modelo mensual): -0.0003094307Es correcto que el resultado obtenido de los modelos lineales muestra un ligero descenso en la altura significativa de las olas (valores negativos para los coeficientes del año), lo cual indica que, en promedio, la altura de las olas ha disminuido anualmente.
Sin embargo , una disminución anual de la altura de las olas (basada en los datos y el modelo lineal) no necesariamente implica una relación directa con el cambio climático. Básicamente, el resultado obtenido podría estar influenciado por factores locales o eventos aislados que afectan la altura de las olas en la región . Aunque el cambio climático afecta a los océanos, este análisis simple no proporciona evidencia suficiente para vincular directamente la tendencia a una causa climática.
Serie de Hs sin tendencia, promedios mensuales, y serie de anomalías
A continuación se presenta la gráfica de los de la serie de tiempo, la gráfica de tendencia, el componente estacional y la gráfica de series de tiempo de anomalías Hs.
ts_decompose(Altura_promedio_ts)Más que sólo representar los promedios mensuales de la serie de tiempo se condiera que es mejor observar también la representación del box plot del mismo, tal como se plateo en el apartado 1.2.
Altura %>%
ggplot(aes(x = factor(meses), y = swh)) +
geom_boxplot()+
theme_classic()
En el siguiente mapa de calor también se puede apreciar esa variabilidad
de la altura de la ola a nivel estacional.
ts_heatmap(Altura_promedio_ts)Series de fourier, transformada de Fourier y periodograma de Welch
A partir de los datos de altura signficante declarados en la parte 2.
Nota: Seee, esta es la parte que lamentablemente menos entendí de series de tiempo. Estuve revisando como resolver las preguntas, pero habían paquetes incompatibles y no entendía del todo bien como realizar los procedimientos o interpretaciones.Soy consciente de que esta será la parte más mhe del trabajo
Transformada rápida de Fourier y realizar la gráfica de periodograma de Welch
A partir de los datos de altura significativa de la ola, primero se calcula el algoritmo de la transformada rápida de Fourier, para ello se emplea Power spectral density estimation using Welch s method.
# Frecuencia de muestreo
sampling_frequency = 12
# Calcular el periodograma usando el método de Welch con spec.pgram()
welch_result = spec.pgram(Altura_promedio$swh, spans = c(5, 5), taper = 0.1, plot = FALSE)
# Extraer la frecuencia y la densidad espectral de potencia
frecuencias <- welch_result$freq * sampling_frequency
potencia <- welch_result$spec
welch_df <- data.frame(Frecuencia = frecuencias, Potencia = potencia)
ggplot(welch_df, aes(x = Frecuencia, y = Potencia)) +
geom_line(color = "blue") +
labs(title = "Periodograma de Welch", x = "Frecuencia (Hz)", y = "Densidad espectral de potencia") +
theme_minimal()
Periodograma de Welch, pero con los datos de anomalías de Hs.
En primer lugar, se extraen los datos de anomalía de Hs.
descomp = decompose(Altura_promedio_ts)
componente_aleatorio <-
data.frame(
Random = na.omit(descomp$random)
)A continuación, se realizá un proceso análogo al anterior, pero con la descomposición de la serie de tiempo en su componente aleatorío, se decir, de anomalías.
# Frecuencia de muestreo
sampling_frequency = 12
# Calcular el periodograma usando el método de Welch con spec.pgram()
welch_result = spec.pgram(componente_aleatorio$Random, spans = c(5, 5), taper = 0.1, plot = FALSE)
# Extraer la frecuencia y la densidad espectral de potencia
frecuencias = welch_result$freq * sampling_frequency
potencia = welch_result$spec
welch_df = data.frame(Frecuencia = frecuencias, Potencia = potencia)
ggplot(welch_df, aes(x = Frecuencia, y = Potencia)) +
geom_line(color = "blue") +
labs(title = "Periodograma de Welch", x = "Frecuencia (Hz)", y = "Densidad espectral de potencia") +
theme_minimal()
Serie de altura significante generada solo por las frencuencias correspondientes al ENSO
No supe como hacer el filtro pas-band (┬┬﹏┬┬)
¯_(ツ)_/¯