How to Extract Point Data from Multiple Subsetted Netcdf Files for Sea Surface Density using terra::rast() and terra::extract functions in R

Description of Data

I downloaded MULTIOBS_GLO_PHY_S_SURFACE_MYNRT_015_013 data for daily sea surface density from 2012 to 2024 from The Copernicus Website. A person can only download 2 years of data at a time; therefore, I have five .nc files ranging over 12 years.

Copernicus Website

https://data.marine.copernicus.eu/product/MULTIOBS_GLO_PHY_S_SURFACE_MYNRT_015_013/download

Desired Output

My aim is to extract the point data from the netcdf files from 2012 to 2024 for the specific survey dates that we did during fieldwork and insert the values as a column in my dataframe.

Issue

I get this error message when I run the line SSSs to select the variable sea surface density from the netcdf files

Error: Not compatible with requested type: [type=character; target=double].

I have successfully extracted chlorophyll-a and sea surface temperature point data from Lm3 AQUA MODIS files using exactly the same method and I don't understand what I'm doing wrong with sea surface salinity.

I can't share the dataframe due to ownership issues

Many thanks if anyone can help

Code

Read in all 5 netcdf files for sea surface density

#SS1: January 1st 2012 to December 31st 2014

sss1<-terra::rast("~/Documents/cmems_obs-mob_glo_phy-sss_my_multi_P1D_1729236771868.nc")
plot(sss1)
sss1<-subset(sss1, 1:1096)

#SS2: January 1st 2015 to December 31st 2017

ss2<-terra::rast("~/Documents/cmems_obs-mob_glo_phy-sss_my_multi_P1D_1729236939325.nc")

#Subset to extract sea surface density only 
ss2<-subset(ss2, 1:1096)

#SS3: January 1st 2018 to December 31st 2020

ss3<-terra::rast("~/Documents/cmems_obs-mob_glo_phy-sss_my_multi_P1D_1729237129644.nc")

#Subset to extract sea surface density only 
ss3<-subset(ss3, 1:1096)

#SS4: January 1st 2021 to December 16th 2023

ss4<-terra::rast("~/Documents/cmems_obs-mob_glo_phy-sss_my_multi_P1D_1729252813869.nc")

#Subset to extract sea surface density only 
ss4<-subset(ss3, 1:1080)


#SS5: December 17th 2023 to October 12th 2024

ss5<-terra::rast("~/Documents/cmems_obs-mob_glo_phy-sss_nrt_multi_P1D_1729253007917.nc")

#Subset to extract sea surface density only 
ss5<-subset(ss3, 1:301)

Concatenate all the subsetted netcdf files from Jan 2012 to Oct 12th 2024 together

total_sss<-c(sss1, ss2, ss3, ss4, ss5)

Extract the Point Data

#Select the sea surface density variable from the netcdf files 
SSSs <- terra::rast(total_sss, "dos")

#Make the dataframe a spatial object of class = "sf" with a CRS of 4326
Ds_Points <- st_as_sf(x=MyDf, 
                      coords = c("Longitude_E_DD", "Latitude_N_DD"), 
                      crs = 4326)
    
#Extract the SSS values from the netcdf files 
SSSs_Data <- terra::extract(SSSs, Ds_Points)

#Add the SSS results to the dataframe with the function cbind 
MyDf <- cbind(MyDf, SSS = SSSs_Data$dos)

This is the information stored in 'total_sss'

class       : SpatRaster 
dimensions  : 189, 154, 4669  (nrow, ncol, nlyr)
resolution  : 0.125, 0.125  (x, y)
extent      : 28.75, 48, 8.625, 32.25  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 
sources     : cmems_obs-mob_glo_phy-sss_my_multi_P1D_1729236771868.nc:dos  (1096 layers) 
              cmems_obs-mob_glo_phy-sss_my_multi_P1D_1729236939325.nc:dos  (1096 layers) 
              cmems_obs-mob_glo_phy-sss_my_multi_P1D_1729237129644.nc:dos  (1096 layers) 
              ... and 2 more source(s)
varnames    : dos (Sea surface density) 
              dos (Sea surface density) 
              dos (Sea surface density) 
              ...
names       : dos_d~=-0_1, dos_d~=-0_2, dos_d~=-0_3, dos_d~=-0_4, dos_d~=-0_5, dos_d~=-0_6, ... 
unit        :       kg/m3,       kg/m3,       kg/m3,       kg/m3,       kg/m3,       kg/m3, ... 
time        : 2012-01-01 to 2020-12-31 UTC

This is the information stored in SS1

class       : SpatRaster 
dimensions  : 189, 154, 1096  (nrow, ncol, nlyr)
resolution  : 0.125, 0.125  (x, y)
extent      : 28.75, 48, 8.625, 32.25  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 
source      : cmems_obs-mob_glo_phy-sss_my_multi_P1D_1729236771868.nc:dos 
varname     : dos (Sea surface density) 
names       : dos_d~=-0_1, dos_d~=-0_2, dos_d~=-0_3, dos_d~=-0_4, dos_d~=-0_5, dos_d~=-0_6, ... 
unit        :       kg/m3,       kg/m3,       kg/m3,       kg/m3,       kg/m3,       kg/m3, ... 
time        : 2012-01-01 to 2014-12-31 UTC

Solution

Replace:

SSSs <- terra::rast(total_sss, "dos")

with

SSSs <- terra::subset(total_sss, 1:xx)

where xx is number of layers for dos varname or with

SSSs <- total_sss[[1:xx]]

SSSs <- total_sss["dos_d"]

where dos_d is partial mach of names. Please have a look on terra:subset() function. And Robert's answer: How to properly subset data using their variable names in terra?

BTW, to temporary share files you can use google drive/github.