Introduction to {resourcecode}

The goal of {resourcecode} is to provide an easy access to the ResourceCODE hindcast database of sea-states. More information on the database can be found here. We will explain in this vignette the way to retrieve data from this database and some the functionalities offered by this package.

This package depends on data in a data package {resourcecodedata} that is available through a drat repository on GitHub. To use the {resourcecode} package, you will need to install {resourcecodedata} on your computer, using the following install.packages function (and later update it using the update.packages function):

install.packages("resourcecodedata",
  repos = "https://resourcecode-project.github.io/drat/",
  type = "source"
)

Once the packages are installed, we can proceed as usual:

library(resourcecodedata)
library(resourcecode)
library(ggplot2)

Database configuration: nodes, bathymetry…

There are two grids in the ResourceCODE hindcast database: the full, dense, FIELD grid, which is the grid where the numerical WWIII model is run and where sea-state parameters and only 1D spectra are available. The second grid, SPEC, is coarser, but contains the full 2D spectral data at each node.

The variable rscd_field contains the coordinates of the FIELD grid, along with the depth of these points and d50, the bottom sediment type.

str(rscd_field)
#> 'data.frame':    328030 obs. of  5 variables:
#>  $ node     : int  1 2 3 4 5 6 7 8 9 10 ...
#>  $ longitude: num  6.92 6.91 6.92 6.93 6.95 ...
#>  $ latitude : num  62.9 62.9 62.9 62.9 62.9 ...
#>  $ depth    : num  3 3 3 3 3 3 3 3 3 3 ...
#>  $ d50      : num  0 1.32 1.32 1.32 1.32 ...
#>  - attr(*, ".internal.selfref")=<externalptr>
head(rscd_field)
#>   node longitude latitude depth   d50
#> 1    1  6.922498 62.91467     3 0.000
#> 2    2  6.910178 62.91301     3 1.322
#> 3    3  6.920552 62.90438     3 1.322
#> 4    4  6.933673 62.89938     3 1.322
#> 5    5  6.945979 62.90012     3 1.322
#> 6    6  6.958059 62.89899     3 1.322

At each location, a large number of sea-state parameters are available, described in rscd_variables:

str(rscd_variables)
#> 'data.frame':    88 obs. of  3 variables:
#>  $ name    : chr  "longitude" "latitude" "tri" "MAPSTA" ...
#>  $ longname: chr  "longitude" "latitude" "tri" "status map" ...
#>  $ unit    : chr  "degree_east" "degree_north" "" "1" ...
head(rscd_variables)
#>        name         longname         unit
#> 1 longitude        longitude  degree_east
#> 2  latitude         latitude degree_north
#> 3       tri              tri             
#> 4    MAPSTA       status map            1
#> 5       dpt            depth            m
#> 6      ucur eastward current        m s-1

The available location can be plotted on a map using the following code. It can be noticed here that the variables rscd_coastline and rscd_islands contain the path of the coast line and the islands, respectively.

lim_lon <- c(-5.25, -4.25)
lim_lat <- c(47.75, 48.75)
field_bzh <- ggplot(rscd_field, aes(x = longitude, y = latitude)) +
  geom_point(size = .1, col = "lightblue") +
  geom_path(data = rscd_coastline, linewidth = .2) +
  geom_path(data = rscd_islands, aes(group = .data$ID), linewidth = .2) +
  coord_sf(xlim = lim_lon, ylim = lim_lat, expand = FALSE, crs = sf::st_crs(4326)) +
  theme_void()
field_bzh

The spectral grid, or SPEC grid is coarser, with more than 24,000 nodes where the full 2D spectra are available. The grid can be added to the previous plot to see the differences in the spatial coverage.

str(rscd_spectral)
#> 'data.frame':    24276 obs. of  5 variables:
#>  $ longitude: num  1.5 -5.5 -6 -6.5 -7 4.5 4 3.5 3 2.5 ...
#>  $ latitude : num  52 52 52 52 52 52.5 52.5 52.5 52.5 52.5 ...
#>  $ name     : chr  "E001500N52000" "W005500N52000" "W006000N52000" "W006500N52000" ...
#>  $ depth    : num  15.4 89.6 95.8 52.7 54.7 ...
#>  $ d50      : num  0.737 -0.8265 -1.585 0.737 0.0346 ...
#>  - attr(*, ".internal.selfref")=<externalptr>
head(rscd_spectral)
#>   longitude latitude          name    depth         d50
#> 1       1.5     52.0 E001500N52000 15.41327  0.73700004
#> 2      -5.5     52.0 W005500N52000 89.64895 -0.82650555
#> 3      -6.0     52.0 W006000N52000 95.83614 -1.58500008
#> 4      -6.5     52.0 W006500N52000 52.67409  0.73700004
#> 5      -7.0     52.0 W007000N52000 54.67464  0.03464983
#> 6       4.5     52.5 E004500N52500 14.18527  0.73700004
field_bzh + geom_point(data = rscd_spectral, col = "orange", size = .1)

When one is interested in some location, we also provide helpers to find the closest point in each of the grid, via the closest_point_FIELD() and closest_point_SPEC() functions.

Getting data

We have two different sources of data, the FIELD grid, and the SPEC grid. We will cover each use case here.

For example, if one is interested in the time series of sea-state parameters at some location, the following code can be adapted:

point_of_interest <- c(longitude = -4.6861533, latitude = 48.3026514)
node <- closest_point_field(point_of_interest)
node
#> $point
#> [1] 134940
#> 
#> $distance
#> [1] 288.974
ts <- get_parameters(node = node$point, parameters = c("hs", "tp", "dp", "cge"))
ggplot(tidyr::pivot_longer(ts, -1), aes(x = time, y = value, col = name)) +
  geom_line() +
  coord_cartesian(expand = FALSE) +
  facet_wrap(~name, ncol = 2, scales = "free_y") +
  scale_x_datetime(name = NULL, date_breaks = "month") +
  scale_y_continuous(name = NULL) +
  theme_minimal() +
  theme(
    legend.position = "none",
    axis.text.x = element_text(angle = 60, hjust = 1)
  )

node_spectral_grid <- closest_point_spec(point_of_interest)

Then one can extract the 1D spectrum time-series and see the contents of these downloaded data.

spec_1d <- get_1d_spectrum(node_spectral_grid$point, start = "1994-01-01", end = "1994-02-28")
str(spec_1d)
#> List of 12
#>  $ longitude : num -4.68
#>  $ latitude  : num 48.3
#>  $ frequency1: num [1:36(1d)] 0.0339 0.0356 0.0392 0.0431 0.0474 ...
#>  $ frequency2: num [1:36(1d)] 0.0356 0.0392 0.0431 0.0474 0.0521 ...
#>  $ ef        : num [1:36, 1:1416] 0.000526 0.0028 0.009767 0.145907 1.126401 ...
#>   ..- attr(*, "dimnames")=List of 2
#>   .. ..$ : NULL
#>   .. ..$ : NULL
#>  $ th1m      : num [1:36, 1:1416] 202 203 205 226 230 ...
#>   ..- attr(*, "dimnames")=List of 2
#>   .. ..$ : NULL
#>   .. ..$ : NULL
#>  $ th2m      : num [1:36, 1:1416] 202 203 205 180 180 ...
#>   ..- attr(*, "dimnames")=List of 2
#>   .. ..$ : NULL
#>   .. ..$ : NULL
#>  $ sth1m     : num [1:36, 1:1416] 13.5 13.9 14.5 16.6 14.7 ...
#>   ..- attr(*, "dimnames")=List of 2
#>   .. ..$ : NULL
#>   .. ..$ : NULL
#>  $ sth2m     : num [1:36, 1:1416] 26.1 26.8 28 32.2 28.7 ...
#>   ..- attr(*, "dimnames")=List of 2
#>   .. ..$ : NULL
#>   .. ..$ : NULL
#>  $ freq      : num [1:36(1d)] 0.0339 0.0373 0.041 0.0451 0.0496 ...
#>  $ forcings  : tibble [1,416 × 14] (S3: tbl_df/tbl/data.frame)
#>   ..$ time  : POSIXct[1:1416], format: "1994-01-01 00:00:00" "1994-01-01 01:00:00" ...
#>   ..$ dpt   : num [1:1416] 27 27.5 29 30.5 31.5 32.5 32.5 32 31 29.5 ...
#>   ..$ wnd   : num [1:1416] 12.2 11.6 10.8 10.9 10.8 ...
#>   ..$ wnddir: num [1:1416] 301 299 298 295 294 ...
#>   ..$ cur   : num [1:1416] 0.1 0.2 0.3 0.4 0.3 ...
#>   ..$ curdir: num [1:1416] 325 27.5 41.2 48.6 55.9 ...
#>   ..$ hs    : num [1:1416] 3.02 3.11 3.23 3.35 3.43 ...
#>   ..$ fp    : num [1:1416] 0.0672 0.067 0.0671 0.0672 0.0727 ...
#>   ..$ f02   : num [1:1416] 0.132 0.129 0.124 0.121 0.12 ...
#>   ..$ f0m1  : num [1:1416] 0.0905 0.0882 0.0864 0.0858 0.0862 ...
#>   ..$ th1p  : num [1:1416] 234 234 235 237 238 ...
#>   ..$ sth1p : num [1:1416] 13.4 13.3 13.1 12.9 12.7 ...
#>   ..$ dir   : num [1:1416] 239 239 239 240 241 ...
#>   ..$ spr   : num [1:1416] 19.6 19.5 18.8 18 17.4 ...
#>  $ station   : chr "W004679N48304"

The same applies to 2D spectra.

spec_2d <- get_2d_spectrum(node_spectral_grid$point, start = "1994-01-01", end = "1994-02-28")
str(spec_2d)
#> List of 9
#>  $ longitude : num -4.68
#>  $ latitude  : num 48.3
#>  $ frequency1: num [1:36(1d)] 0.0339 0.0356 0.0392 0.0431 0.0474 ...
#>  $ frequency2: num [1:36(1d)] 0.0356 0.0392 0.0431 0.0474 0.0521 ...
#>  $ efth      : num [1:36, 1:36, 1:1416] 0.00024 0.000629 0.000872 0.000705 0.000423 ...
#>   ..- attr(*, "dimnames")=List of 3
#>   .. ..$ : NULL
#>   .. ..$ : NULL
#>   .. ..$ : NULL
#>  $ freq      : num [1:36(1d)] 0.0339 0.0373 0.041 0.0451 0.0496 ...
#>  $ dir       : num [1:36(1d)] 0 10 20 30 40 50 60 70 80 90 ...
#>  $ forcings  : tibble [1,416 × 6] (S3: tbl_df/tbl/data.frame)
#>   ..$ time  : POSIXct[1:1416], format: "1994-01-01 00:00:00" "1994-01-01 01:00:00" ...
#>   ..$ dpt   : num [1:1416] 27 27.5 29 30.5 31.5 32.5 32.5 32 31 29.5 ...
#>   ..$ wnd   : num [1:1416] 12.2 11.6 10.8 10.9 10.8 ...
#>   ..$ wnddir: num [1:1416] 301 299 298 295 294 ...
#>   ..$ cur   : num [1:1416] 0.1 0.2 0.3 0.4 0.3 ...
#>   ..$ curdir: num [1:1416] 145 208 221 229 236 ...
#>  $ station   : chr "W004679N48304"

In addition, it is possible to plot the wave elevation directional spectra for any given time, which can be specified by the time index or directly the date:

plot_2d_specta(spec_2d, "1994-01-15 18:00")