0

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 ?

Robert Hijmans
  • 40,301
  • 4
  • 55
  • 63
Alice Hobbs
  • 1,021
  • 1
  • 15
  • 31
  • the dolphin data don't have time stamps? just asking as the AQUA are time based. you might also clip your SSTs to your dolphin sighting extent. – Chris Apr 10 '22 at 21:38
  • `terra::extract` has `layer` and `method` as arguments, you have 59 layers, and the error in the second approach points to `i`, meaning which one of those 59 should be extracted from. It could, dutifully, extract point in cell method 'simple' from each of 59, but that would likely make it difficult for you to predict or expect to find dophins where/when based on temp. On dog walks, I see dolphin heading north in May, last sightings south early September. But my dolphins aren't transmitting. – Chris Apr 10 '22 at 22:45
  • I added a timestamp column to the .csv file and I followed the steps that I have shown above and the output is still 59 columns of NAs. I feel very confused as I am a complete novice to this type of analysis – Alice Hobbs Apr 11 '22 at 04:38
  • Me too. But I find, if I simplify a problem to just one case, a sighting say, at the same time as AQUA[1] or [42], I have a shot at generalizing to all sightings against the stack. Still think it useful to crop your SSTs to sighting extents as there's an awful lot of 'dolphin not there' data to lug around. And think of this as play. – Chris Apr 11 '22 at 04:45
  • I don't usually play with factors and you may find dolphin sightings as factor troublesome, though I can't say why, except I was happy when stringsAsFactors = FALSE became default. Probably easier, overall, that your dates match format, but, ok, your cols are 5 years of months...head explodes. – Chris Apr 11 '22 at 04:55
  • Thanks, Chris, I really do appreciate your advice. I understand what you mean when you say select a sighting that is at the same time as the AQUA, but I don't understand what you mean by I have a shot at generalizing to all sightings against the stack. I am a complete novice, so could you please explain further? Thanks for your patience – Alice Hobbs Apr 11 '22 at 04:56
  • My head is exploding too. I managed to extract the data that I needed from another shapefile using this method, with no problem. For some reason, it is not working for this shapefile and I don't understand what the difference is, apart from the other shapefile had 1 feature (with 650 ID's), and the shapefile in the example has 650 features, which are my ID's. – Alice Hobbs Apr 11 '22 at 04:58
  • So reorganize the question a little. `dput(head(your_cvs_object))` and paste in `structure(...)`, then scratch about for SSTs that bracket the 6 observations and perhaps have link to them...then decide if a dolphin sample belongs to the temps of the SST prior, or after. But, probably first put, 'before, I did this and everything went swimmingly. And now, ...' And then someone who really knows what they are doing will come along and speak intelligently resolving this. That way you'll have a complete MRE. – Chris Apr 11 '22 at 05:09
  • You're down to small differences and looking at it till you're blind, but your dolphin $ values are upper case `L`, not lower `l`, details I frustratingly often overlook.. – Chris Apr 12 '22 at 12:48
  • And if you`dput(head(d))`, we can track down our own AQUA. – Chris Apr 12 '22 at 16:34

1 Answers1

2

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.

Robert Hijmans
  • 40,301
  • 4
  • 55
  • 63
  • Thank you, Robert, I was really happy about your solution. I am starting to think that the problem might be in the AQUA/MODIS files themselves and not the code. I ran an is.na analysis and all the columns stated 'FALSE'. After spending a week configuring why, I then zoomed into where my administration boundary file, the point data, and the AQUA/MODIS files intersect in QGIS, and I realised that the AQUA/MODIS files do not cover the inner coastal areas. Therefore, a quarter of my points are not overlaying onto the raster files, so it's impossible to extract the data. Hence the NAs. – Alice Hobbs Apr 12 '22 at 11:55
  • Please see the output from the example code that you provided above. I am not an expert, but I think the NULL values for the points that do not overlay on the SST raster might be affecting the output. Do you think this assumption is correct? – Alice Hobbs Apr 12 '22 at 11:58
  • I have edited your question for clarity. There was a lot of mix-up making it very difficult to follow what you actually do (as yo try many things). Have a look at that, and again at my answer. The best course of action is to *not* use sp (no SpatialPoints) and raster (no stack); only terra (`rast`, `vect` and `extract`). – Robert Hijmans Apr 12 '22 at 16:03