I have a data.frame with three variables (ID, longitude, and latitude) for dolphin observations.
I also have 59 AQUA/MODIS netCDF files with sea surface temperature (SST). I would like to extract the raster data for SST for the dolphin observation locations.
Dolphin data
d <- read.csv("dolphins.csv")
str(d)
# 'data.frame': 650 obs. of 4 variables:
#$ ID : num 1 2 3 4 5 6 7 8 9 10 ...
#$ Date : Factor w/ 73 levels "1/24/17","1/24/18",..: 67 67 67 8 8 8 8 8 8 8 ...
#$ Latitude : num 42 42 42 42 42 ...
#$ Longitude: num 19.1 19.1 19.1 19.1 19.1 ...
SST rasters
library(terra)
filenames = list.files('Ocean_ColorSST_2016_2021',pattern='*.nc',full.names=TRUE)
SSTs <- rast(filenames, "sst")
SSTs
#class : SpatRaster
#dimensions : 4320, 8640, 59 (nrow, ncol, nlyr)
#resolution : 0.04166667, 0.04166667 (x, y)
#extent : -180, 180, -90.00001, 90 (xmin, xmax, ymin, ymax)
#coord. ref. : lon/lat WGS 84
#sources : AQUA_MODIS.20160901_20160930.L3m.MO.SST.sst.4km.nc:sst
# AQUA_MODIS.20161001_20161031.L3m.MO.SST.sst.4km.nc:sst
# AQUA_MODIS.20161101_20161130.L3m.MO.SST.sst.4km.nc:sst
# ... and 56 more source(s)
#varnames : sst (Sea Surface Temperature)
# sst (Sea Surface Temperature)
# sst (Sea Surface Temperature)
# ...
#names : sst, sst, sst, sst, sst, sst, ...
#unit : degree_C, degree_C, degree_C, degree_C, degree_C, degree_C, ...
Create a SpatialPointsDataFrame
library(sp)
points_spdf_M <- d
coordinates(points_spdf_M) <- ~ Latitude + Longitude
crs(points_spdf_M) <- crs(SSTs)
points_spdf_M
#class : SpatialPoints
#features : 650
#extent : 41.5978, 42.67778, 18.24337, 19.99933 (xmin, xmax, ymin, ymax)
#crs : +proj=longlat +datum=WGS84 +no_defs
Extract raster data from dolphin IDs returns NAs
library(raster)
ncin_SST <- stack(SSTs)
Extract_SST_M <- extract(ncin_SST, points_spdf_M)
head(Extract_SST_M)
#sst.1 sst.2 sst.3 sst.4 sst.5 sst.6 sst.7 sst.8 sst.9 sst.10 sst.11 sst.12 sst.13 sst.14 sst.15 sst.16 sst.17 sst.18
# [1,] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
# [2,] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
# [3,] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Another method that I tried
v <- vect(d, geom=c("Longitude", "Latitude"))
e <- extract(SSTs, v)
head(e)
Output ?
You create a SpatialPointsDataFrame like this
coordinates(Final_M_Points) <- ~ Latitude + Longitude
Where it should be
coordinates(Final_M_Points) <- ~ Longitude + Latitude
Your workflow simplified --- using only terra
library(terra)
filenames <- list.files("---")
SSTs <- rast(filenames, "sst")
d <- read.csv("file.csv")
v <- vect(d, geom=c("Longitude", "Latitude"))
e <- extract(SSTs, v)
It is better to use terra for this, and avoid sp and raster. In this case you would have been warned about the reversed lon and lat.