TutoRial: Impacto del Cambio Climático en la Temperatura del Agua y la Acuicultura en el Mediterráneo

Autor/a

Javier Atalah

Fecha de publicación

29 de noviembre de 2023

1 Cargar los paquetes

Installar con install.packages() si no los tienes.

library(tidyverse) # para gráficas y manipulación de datos
library(sf) # para analysisi espaciales
library(stars) # para leer fichero .nc
library(rnaturalearth) # para los mapas
theme_set(theme_minimal()) # definir el tema para las gráficas

2 Sitios de interés

En primer lugar, defino el área geográfica de interés para el estudio. En este caso, nos interesa las zonas donde se lleva a cabo la acuicultura de peces. Como ejemplo, he elegido el área de San Pedro del Pinatar en la Región de Murcia, donde existen granjas de lubina y dorada.

Este código en R está utilizando la librería sf para manejar datos espaciales.

  1. La variable bbox_s1 se crea para contener el resultado del cálculo del bounding box.

  2. st_bbox es una función que crea un objeto de tipo “bounding box” a partir de las coordenadas especificadas.

  3. Las coordenadas se proporcionan en forma de parámetros con nombres, como xmin, xmax, ymin y ymax, que representan los valores mínimos y máximos de longitud y latitud respectivamente.

  4. El argumento crs especifica el sistema de referencia de coordenadas (en este caso, 4326, que es el sistema de coordenadas geográficas de referencia WGS 84).

  5. st_as_sfc() se utiliza para convertir el bounding box en una geometría simple en 2D (SFC), que es la forma en que sf maneja las geometrías.

bbox_s1 <-
  st_bbox(c(
    xmin = -0.64,
    xmax = -0.51,
    ymin = 37.8,
    ymax = 37.9
  ),
    crs = st_crs(4326)) %>% 
  st_as_sfc()

bbox_s1
Geometry set for 1 feature 
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -0.64 ymin: 37.8 xmax: -0.51 ymax: 37.9
Geodetic CRS:  WGS 84

Este fragmento de código en R utiliza la librería ggplot2 y sf para crear un gráfico:

  1. Se crea un objeto spain que representa el mapa de España, utilizando la función ne_countries de la librería rnaturalearth.

  2. ggplot() inicializa el gráfico.

  3. geom_sf() se utiliza para agregar capas al gráfico, en este caso, se están añadiendo dos capas de objetos de la clase sf para representar el bounding box y el mapa de España.

  4. En la primera capa, bbox_s1 se representa como un polígono con relleno de color rojo oscuro y transparencia del 30%.

  5. En la segunda capa, se representa el mapa de España.

  6. coord_sf() se utiliza para establecer los límites de los ejes x e y del gráfico. Aquí, se establecen los límites de longitud entre -1 y -0.4, y los límites de latitud entre 37.5 y 38.

spain <- ne_countries(scale = "large", returnclass = "sf", country = 'spain')

ggplot() + 
  geom_sf(data = spain)

last_plot() + 
  geom_sf(data = bbox_s1, fill = "darkred",alpha = .3) + 
  coord_sf(xlim =  c(-1, -.4), ylim = c(37.5, 38))

3 Datos Copernicus

En primer lugar, se recopilan los datos de temperatura entre 0 y 40 metros de profundidad para las áreas de interés. En este ejemplo, he seleccionado un área cerca de la costa de San Pedro del Pinatar. Utilizo código Python para acceder a los datos a través de Motuclient. Aquí puedes cambiar mis credenciales de Copernicus por las tuyas.

Nota

Sí, en un documento Quarto, puedes ejecutar código Python utilizando el encabezado {python} para indicar que es código de Python. Asegúrate de tener el entorno de Python configurado en tu entorno Quarto para poder ejecutar el código. Alternativamente, en RStudio, puedes abrir un nuevo script de Python desde la pestaña “Archivo nuevo” y ejecutarlo en el entorno de Python de RStudio.

Recuerda que debes tener instalado Python y las bibliotecas necesarias en tu entorno local para poder ejecutar el código sin problemas.

import motuclient

!python -m motuclient --motu https://my.cmems-du.eu/motu-web/Motu --service-id MEDSEA_MULTIYEAR_PHY_006_004-TDS --product-id med-cmcc-tem-rean-d --longitude-min -0.64 --longitude-max -0.51 --latitude-min 37.8 --latitude-max 36.9 --date-min "1987-07-01 00:00:00" --date-max "2020-07-01 23:59:59" --depth-min 0 --depth-max 40 --variable thetao --out-dir data --out-name thetao_1987_2020_0_40m_san_pedro_pinatar.nc --user jatalah --pwd Estela2012

4 Data exploration

Luego que me he bajado el fichero .nc en la carpeta data, los puedo leer usando la funcion read_ncdf del paquete stars . Ademas se muestra el resumen del objeto R Stars con 4 dimensiones y 1 atributo. El atributo, en este caso, es la temperatura del agua (thetao) en grados Celsius (°C). El resumen proporciona información estadística sobre el conjunto de datos, como los valores mínimos, primer cuartil, mediana, media, tercer cuartil y máximo.

Además, se proporciona información sobre las dimensiones del conjunto de datos, que incluyen la longitud (lon), latitud (lat), profundidad (depth) y tiempo (time). Cada dimensión tiene un rango específico, y se proporciona información adicional sobre el desplazamiento y el paso en cada dimensión.

sitio1 <-
  read_ncdf(
    'data/thetao_1987_2020_0_40m_san_pedro_pinatar.nc',
    var = 'thetao',
    make_units = F
  )

sitio1
stars object with 4 dimensions and 1 attribute
attribute(s), summary of first 1e+05 cells:
            Min.  1st Qu.   Median    Mean  3rd Qu.     Max.
thetao  14.18145 18.30357 21.29844 20.9303 23.61457 26.16651
dimension(s):
      from    to                  offset  delta  refsys
lon      1     4                      NA     NA  WGS 84
lat      1    23                      NA     NA  WGS 84
depth    1    13                      NA     NA      NA
time     1 12055 1987-07-01 12:00:00 UTC 1 days POSIXct
                        values x/y
lon        [4] -0.625,...,-0.5 [x]
lat   [23] 36.8958,...,37.8125 [y]
depth [13] 1.01824,...,37.8522    
time                      NULL

Los datos se convierten en un tibble para poder manipularlos más fácilmente usando tidyverse

s1 <- 
  sitio1 %>% 
  as_tibble()

s1
# A tibble: 14,417,780 × 5
         lon       lat     depth time                thetao
   <dbl[1d]> <dbl[1d]> <dbl[1d]> <dttm>               <dbl>
 1    -0.625      36.9      1.02 1987-07-01 12:00:00   20.5
 2    -0.583      36.9      1.02 1987-07-01 12:00:00   20.5
 3    -0.542      36.9      1.02 1987-07-01 12:00:00   20.5
 4    -0.5        36.9      1.02 1987-07-01 12:00:00   20.5
 5    -0.625      36.9      1.02 1987-07-01 12:00:00   20.5
 6    -0.583      36.9      1.02 1987-07-01 12:00:00   20.5
 7    -0.542      36.9      1.02 1987-07-01 12:00:00   20.5
 8    -0.5        36.9      1.02 1987-07-01 12:00:00   20.5
 9    -0.625      37.0      1.02 1987-07-01 12:00:00   20.5
10    -0.583      37.0      1.02 1987-07-01 12:00:00   20.5
# ℹ 14,417,770 more rows

Se ve que hay datos temperature (thetao) para cada profundidad y tiempo, con un total de 14,417,780 filas. Para simplificar los datos, calculo la temperatura promedio para cada fecha y profundidad.

s1_mean <- 
  s1 %>% 
  group_by(time, depth) %>%
  summarise(thetao = mean(thetao), .groups = 'drop')
  
s1_mean
# A tibble: 156,715 × 3
   time                    depth thetao
   <dttm>              <dbl[1d]>  <dbl>
 1 1987-07-01 12:00:00      1.02   20.8
 2 1987-07-01 12:00:00      3.17   20.7
 3 1987-07-01 12:00:00      5.46   20.5
 4 1987-07-01 12:00:00      7.92   20.4
 5 1987-07-01 12:00:00     10.5    20.3
 6 1987-07-01 12:00:00     13.3    20.1
 7 1987-07-01 12:00:00     16.3    19.7
 8 1987-07-01 12:00:00     19.4    19.3
 9 1987-07-01 12:00:00     22.7    18.7
10 1987-07-01 12:00:00     26.2    18.1
# ℹ 156,705 more rows

Ahora podemos visualizar la serie de tiempo de temperatura para cada una de las 14 profundidades en el sitio 1. La gráfica muestra la estacionalidad de los datos a través de la serie de tiempo y los cambios de temperatura con la profundidad.

ggplot(s1_mean) + 
  geom_path(aes(time, thetao, color = depth)) +
  scale_color_viridis_c(option = 'B')

Ahora se pueden filtrar los datos de acuerdo con la hipótesis de trabajo. Por ejemplo, buscar la profundidad máxima para cada año en que se registraron temperaturas mayores a 25°C. Para visualizar las tendencias temporales, añadí una curva de tipo loess utilizando geom_smooth. En la gráfica se puede observar que la isoterma de 25°C generalmente alcanza una profundidad máxima de alrededor de 10 metros; sin embargo, en los últimos años esta profundidad rondó los 20 metros.

s1_mean %>% 
  filter(thetao>25) %>% 
  group_by(year(time)) %>% 
  slice_max(depth, n = 1) %>% 
  ggplot(aes(time, depth)) + 
  geom_point() +
  geom_smooth() +
  labs( y = "Profundidad máxima (m) de la isoterma de 25°C", x = "Año")

Otro ejemplo es encontrar el primer día del año para cada año en que se registraron temperaturas mayores a 25°C. La gráfica muestra una tendencia de que las temperaturas mayores a 25°C se alcanzan cada vez más temprano en el año, con una diferencia de alrededor de 20 días más temprano entre el inicio y el final de la serie de tiempo.

s1_mean %>% 
  filter(thetao>25) %>% 
  group_by(year = year(time)) %>% 
  mutate(jday = yday(time)) %>% 
  slice_min(jday, n = 1) %>% 
  ggplot(aes(year, jday )) + 
  geom_point() +
  geom_smooth() +
  labs( y = "Día juliano", x = "Año")

5 Obtención de datos a traves del paquete CopernicusMarine

El paquete CopernisuMarine permite seleccionar y descargar datos Copernicus desde sin necesidad de software externo (por ejmeplo Python).

CRAN version

install.packages("CopernicusMarine")

Development version on github

devtools::install_github('pepijn-devries/CopernicusMarine')

Este código en R utiliza las bibliotecas CopernicusMarine y leaflet para visualizar datos del Servicio Marino de Copernicus en un mapa interactivo utilizando Leaflet. Se crea un mapa interactivo que muestra datos de temperatura del agua del mar en la región especificada utilizando datos del Servicio Marino de Copernicus.

Primero cargo el paquete e introduzco mis credenciales de Copernicus.

library(CopernicusMarine)
options(CopernicusMarine_uid = "jatalah")
options(CopernicusMarine_pwd = "*********")
library(leaflet)

leaflet::leaflet() %>%
  leaflet::setView(lng = -8.8, lat = 42.3, zoom = 4) %>%
  leaflet::addProviderTiles("Esri.WorldImagery") %>%
  addCopernicusWMSTiles(
    product     = "MEDSEA_MULTIYEAR_PHY_006_004",
    layer       = "med-cmcc-tem-rean-d",
    variable    = "thetao"
  )

Se pueden inspeccionar los detalles de la capa de interés en este caso, MEDSEA_MULTIYEAR_PHY_006_004.

copernicus_product_details(
  product  = "MEDSEA_MULTIYEAR_PHY_006_004",
  layer    = "med-cmcc-tem-rean-d",
  variable = "thetao"
)
$id
[1] "med-cmcc-tem-rean-d/thetao"

$subdatasetId
[1] "med-cmcc-tem-rean-d"

$subdatasetTitle
[1] "Potential temperature (3d) - daily mean"

$variableId
[1] "thetao"

$variableTitle
[1] "Sea water potential temperature"

$units
[1] "degrees_C"

$bbox
[1] -6.00000 30.18750 36.29167 45.97917

$tValues
[1] "1987-01-01T12:00:00.000Z/2021-06-30T12:00:00.000Z/P1D"

$zValues
[1] "-1.0182366371154785,-3.1657474040985107,-5.464963436126709,-7.920377254486084,-10.536603927612305,-13.318384170532227,-16.270586013793945,-19.398210525512695,-22.706392288208008,-26.20039939880371,-29.885643005371094,-33.76767349243164,-37.85219192504883,-42.14503860473633,-46.6522102355957,-51.379859924316406,-56.334285736083984,-61.52195739746094,-66.94949340820312,-72.62368774414062,-78.55149841308594,-84.74004364013672,-91.1966323852539,-97.92872619628906,-104.94397735595703,-112.25020599365234,-119.85543060302734,-127.76783752441406,-135.9958038330078,-144.5478973388672,-153.43284606933594,-162.6596221923828,-172.2373504638672,-182.17535400390625,-192.48313903808594,-203.17044067382812,-214.24716186523438,-225.72340393066406,-237.60946655273438,-249.9158477783203,-262.6532287597656,-275.83251953125,-289.46478271484375,-303.5613098144531,-318.133544921875,-333.1931457519531,-348.751953125,-364.82196044921875,-381.4154357910156,-398.5447082519531,-416.2223205566406,-434.4610595703125,-453.2737731933594,-472.6734924316406,-492.6734619140625,-513.2869873046875,-534.527587890625,-556.4088745117188,-578.944580078125,-602.1486206054688,-626.034912109375,-650.6175537109375,-675.9107055664062,-701.9286499023438,-728.6856079101562,-756.196044921875,-784.4743041992188,-813.5348510742188,-843.3921508789062,-874.0606689453125,-905.5548095703125,-937.8890991210938,-971.077880859375,-1005.135498046875,-1040.0762939453125,-1075.914306640625,-1112.6636962890625,-1150.33837890625,-1188.9521484375,-1228.518798828125,-1269.0517578125,-1310.564208984375,-1353.0693359375,-1396.5799560546875,-1441.108642578125,-1486.6678466796875,-1533.2694091796875,-1580.9251708984375,-1629.6466064453125,-1679.44482421875,-1730.330322265625,-1782.3135986328125,-1835.404541015625,-1889.6126708984375,-1944.9471435546875,-2001.4166259765625,-2059.029052734375,-2117.79248046875,-2177.714111328125,-2238.80029296875,-2301.0576171875,-2364.49169921875,-2429.107666015625,-2494.91015625,-2561.903076171875,-2630.08984375,-2699.4736328125,-2770.056640625,-2841.8408203125,-2914.826904296875,-2989.015869140625,-3064.407470703125,-3141.00146484375,-3218.796142578125,-3297.790283203125,-3377.9814453125,-3459.3662109375,-3541.94189453125,-3625.703857421875,-3710.6474609375,-3796.76806640625,-3884.0595703125,-3972.51611328125,-4062.13037109375,-4152.89599609375,-4244.80419921875,-4337.84765625,-4432.017578125,-4527.30419921875,-4623.69873046875,-4721.19140625,-4819.77099609375,-4919.42724609375,-5020.1494140625,-5121.92578125,-5224.74462890625,-5328.59375,-5433.4609375,-5539.33349609375,-5646.19921875,-5754.0439453125"

$zUnits
[1] "m"

$valueMin
[1] -3.15

$valueMax
[1] 36.85

$colormapId
[1] "rainbow"

$logScale
[1] FALSE

$hasVectorStyle
[1] FALSE

$wmsUrl
[1] "https://my.cmems-du.eu/thredds/wms/med-cmcc-tem-rean-d"

$subsetUrl
[1] "https://my.cmems-du.eu/motu-web/Motu"

$subsetAttrs
$subsetAttrs$depth
$subsetAttrs$depth$units
[1] "m"

$subsetAttrs$depth$min
[1] 1.018237

$subsetAttrs$depth$max
[1] 5754.044

$subsetAttrs$depth$step
NULL


$subsetAttrs$latitude
$subsetAttrs$latitude$units
[1] "degrees_north"

$subsetAttrs$latitude$min
[1] 30.1875

$subsetAttrs$latitude$max
[1] 45.97917

$subsetAttrs$latitude$step
[1] 0.04166667


$subsetAttrs$longitude
$subsetAttrs$longitude$units
[1] "degrees_east"

$subsetAttrs$longitude$min
[1] -6

$subsetAttrs$longitude$max
[1] 36.29167

$subsetAttrs$longitude$step
[1] 0.04166667


$subsetAttrs$time
$subsetAttrs$time$units
[1] "milliseconds since 1970-01-01"

$subsetAttrs$time$min
[1] 536500800000

$subsetAttrs$time$max
[1] 1.625054e+12

$subsetAttrs$time$step
NULL



$subsetVariableIds
[1] "thetao"

Luego me bajo los datos de un cuadrado al frente de la costa de Murcia

copernicus_download_motu(
  destination = "data/murcia_test.nc",
  product = "MEDSEA_MULTIYEAR_PHY_006_004",
  layer = "med-cmcc-tem-rean-d",
  variable = "thetao",
  output = "netcdf",
  region = c(-0.64, 37.8, -0.51, 37.9),
  timerange = c("2020-11-01", "2020-11-02"),
  verticalrange = c(0, 50),
  overwrite = TRUE
)

Leo los datos con el paquete stars y los visualizo con ggplot

murcia_test <- stars::read_ncdf('data/murcia_test.nc', var = 'thetao', make_units = F)

ggplot() +
  geom_stars(data = murcia_test, sf = T) +
  theme_void() +
  scale_fill_viridis_c(option = "H", name = "ºC") +
  facet_wrap(~round(depth, 1))

6 Varias areas de interes (AOI)

Preparo un pipeline para bajr varias AOI a laz y visulizarlas

aoi <- 
  tibble(
    xmin = c(-0.64, 21.019739336599876, 33.22144124751825),
    ymin = c(37.8, 38.33459744736376, 34.65349160442335),
    xmax = c(-0.51, 21.107624652007036, 33.27774617915888),
    ymax = c(37.9, 38.41161035261746, 34.7056766630171),
    dest = c("murcia", "greece", "creete")
  )

Creo una función para bajar los datos en diferentes AOI donde se desarrolla acuicultura de peces. Las áreas se pueden seleccionar desde un mapa o Ocean Data View para poder ver donde hay datos disponibles.

motu_dowload_func <- function(xmin, ymin, xmax, ymax, dest) {
    copernicus_download_motu(
      destination = paste0("data/", dest, ".nc"),
      product = "MEDSEA_MULTIYEAR_PHY_006_004",
      layer = "med-cmcc-tem-rean-d",
      variable = "thetao",
      output = "netcdf",
      region = c(xmin, ymin, xmax, ymax),
      timerange = c("2020-11-01", "2020-11-02"),
      verticalrange = c(0, 50),
      overwrite = TRUE
    )
  }

Aplico la funcion a cada fila del tibble que contiene las AOI.

pmap(aoi, motu_dowload_func)
[[1]]
[1] TRUE

[[2]]
[1] TRUE

[[3]]
[1] TRUE

Una vez descargados los datos, puedo leerlos con una función creada para ello.

nc_read_func <- function(path){stars::read_ncdf(path, var = 'thetao', make_units = F)}

murcia <- nc_read_func('data/murcia.nc')
greece <-nc_read_func('data/greece.nc')
creete <- nc_read_func('data/creete.nc')

Finalmente los puedo visualizar los datos para cada region.

nc_plot_func <- function(data) {
  ggplot() +
    geom_stars(data = data, sf = T) +
    theme_void() +
    scale_fill_viridis_c(option = "H", name = "ºC") +
    facet_wrap( ~ round(depth, 1))
}

nc_plot_func(murcia)

nc_plot_func(greece)

nc_plot_func(creete)

7 Tarea

Replicar estos ejemplos para otras áreas del Mediterráneo donde se practica la acuicultura de peces.

Tip

Eligir otras areas utilzando el mapa en este articulo:

Trujillo, P., Piroddi, C., & Jacquet, J. (2012). Fish farms at sea: the ground truth from Google Earth. PLoS One, 7(2), e30546. https://doi.org/10.1371/journal.pone.0030546