NZ Copernicus

Author

Javier Atalah

Published

March 28, 2023

Code

library(tidyverse)
library(sf)
library(raster)
library(mapview)
library(RStoolbox)
library(reticulate)
library(tidync)
library(rnaturalearth)
library(s2)
conflicted::conflicts_prefer(dplyr::filter)
sf_use_s2(FALSE)

1 NZ Coastal area 50 km buffer

Code

nz_high <- 
  ne_countries(country = 'new zealand', scale = "large", returnclass = "sf") %>% 
  st_crop(
    xmin = 165,
    xmax = 180,
    ymin = -48,
    ymax = -32
  )

costal_area <- 
  s2_buffer_cells(
    s2_geog_from_wkb(st_as_binary(st_geometry(nz_high)), check = FALSE),
    distance = 5e4,
    max_cells = 1e4
  ) %>%
  st_as_sf() %>% 
  st_difference(nz_high)

bbox <- st_bbox(extent(nz_high))

ggplot() +
  geom_sf(data = costal_area, fill = "lightblue", col = "gray70") +
  geom_sf() +
  coord_sf() +
  theme_minimal()

2 Bathymetry

Download bathymetry data from GEBCO

Code

bathy <-
  raster('data/gebco_2022_n-31.2012_s-50.0098_w161.8945_e181.4062.tif') %>% 
  crop(extent(costal_area)) %>% 
  mask(costal_area) %>%
  fortify(maxpixels = 9e9) %>%
  drop_na() %>%
  as_tibble() %>%
  rename(depth = gebco_2022_n.31.2012_s.50.0098_w161.8945_e181.4062) %>% 
  filter(depth<0) %>% 
  mutate(depth = abs(depth))

Then we can visualise it with ggplot.

Code

ggplot() +
  geom_raster(data = bathy, aes(x, y, fill = depth)) +
  scale_fill_viridis_b(option = 'D', name = "Depth", trans = 'log10') +
  geom_contour(data = bathy, aes(x, y, z = depth), color = 'gray50') +
  coord_sf(expand = FALSE) +
  labs(x = NULL, y = NULL) +
  theme_minimal()

3 SST Copernicus data

First, I download SST data using motuclient for Python.

Code

import motuclient

!python -m motuclient --motu https://my.cmems-du.eu/motu-web/Motu --service-id SST_GLO_SST_L4_REP_OBSERVATIONS_010_011-TDS --product-id METOFFICE-GLO-SST-L4-REP-OBS-SST --longitude-min 166 --longitude-max 180 --latitude-min -48 --latitude-max -33 --date-min "2021-12-31 12:00:00" --date-max "2021-12-31 12:00:00" --variable analysed_sst --out-dir data --out-name  sst_dec_21.nc --user jatalah --pwd **********

The nc file is first converted into a tibble, then into an sf object, so that it can be clipped based on the polygon of the coast of NZ.

Code

# another option
# sst_dat <-  raster("data/sst_dec_21.nc",
#             varname = "analysed_sst") %>% 
#   raster::as.data.frame(xy = TRUE) %>% 
#   rename(sst = analysed.sea.surface.temperature) %>% 
#   st_as_sf(coords = c("x", "y"),
#            crs = 4326,
#            agr = "constant")


sst_dat <- tidync('data/sst_dec_21.nc') 

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

# A tibble: 4 × 4
  grid     ndims variable     nvars
  <chr>    <int> <chr>        <int>
1 D2,D1,D0     3 analysed_sst     1
2 D0           1 time             1
3 D1           1 lat              1
4 D2           1 lon              1

Code

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

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

Code

sst_dat_t <- 
sst_dat %>% 
  hyper_tibble() %>% 
  mutate(time = lubridate::ymd_hms("1981-01-01 00:00:00", tz = "UTC") + time,
         sst = analysed_sst -273,15) 


sst_dat_sf <- 
  sst_dat_t %>%
  st_as_sf(coords = c("lon", "lat"),
           crs = 4326,
           agr = "constant") %>% 
  st_crop(extent(costal_area)) %>%
  st_intersection(costal_area)

ggplot() +
  geom_sf(data = sst_dat_sf, aes(color = sst )) +
  scale_color_viridis_c(name = "°C", option = "H") +
  geom_sf(
    fill = grey(0.9),
    color = grey(0.6),
    lwd = 0.2,
    data = nz_high
  ) +
  coord_sf() +
  labs(y = "Latitude", x = 'Longitude') +
  theme_minimal()

4 Long-term trend in SST

Code

!python -m motuclient --motu https://my.cmems-du.eu/motu-web/Motu --service-id SST_GLO_SST_L4_REP_OBSERVATIONS_010_011-TDS --product-id METOFFICE-GLO-SST-L4-REP-OBS-SST --longitude-min 171 --longitude-max 175 --latitude-min -42 --latitude-max -40 --date-min "2000-01-01 00:00:00" --date-max "2022-01-01 00:00:00" --variable analysed_sst --out-dir nc_files --out-name sst_all_marl.nc --user jatalah --pwd Estela2012

Code

marl_bbox <- 
st_bbox(c(xmin = 172.7,
          xmax = 174.5,
          ymin = -41.3,
          ymax = -40.4))

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

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

Code

sst_raw <- 
  tidync('nc_files/sst_all_marl.nc') %>% 
  hyper_tibble() %>% 
  mutate(time = lubridate::ymd_hms("1981-01-01 00:00:00", tz = "UTC") + time,
         year = year(time),
         month = month(time),
         sst = analysed_sst -273,15) 
  
sum_sst <-
  sst_raw %>%
  group_by(lon, lat, year, month) %>%
  summarise(sst = mean(sst)) # calcular medias mensuales

For the analysis of SST trends, temperatures from the month of February were used as an indicator of the warmest month.

The figure shows Sen’s slopes (increase in degrees per year per pixel) as an indicator of the trend of SST changes over time in the study area.

Code

sum_sst %>% 
  dplyr::filter(month ==1) %>% 
  group_by(lat,lon) %>% 
  nest() %>% 
  mutate(sen = map(data, ~er.helpers::sen_slope(.x$sst))) %>%
  dplyr::select(sen) %>%
  unnest(cols = c('sen')) %>%
  ungroup() %>% 
  st_as_sf(coords = c("lon", "lat"), crs = 4326, agr = "constant") %>%
  ggplot() +
  geom_sf(aes(color = slope ), size = 2.7) +
  # scale_color_viridis_c(option = 'A', name = "Trend") +
  scale_color_viridis_c(name = "°C", option = "H") +
  geom_sf(
    fill = grey(0.9),
    color = grey(0.6),
    lwd = 0.2,
    data = nz_high
  ) +
  coord_sf(xlim = c(170, 175),
           ylim = c(-42, -40),
           expand = TRUE) +
  labs(y = "Latitude", x = 'Longitude') +
  theme_minimal()