Usando Python desde R en el contexto de Copernicus

Interfaz de R a Python usando Reticulate

El paquete reticulate proporciona una interfaz R para módulos, clases y funciones de Python.

https://rstudio.github.io/reticulate/

El primer paso es instalar Python a traves de por ejemplo Anaconda. https://www.anaconda.com/products/distribution

Luego instalar el paquete reticulate desde CRAN usando install.packages('reticulate')

Version de Python

De manera predeterminada, reticulate usa la versión de Python que se encuentra en el PATH Sys.which("python")

La función use_python() te permite especificar una versión alternativa, por ejemplo:

library(reticulate)

use_python("/usr/local/bin/python")

quit
library(reticulate)

Funciones básicas usando Python chunks

a = "Hola " + "Mundo"
print(a)

Hola Mundo

Ejemplo de una gráfica Python Matplotlib

import matplotlib.pyplot as plt
import numpy as np

# Data for plotting
t = np.arange(0.0, 2.0, 0.01)
s = 1 + np.sin(2 * np.pi * t)

fig, ax = plt.subplots()
ax.plot(t, s)

[<matplotlib.lines.Line2D object at 0x000001F65ADF9CA0>]

ax.set(xlabel='time (s)', ylabel='voltage (mV)')

[Text(0.5, 0, 'time (s)'), Text(0, 0.5, 'voltage (mV)')]

ax.grid()

plt.show()

Acceder a Copernicus Copernicus Marine Service Catalogue

Si no tiene motuclient en Python, deberá instalarlo usando:

python -m pip install motuclient==1.8.4 --no-cache-dir

Deberás reemplazar tus credenciales de inicio de sesión de Copernicus en el código de abojo y correlo en un chuck de Python.

También puedes cambiar el nombre del archivo y el directorio de salida según sus necesidades.

import motuclient

!python -m motuclient --motu https://nrt.cmems-du.eu/motu-web/Motu --service-id MEDSEA_ANALYSISFORECAST_PHY_006_013-TDS --product-id med-cmcc-tem-an-fc-h --longitude-min -17.2917 --longitude-max 36.2917 --latitude-min 30.1875 --latitude-max 45.9792 --date-min "2022-06-18 23:30:00" --date-max "2022-06-18 23:30:00" --depth-min 1.0182 --depth-max 1.0183 --variable thetao --out-dir nc_files --out-name temp_18_june.nc --user jatalah --pwd ***********

Convertir fichero NetCDF a tibble

Primero se instala el paquete tidync usando:

install.packages("tidync")

Ver tutorial en este blog https://ropensci.org/blog/2019/11/05/tidync/

quit
library(tidync)

sst_june_22 <- tidync('nc_files/temp_18_june.nc')

sst_june_22

Data Source (1): temp_18_june.nc ...

Grids (5) <dimension family> : <associated variables> 

[1]   D3,D2,D1,D0 : thetao    **ACTIVE GRID** ( 489060  values per variable)
[2]   D0          : time
[3]   D1          : depth
[4]   D2          : lat
[5]   D3          : lon

Dimensions 4 (all active): 
  
  dim   name  length         min         max start count        dmin        dmax unlim coord_dim 
  <chr> <chr>  <dbl>       <dbl>       <dbl> <int> <int>       <dbl>       <dbl> <lgl> <lgl>     
1 D0    time       1 64409730    64409730        1     1 64409730    64409730    FALSE TRUE      
2 D1    depth      1        1.02        1.02     1     1        1.02        1.02 FALSE TRUE      
3 D2    lat      380       30.2        46.0      1   380       30.2        46.0  FALSE TRUE      
4 D3    lon     1287      -17.3        36.3      1  1287      -17.3        36.3  FALSE TRUE      

Explorar metadata

ncmeta::nc_grids('nc_files/temp_18_june.nc') %>% tidyr::unnest(cols = c(variables))

ncmeta::nc_vars('nc_files/temp_18_june.nc') 

Obten la información de las unidades de tiempo y luego utilízala para convertir los valores sin procesar en una fecha y hora

ncmeta::nc_atts('nc_files/temp_18_june.nc', "time")  %>%
  dplyr::select(value) %>% 
  flatten_df() %>% 
  pull(units)

[1] "minutes since 1900-01-01 00:00:00"

temp_june_22_t <- 
  sst_june_22 %>% 
  hyper_tibble() %>% 
  mutate(time = lubridate::ymd_hms("1900-01-01 00:00:00", tz = "UTC") + time * 60)

head(temp_june_22_t)

Grafica usando ggplot

# get countries data for the map
library(rnaturalearth)
library(rnaturalearthdata)
countries <- ne_countries(scale = "medium", returnclass = "sf")

# plot 
ggplot() +
  # add raster layer
  geom_tile(aes(x = lon, y = lat, fill = thetao), data = temp_june_22_t) +
  # define color palette of raster layer
  scale_fill_viridis_c(option = 'A', name = "°C") +
  # add countries layers
  geom_sf(
    fill = grey(0.9),
    color = grey(0.6),
    lwd = 0.2,
    data = countries
  ) +
  # define spatial extent
  coord_sf(
    xlim = range(temp_june_22_t$lon),
    ylim = range(temp_june_22_t$lat),
    expand = F,
    ndiscr = 500
  ) +
  labs(subtitle = "SST -  Sábado 18 Junio 2022",
       y = "Latitude", x = 'Longitude') +
  theme_minimal() 

Gráfica Python

import matplotlib.pyplot as plt
import xarray as xr

C:\Users\javiera\MINICO~1\lib\site-packages\xarray\backends\cfgrib_.py:27: UserWarning: Failed to load cfgrib - most likely there is a problem accessing the ecCodes library. Try `import cfgrib` to get the full error message
  warnings.warn(

import warnings
warnings.filterwarnings('ignore')


# satellite observation dataset
path_sat = 'nc_files/temp_18_june.nc'

# Open the files and store them in a Python variable
sat = xr.open_dataset(path_sat)
print(sat)

<xarray.Dataset>
Dimensions:  (depth: 1, time: 1, lat: 380, lon: 1287)
Coordinates:
  * depth    (depth) float32 1.018
  * lon      (lon) float32 -17.29 -17.25 -17.21 -17.17 ... 36.21 36.25 36.29
  * time     (time) datetime64[ns] 2022-06-18T23:30:00
  * lat      (lat) float32 30.19 30.23 30.27 30.31 ... 45.85 45.9 45.94 45.98
Data variables:
    thetao   (time, depth, lat, lon) float32 ...
Attributes:
    title:                           Potential Temperature (3D) - Hourly Mean
    FROM_ORIGINAL_FILE__field_type:  hourly_mean_centered_at_time_field
    source:                          MFS EAS6
    institution:                     Centro Euro-Mediterraneo sui Cambiamenti...
    contact:                         servicedesk.cmems@mercator-ocean.eu
    references:                      Clementi, E., Aydogdu, A., Goglio, A. C....
    comment:                         Please check in CMEMS catalogue the INFO...
    Conventions:                     CF-1.0
    bulletin_date:                   20220615
    bulletin_type:                   forecast
    _CoordSysBuilder:                ucar.nc2.dataset.conv.CF1Convention
    history:                         Data extracted from dataset http://local...

sat.thetao.sel(time='2022-06-18').plot()

<matplotlib.collections.QuadMesh object at 0x000001BB951A30D0>

plt.show()

Plot vertical profiles

Ahora usamos los datos del modelo para trazar perfiles verticales de clorofila.

Primero leemos los datos y comprobamos las unidades de tiempo para corregirlos.

path <- 'nc_files/mod_bio.nc'
mod_bio_dat <- tidync(path)

ncmeta::nc_grids(path) %>% tidyr::unnest(cols = c(variables))

ncmeta::nc_vars(path)

ncmeta::nc_atts(path, "time")  %>%
  dplyr::select(value) %>% 
  flatten_df() %>% 
  pull(units)

[1] "seconds since 1970-01-01 00:00:00"

Luego definimos un localidad y dos fechas para la gráfica usando la función hyper_filter del paquete tidync.

mod_bio_dat_t <- 
mod_bio_dat %>% 
  hyper_filter(longitude  = between(longitude, 5, 5), 
               latitude = between(latitude, 42.02, 42.03)) %>% 
  hyper_tibble() %>% 
  mutate(time = lubridate::ymd_hms("1970-01-01 00:00:00", tz = "UTC") + time,
         month = lubridate::month(time))

We can now plot the vertical profile using ggplot

ggplot(mod_bio_dat_t, aes(depth, chl, color = factor(month))) +
  geom_path () +
  coord_flip() +
  scale_x_reverse() +
  labs(x= 'Depth (m)', y = Chl-italic(a)~(mg/m^3), parse = T) +
  theme_minimal() +
  scale_color_discrete(name = "month")

Trazar cortes verticales

Este script se utiliza para trazar la temperatura en un corte vertical a lo largo de una línea longitudinal.

Primero definimos las coordenadas del corte vertical usando la función hyper_filter.

Luego lo convertimos en un tibble usando hyper_tibble y convertimos el formato de hora.

## plot vertical slice
mod_bio_dat_slice <-
  mod_bio_dat %>%
  hyper_filter(
    longitude  = between(longitude, 13.2, 16.5),
    latitude = between(latitude, 43.68, 43.7)
  ) %>%
  hyper_tibble() %>%
  mutate(
    time = lubridate::ymd_hms("1970-01-01 00:00:00", tz = "UTC") + time,
    month = lubridate::month(time)
  )

Ahora usamos geom_tile del paquete ggplot para visualizar el perfil vertical de temperatura a lo largo del corte longitudinal.

ggplot(mod_bio_dat_slice,aes(depth, longitude, fill  = chl ))+
  geom_tile(width = 8) +
  coord_flip() +
  scale_x_reverse() +
  labs(x = 'Depth (m)',
       y = "Longitude",
       parse = T) +
  theme_minimal() +
  scale_fill_viridis_c(name = Chl - italic(a) ~ (mg / m ^ 3))