Extract time series subset from NetCDF for many spatial points using R [closed]

Question

I have been following mdsumner's response on Extract time series from .nc files for many locations simultaneously which is a similar question, but am having some trouble. Here is my code:

library(raster)
b <- brick("http://thredds.northwestknowledge.net:8080/thredds/ncss/grid/MACAV2/bcc-csm1-1/macav2metdata_tasmax_bcc-csm1-1_r1i1p1_rcp85_2036_2040_CONUS_daily.nc?&var=air_temperature&north=47.6200&south=45.9300&west=-121.6100&east=-119.0900&temporal=all&accept=netcdf&point=false",
           varname = "air_temperature")

pts <- read.csv("/predpts_latlong.csv",
                stringsAsFactors = FALSE)
# The above file path was removed for privacy - file linked below

ts <- extract(b, cbind(pts$LONGDD, pts$LATDD))

The result is a table with all of my lat/long points filled with null values.

Here is a link to my predpts_latlong.csv file: https://drive.google.com/file/d/1Gc2djozljTbYIR5XvcMgDZ1lJfjMTWG5/view?usp=drive_link

A few other questions: Is there a way to select a specific subset of the days within the NetCDF file? I'm only looking for four weeks out of every year in the NetCDF. Specifically, June 10th - July 7th for each year.

I have unique IDs for every lat/long point in my predpts_latlong.csv (not currently in the file). How can I add these onto the resulting data table?

Also, I need to extract this data from about 100 netCDFs. Is there a way to streamline this process for multiple netCDF files?

Answer 1

I would recommend using the terra library which is the modern analog to raster. Here is a simple workflow that will read your data, coerce to spatial classes and extract the raster values at your points.

Add libraries and read NetCDF raster data.

library(terra)
library(sf)

b <- rast("http://thredds.northwestknowledge.net:8080/thredds/ncss/grid/MACAV2/bcc-csm1-1/macav2metdata_tasmax_bcc-csm1-1_r1i1p1_rcp85_2036_2040_CONUS_daily.nc?&var=air_temperature&north=47.6200&south=45.9300&west=-121.6100&east=-119.0900&temporal=all&accept=netcdf&point=false")

I do not have access to your points so, I just create some dummy data that would look like your csv file. The coordinate columns here are "x","y" so rename accordingly. After reading the flat file, I coerce into a spatial sf POINT object. If your point data is not in the same coordinate system as your raster, then things will not match and NA's (NaN) will be produced.

Note; when plotting the sf object I call st_geometry to only plot the geometry with no attributes. The index when plotting the raster is in reference to a specific layer so, b[[1:10]] would be the first 10 layers.

pts <- as.data.frame(spatSample(b[[1]], 50, xy=TRUE, values=FALSE)) 
  pts$ID <- 1:nrow(pts)
    head(pts)

pts <- st_as_sf(pts,  coords = c("x", "y"), crs = st_crs(crs(b)), 
                agr = "constant")

plot(b[[1]])
  plot(st_geometry(pts), pch=20, col="black", add=TRUE)

Now, you can extract your raster values. This will pull layer-by-layer values so the number of columns in the resulting matrix is equal to the number of layers (plus an ID column) so, in this case 1873. Once you have the data extracted you can cbind it to your points. Since geometry is sticky in sf objects, you do not need to worry about it becoming a different object class, which is normally the case.

bdat <- extract(b, pts)
pts <- cbind(pts, bdat[,-1]) 

plot(pts["air_temperature_1"], pch=20, cex=2)

As far as indexing by date, it would have been helpful if you had provided some information regarding the temporal period (year range) and sampling interval (ie., daily, monthly, yearly). Thankfully, it looks like it is encoded in the ncdf file (2036-01-01 to 2040-12-31). There is a function in the spatialEco package that will create a date object representing a sequence of from-to by time step.

Here we create a dates object of 2036-01-01 to 2040-12-31 by day (n=1827).

library(spatialEco)   

dates <- spatialEco::date_seq("2036-01-01", "2040-12-31", step = "day")

A quick way to query would be to create an index by day of year or week. For the dates you specified it would be weeks 24-27 (week 24 - June 10, 2024 through June 16, 2024, ..., week 27 - July 1, 2024 through July 7, 2024). We can coerce our daily dates to weeks or day of year (for each daily) using strftime and then use which to create the index query. Then we simply pass the index to the raster object using a double bracket.

doy <- strftime(dates, format = "%j")
  ( didx <- which(doy >= 161 & doy <= 188) )

weeks <- strftime(dates, format = "%V")
  ( widx <- which(weeks %in% c(24:27)) )

Before sub-setting your raster layers, you can verify the index by using it on the weeks or dates vector.

weeks[didx]
dates[didx]

b.sub <- b[[didx]]

If we wanted to create a more specific query (say slightly different growing season for each year) we can create a series of yearly queries and concatenate them together. Here is an example that replicates the above index.

( idx <- c(which(dates %in% date_seq("2036-04-10", "2036-05-7", step = "day")),
  which(dates %in% date_seq("2037-04-10", "2037-05-7", step = "day")),
  which(dates %in% date_seq("2038-04-10", "2038-05-7", step = "day")),
  which(dates %in% date_seq("2039-04-10", "2039-05-7", step = "day")),
  which(dates %in% date_seq("2040-04-10", "2040-05-7", step = "day"))) )

To process multiple files, look at for or lapply. Create a vector of your URL's, and iterate through them writing the resulting spatial objects, storing them in a list or appending them.