trajectory query example

[mdsumner]

Example data set provided at https://github.com/appelmar/gdalcubes_R/issues/26

## function to return data frame of 
## fullname (file path)
## date (date-time of the band in the file)
## band (number of the band within this file)
weather_files <- function(...) {
  ## the location of the files/s here is your job (as is caching this file list)
  files <- "weather_sample.nc"
  timefun <- function(file) RNetCDF::var.get.nc(RNetCDF::open.nc(file), "time")
  ## MDS I personally don't like to automate this step, the metadata is often
  ## insufficient
  #ncmeta::nc_att(files[1], "time", "units")$value$units
  ## [1] "hours since 1900-01-01 00:00:00.0"
  ltime <- lapply(files, timefun)
  lens <- lengths(ltime)  ## so we expand file name out for every band
  time <- unlist(ltime)
  data.frame(date = as.POSIXct("1900-01-01 00:00:00", tz = "UTC") + time * 3600, 
             fullname = rep(normalizePath(files), lens), 
             band = unlist(lapply(lens, seq_len)), stringsAsFactors = FALSE)
}

## function to read from netcdf file for a given date-time
## only finds the nearest, and will error if out of bounds
## 12hourly is a hardcoded string
read_weather <- function(date, varname= c("mwp", "swh"), ..., time.resolution = "12hourly", returnfiles = FALSE) {
  files <- weather_files()
  if (returnfiles) return(files)
  var <- match.arg(varname)
  if (missing(date)) date <- files$date[1]
  date <- as.POSIXct(date)
  stopifnot(length(date) == 1)  ## makes life easier atm
  idx <- which.min(abs(files$date - date))
  ## here we need to stop() if the date isn't within the bounds of the files
  if (date < min(files$date)) stop()
  if (date > max(files$date)) stop()
  raster::setZ(raster::setExtent(raster::raster(files$fullname[idx], varname = var, band = files$band[idx]), 
                    raster::extent(0, 360, -90, 90)), files$date[idx])
}

## read the track data
ais <- readr::read_csv("ais_sample.csv")
## query the read function with lon, lat, date-time
out <- raadtools::extract(read_weather, ais[c("LONGITUDE", "LATITUDE", "TIME")])

out is a vector of the mwp values at the right time and location (within the limits provided in the netcdf file)