Código
library(tidyverse, quietly = T)
library(tidync)
library(knitr)Copernicus Marine Service Mediterranean Sea training
A pipeline to extract and process Copernicus NetCDF data in a tidy format.
The script is based on the library tydinc. For details of the and tutorial see this blog.
The examples are adapted from Copernicus Marine Training Workshop for the Mediterranean sea - Edition 2022.
First download the SST data (or whatever data you are interested in) from the Copernicus Marine Service Catalogue using Python and the code below. You will need to replace your Copernicus login credentials in the code. You can also change the file name and output directory to your needs.
Python can be accessed through Rstudio using the reticulate library.
If you don’t have motuclient in Python, you will need to install it using:
!python -m pip install motuclient==1.8.4 --no-cache-dir
Código
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 *******1 Import data with tidync and check variables
Código
sst_may_22 <- tidync('nc_files/temp_18_june.nc')
ncmeta::nc_grids('nc_files/temp_18_june.nc') %>%
tidyr::unnest(cols = c(variables)) %>%
kable()| grid | ndims | variable | nvars |
|---|---|---|---|
| D3,D2,D1,D0 | 4 | thetao | 1 |
| D0 | 1 | time | 1 |
| D1 | 1 | depth | 1 |
| D2 | 1 | lat | 1 |
| D3 | 1 | lon | 1 |
Código
ncmeta::nc_vars('nc_files/temp_18_june.nc') %>%
kable()| id | name | type | ndims | natts |
|---|---|---|---|---|
| 0 | depth | NC_FLOAT | 1 | 10 |
| 1 | thetao | NC_FLOAT | 4 | 7 |
| 2 | lon | NC_FLOAT | 1 | 8 |
| 3 | time | NC_DOUBLE | 1 | 9 |
| 4 | lat | NC_FLOAT | 1 | 8 |
2 Obtain the time units information and then use it to convert the raw value into a date-time
Código
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"Código
sst_may_22_t <-
sst_may_22 %>%
hyper_tibble() %>%
mutate(time = lubridate::ymd_hms("1900-01-01 00:00:00", tz = "UTC") + time * 60)3 Plot the data
Código
# 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 = sst_may_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(sst_may_22_t$lon),
ylim = range(sst_may_22_t$lat),
expand = F,
ndiscr = 500
) +
labs(title = "Sea Surface Temperature (SST) - Sunday 22 May 2022",
y = "Latitude") +
theme_minimal() 
4 Vertical profiles
We now use the model data to plot some chlorophyll vertical profiles, to visualize how the chlorophyll concentration changes with depth.
First we read the data and check the time units to correct them.
Código
path <- 'nc_files/mod_bio.nc'
mod_bio_dat <- tidync(path)
ncmeta::nc_grids(path) %>%
tidyr::unnest(cols = c(variables)) %>%
kable()| grid | ndims | variable | nvars |
|---|---|---|---|
| D3,D2,D1,D0 | 4 | chl | 1 |
| D0 | 1 | time | 1 |
| D1 | 1 | depth | 1 |
| D2 | 1 | latitude | 1 |
| D3 | 1 | longitude | 1 |
Código
ncmeta::nc_vars(path) %>%
kable()| id | name | type | ndims | natts |
|---|---|---|---|---|
| 0 | depth | NC_FLOAT | 1 | 9 |
| 1 | chl | NC_FLOAT | 4 | 7 |
| 2 | latitude | NC_FLOAT | 1 | 7 |
| 3 | time | NC_DOUBLE | 1 | 9 |
| 4 | longitude | NC_FLOAT | 1 | 7 |
Código
ncmeta::nc_atts(path, "time") %>%
dplyr::select(value) %>%
flatten_df() %>%
pull(units)[1] "seconds since 1970-01-01 00:00:00"Then we define a location and two dates to plot using the hyper_filter function of the tidync library.
Código
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
Código
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")
5 Vertical slices
This scrip is used to plot temperature on a vertical slice along a line.
First we define the coordinates of the vertical slice using the tidync function hyper_filter.
Then we convert it into a tibble suing hyper_tibble and convert the time format.
Código
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)
)Now we geom_tile from the ggplot library to visualize the temperature vertical profile along the vertical slice.
Código
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))