Read a .nc file and get variables at certain locations

Question

I have two data files. One of them contains locations with certain spots.

data_locs <- data.frame(
  Year = c(2000, 2001, 1993, 1993, 1994, 1995),
  Month = c(1, 1, 2, 9, 9, 9),
  Day = c(24, 22, 4, 9, 21, 27),
  Long = c(11.3667, 11.6167, -0.2833, -8.9848, -8.9833, -9.3048),
  Lat = c(57.6833, 57.5500, 60.0000, 42.2412, 42.2428, 41.9963),
  Depth = c(68, 29, 141, 99, 99, 170)
)

And on the other side I have a .nc file with data from Copernicus with some variables: temperature, salinity, ocean mix layer thickness. I downloaded an example from Copernicus with reduced spatial and temporal extents that can be downloaded from my github.

What I want to do is extract the values of the raster matching the date (year/month/day), location and depth the get the projected values of the variables. However, all examples I can find online are using ncvar_get to get grided values for latitude, longitude, etc. and then looking for closes locations instead of using the capabilities of the rasters.

It may be a problem of concept but isn't it better to use the projected? How do I extract the values of the locations in data_locs?

Answer 1

I would use {terra} package for it. A kind of workflow might be:

I have changed your data.frame a bit to have a date in the same period as the layers in provided .nc file, the same applies to coordinates (they are within extend of .nc. I have added date variable in the same format as is stored in .nc file (read in by terra:rast() function).

data_locs <- data.frame(
  Year = c("2018", "2000", "2001"),
  Month = c("08", "01", "01"),
  Day = c("12", "24", "22"),
  Long = c(-6, -5, -7),
  Lat = c(45, 46, 47),
  Depth = c(68, 29, 141)
) |>
  dplyr::mutate(date = as.POSIXct(paste0(Year, "-", Month, "-", Day), tz="UTC"))

Let's read the .nc file:

r <- terra::rast("~/Downloads/cmems_mod_glo_phy_my_0.083deg_P1D-m_multi-vars_8.00W-5.00W_45.00N-47.00N_0.49-5727.92m_2018-08-12-2018-08-17.nc")
r
#> class       : SpatRaster 
#> dimensions  : 25, 37, 618  (nrow, ncol, nlyr)
#> resolution  : 0.08333333, 0.08333333  (x, y)
#> extent      : -8.041667, -4.958333, 44.95833, 47.04167  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84) 
#> sources     : cmems_mod_glo_phy_my_0.083deg_P1D-m_multi-vars_8.00W-5.00W_45.00N-47.00N_0.49-5727.92m_2018-08-12-2018-08-17.nc:bottomT  (6 layers) 
#>               cmems_mod_glo_phy_my_0.083deg_P1D-m_multi-vars_8.00W-5.00W_45.00N-47.00N_0.49-5727.92m_2018-08-12-2018-08-17.nc:mlotst  (6 layers) 
#>               cmems_mod_glo_phy_my_0.083deg_P1D-m_multi-vars_8.00W-5.00W_45.00N-47.00N_0.49-5727.92m_2018-08-12-2018-08-17.nc:thetao  (300 layers) 
#>               ... and 2 more sources
#> varnames    : bottomT (Sea floor potential temperature) 
#>               mlotst (Density ocean mixed layer thickness) 
#>               thetao (Temperature) 
#>               ...
#> names       : bottomT_1, bottomT_2, bottomT_3, bottomT_4, bottomT_5, bottomT_6, ... 
#> unit        : degrees_C, degrees_C, degrees_C, degrees_C, degrees_C, degrees_C, ... 
#> time        : 2018-08-12 to 2018-08-17 UTC

618 layers spread by variables, depth and time.

terra::time(r)[1:10]
#>  [1] "2018-08-12 UTC" "2018-08-13 UTC" "2018-08-14 UTC" "2018-08-15 UTC"
#>  [5] "2018-08-16 UTC" "2018-08-17 UTC" "2018-08-12 UTC" "2018-08-13 UTC"
#>  [9] "2018-08-14 UTC" "2018-08-15 UTC"

Let's convert the data frame to terra vect() object (I have assumed the CRS of .nc and your data.frame are the same, if they are different, use terra::project() function to mach them):

v <- terra::vect(data_locs, geom =c("Long", "Lat"), crs = terra::crs(r))

Now we have to: limit the rasters only to those which match the date, the depth and extract the values. As the depth in the layers doesn't match those given in data.frame either you would stick to neighbors, either you can use terra::approximate() function to estimate the values on dept 68 with surrounding layers.

r[[which(terra::time(r) == v[1, ]$date)]] |>
  terra::subset(subset = "thetao_depth=65.807266_1") |>
  terra::extract(v)
#>   ID thetao_depth=65.807266_1
#> 1  1                 12.42308
#> 2  2                 12.32127
#> 3  3                 12.36888

The values 12.423, 12.321 and 12.368 corresponds to 3 points from data.frame from depth = 65.8 (close to requested 68) for potential water temperature.

To get names of all layers use terra:names(r) and grepl() for requested variable/depth.

And the plot on the end because spatial analysis without plot is not an analysis :).

r[[which(terra::time(r) == v[1, ]$date)]] |>
  terra::subset(subset = "thetao_depth=65.807266_1") |>
  terra::plot()

^{Created on 2024-11-21 with reprex v2.1.1.9000}