copernicus altimetry changes complete, now readcurr() is compatible w…

…ith 'read_<alt>_daily' family of functions

NAMESPACE

@@ -34,11 +34,13 @@ export(read_amsr_ice)
 export(read_amsre_ice)
 export(read_ccmp)
 export(read_cersat_ice)
+export(read_dir_daily)
 export(read_err_daily)
 export(read_geoid)
 export(read_leads_clim)
 export(read_leads_clim_north)
 export(read_leads_clim_south)
+export(read_mag_daily)
 export(read_oc_sochla)
 export(read_rema_tiles)
 export(read_sla_daily)

NEWS.md

@@ -1,8 +1,17 @@
 # raadtools dev
 
+* Converting read and files functions to versions that use GDAL more tightly, in vapour there is a function `gdal_raster_data()` which 
+can be given an extent, dimension, resolution, crs, or combination of these (or none of them) and will return the data read by the warper api. The data returned has attributes extent, dimension, projection (crs) so we know what the properties of the data are (so a user can give
+partial or missing specifications of what 'xylim' is (an extent in the native projection this function uses), and we can maintain the behaviour of raadtools up to this point. But, it's faster and simpler.)  Improvements include no longer needing dateline logic, so we can give c(0, 220, -60, 10) or any other longlat extent and gdal sorts it out. raadtools is currently not using any reprojection so this only applies for longlat sources. 
+
+The changes involve using 'vrt_dsn' which is an augmented version of 'fullname', and makes it useable in other contexts. Sometimes we have combinations like "ugos_vrt_dsn" and "vgos_vrt_dsn", but we're working to remove those inconsistencies and simply provide 1 variable for 1 or many time steps as a new basis. 
+
+So far updated are readsst, readcurr and its siblings read_ugos_daily, read_vgos_daily, read_adt_daily, read_err_daily, read_ugosa_daily, read_vgosa_daily, read_sla_daily, and derived versions read_mag_daily and read_dir_daily (trigonometric combinations of read_ugos and read_vgos). This trigonometry on u/v is part of the reason we are avoiding reprojected output, that will require extra thought because currently (haha) all vector quantities are assumed aligned to their native grid and we aren't messing with that. 
+
+
 * Removed rgeos and maptools. 
 
-* Removed rgdal imports, replaced with raster and reproj. 
+* Removed rgdal imports, replaced with raster and reproj.  This has meant there were cases where we were using SpatialPoints and now we just have a matrix x,y - so there were cases where we needed 'drop = FALSE' for when one point was looked up (and R was flattening to a two element vector, interpreted as cell number - in extract(function, matrix) for example). 
 
 
 *Function `readice()` has been expanded to allow 'hemisphere = "both"' and for 'xylim' to be a full raster grid (terra or raster format). If 'both' is specified the warper is applied to VRT versions of the NSIDC files, which allows them to be combined in one reprojection step. In this case 'xylim' can be specified, to give a projected grid of any form. If not supplied (when hemisphere = 'both') then longlat raster at 0.25 degrees is assumed. 

R/copernicus_altimetry.R

@@ -81,7 +81,24 @@ read_err_daily <- function(date, xylim = NULL, latest = TRUE, returnfiles = FALS
   read_copernicus_daily(date = date, xylim = xylim, latest  = latest, returnfiles = returnfiles, varname = varname, ..., inputfiles = inputfiles)
 }
 
+vlen <- function(x, y) sqrt(x * x + y * y)
+#' @name read_adt_daily
+#' @export
+read_mag_daily <- function(date, xylim = NULL, latest = TRUE, returnfiles = FALSE, ..., inputfiles = NULL) {
+  u <- read_ugos_daily(date = date, xylim = xylim, latest  = latest, returnfiles = returnfiles,  ..., inputfiles = inputfiles)
+  if (returnfiles) return(u)
+  v <- read_vgos_daily(date = date, xylim = xylim, latest  = latest, returnfiles = returnfiles,  ..., inputfiles = inputfiles)
+  vlen(u, v)
+}
 
+#' @name read_adt_daily
+#' @export
+read_dir_daily <- function(date, xylim = NULL, latest = TRUE, returnfiles = FALSE, ..., inputfiles = NULL) {
+  u <- read_ugos_daily(date = date, xylim = xylim, latest  = latest, returnfiles = returnfiles,  ..., inputfiles = inputfiles)
+  if (returnfiles) return(u)
+  v <- read_vgos_daily(date = date, xylim = xylim, latest  = latest, returnfiles = returnfiles,  ..., inputfiles = inputfiles)
+  (90 - atan2(v, u) * 180/pi) %% 360
+}
 
 read_copernicus_daily <- function(date, xylim = NULL, latest = TRUE, returnfiles = FALSE, varname, lon180 = FALSE, ..., inputfiles = NULL) {
   if (is.null(inputfiles)){

R/currents.R

@@ -77,8 +77,8 @@ currentsfiles <- function(time.resolution = c("daily", "weekly"), ...) {
   time.resolution <- match.arg(time.resolution)
   if (time.resolution != "daily") warning("only daily available, no weekly - ignoring 'time.resolution'")
   #raadfiles::altimetry_daily_files()
-  files <- dplyr::inner_join(dplyr::rename(altimetry_daily_ugos_files(), ugos_vrt = vrt_dsn), 
-                             altimetry_daily_vgos_files() |> dplyr::transmute(date, vgos_vrt = vrt_dsn), "date")
+  files <- dplyr::inner_join(dplyr::rename(altimetry_daily_varname_files("ugos"), ugos_vrt = vrt_dsn), 
+                             altimetry_daily_varname_files("vgos") |> dplyr::transmute(date, vgos_vrt = vrt_dsn), "date")
   files
 }
 
@@ -187,13 +187,11 @@ readcurr <- function (date, time.resolution = c("daily"),
     files <- inputfiles
   }
 
-
-  thefun <- read_i
 
-  if (magonly) thefun <- read_i_mag0
-  if (dironly) thefun <- read_i_dir0
-  if (uonly ) thefun <- read_i_u0
-  if (vonly) thefun <- read_i_v0
+  if (magonly) thefun <- read_mag_daily
+  if (dironly) thefun <- read_dir_daily
+  if (uonly ) thefun <- read_ugos_daily
+  if (vonly) thefun <- read_vgos_daily
 
 
 
@@ -206,27 +204,29 @@ readcurr <- function (date, time.resolution = c("daily"),
   nfiles <- nrow(files)
 
   ## prevent reading more than one unless mag/dironly
-  if (nfiles > 1L & !magonly & !dironly & !uonly & !vonly) {
+  if ( !magonly & !dironly & !uonly & !vonly) {
     files <- files[1L,]
+    message("reading u and v as layers is soft deprecated, please use read_ugos_daily() or read_ugos_daily()")
+    message("for reading mag-nitude or dir-ection only, please use read_mag_daily() or read_dir_daily()")
+
     nfiles <- 1L
-    warning("only one time step can be read at once unless one of 'magonly', 'dironly', 'uonly' or 'vonly' is TRUE")
+
+    thefun <- function(date, xylim = NULL, latest = TRUE, returnfiles = FALSE, ..., inputfiles = NULL) {
+      u <- read_ugos_daily(date = date, xylim = xylim, latest  = latest, returnfiles = returnfiles,  ..., inputfiles = inputfiles)
+      if (returnfiles) return(u)
+      v <- read_vgos_daily(date = date, xylim = xylim, latest  = latest, returnfiles = returnfiles,  ..., inputfiles = inputfiles)
+      stack(list(u, v))
+    }
   }
   if ((magonly + dironly + uonly + vonly) > 1) stop("only one of 'magonly', 'dironly', 'uonly' or 'vonly' may be TRUE")
 
-
-
-  print(files$date)
-  m <- vlen(read_i(files$ugos_vrt, grid), read_i(files$vgos_vrt, grid))
-
-  #ximage::ximage(matrix(m, grid$dimension[2], byrow = TRUE))
-  return(NULL)
-  dots <- list(...)
 
 
-
   op <- options(warn = -1)
   on.exit(options(op))
-  r0 <- stack(lapply(files$fullname, thefun, xylim = xylim, lon180 = lon180), filename = filename)
+
+  r0 <- lapply(files$date, thefun, xylim = xylim, lon180 = lon180)
+  if (length(r0) == 1) r0 <- r0[[1]] else r0 <- stack(r0)
   if (nlayers(r0) == nrow(files)) {
     r0 <- setZ(r0, files$date)
   } else {
@@ -246,112 +246,4 @@ readcurr <- function (date, time.resolution = c("daily"),
 r0
 
 
-}
-
-# dimensions:
-#   time = 1 ;
-# lat = 720 ;
-# lon = 1440 ;
-# nv = 2 ;
-# variables:
-#   float time(time) ;
-# time:long_name = "Time" ;
-# time:standard_name = "time" ;
-# time:units = "days since 1950-01-01 00:00:00 UTC" ;
-# time:calendar = "julian" ;
-# time:axis = "T" ;
-# float lat(lat) ;
-# lat:long_name = "Latitude" ;
-# lat:standard_name = "latitude" ;
-# lat:units = "degrees_north" ;
-# lat:bounds = "lat_bnds" ;
-# lat:axis = "Y" ;
-# lat:valid_min = -90 ;
-# lat:valid_max = 90 ;
-# float lat_bnds(nv, lat) ;
-# float lon(lon) ;
-# lon:long_name = "Longitude" ;
-# lon:standard_name = "longitude" ;
-# lon:units = "degrees_east" ;
-# lon:bounds = "lon_bnds" ;
-# lon:axis = "X" ;
-# lon:valid_min = 0 ;
-# lon:valid_max = 360 ;
-# float lon_bnds(nv, lon) ;
-# int crs ;
-# crs:grid_mapping_name = "latitude_longitude" ;
-# crs:semi_major_axis = 6371000 ;
-# crs:inverse_flattening = 0 ;
-# int nv(nv) ;
-# int u(lon, lat, time) ;
-# u:_FillValue = -2147483647 ;
-# u:long_name = "Absolute geostrophic velocity: zonal component" ;
-# u:standard_name = "surface_eastward_geostrophic_sea_water_velocity" ;
-# u:units = "m/s" ;
-# u:scale_factor = 1e-04 ;
-# int v(lon, lat, time) ;
-# v:_FillValue = -2147483647 ;
-# v:long_name = "Absolute geostrophic velocity: meridian component" ;
-# v:standard_name = "surface_northward_geostrophic_sea_water_velocity" ;
-# v:units = "m/s" ;
-# v:scale_factor = 1e-04 ;
-# 
-# // global attributes:
-#   :cdm_data_type = "Grid" ;
-# :title = "DT merged Global Ocean Gridded Absolute Geostrophic Velocities SSALTO/Duacs L4 product" ;
-# :summary = "This dataset contains Delayed Time Level-4 absolute geostrophic velocities products from multi-satellite observations over Global Ocean." ;
-# :comment = "Surface product; Absolute Geostrophic Velocities" ;
-# :time_coverage_resolution = "P1D" ;
-# :product_version = "5.0" ;
-# :institution = "CNES, CLS" ;
-# :project = "SSALTO/DUACS" ;
-# :references = "www.aviso.altimetry.fr" ;
-# :contact = "aviso@altimetry.fr" ;
-# :license = "http://www.aviso.altimetry.fr/fileadmin/documents/data/License_Aviso.pdf" ;
-# :platform = "ERS-1, Topex/Poseidon" ;
-# :date_created = "2014-02-28 13:07:00" ;
-# :history = "2014-02-28 13:07:00:creation" ;
-# :Conventions = "CF-1.6" ;
-# :standard_name_vocabulary = "http://cf-pcmdi.llnl.gov/documents/cf-standard-names/standard-name-table/12/cf-standard-name-table.html" ;
-# :geospatial_lat_min = -90 ;
-# :geospatial_lat_max = 90 ;
-# :geospatial_lon_min = 0 ;
-# :geospatial_lon_max = 360 ;
-# :geospatial_vertical_min = "0.0" ;
-# :geospatial_vertical_max = "0.0" ;
-# :geospatial_lat_units = "degrees_north" ;
-# :geospatial_lon_units = "degrees_east" ;
-# :geospatial_lat_resolution = 0.25 ;
-# :geospatial_lon_resolution = 0.25 ;
-
-
-
-.currentsfiles1 <- function(fromCache = TRUE, ...) {
-  # datadir = getOption("default.datadir")
-  # cachefile <- file.path(datadir, "cache", sprintf("currentsfiles_weekly.Rdata"))
-  # if (fromCache) {
-  #   load(cachefile)
-  #   cfs$fullname <- file.path(datadir, cfs$file)
-  #   return(cfs)
-  # }
-  # 
-  ftx <- .allfilelist()
-  cfiles <- grep("aviso_old", ftx, value = TRUE)
-  cfiles1 <- grep("current", cfiles, value = TRUE)
-  cfiles2 <- grep("merged_madt", cfiles1, value = TRUE)
-  cfiles3 <- grep("nc$", cfiles2, value = TRUE)
-  #data.source = file.path(datadir, "current", "aviso", "upd", "7d")
-  #cfiles <- list.files(data.source, pattern = ".nc$", full.names = TRUE)
-  datepart <- sapply(strsplit(basename(cfiles3), "_"), function(x) x[length(x)-1])
-  currentdates <- timedateFrom(as.Date(strptime(datepart, "%Y%m%d")))
-
-  cfs <- data.frame(fullname = cfiles3,  date = currentdates, stringsAsFactors = FALSE)
-  cfs <- cfs[diff(cfs$date) > 0, ]
-
-  ## drop duplicates, this should prefer upd to nrt
-  cfs <- cfs[!duplicated(cfs$date), ]
-  #save(cfs, file = cachefile)
-  #cfs$fullname <- file.path(datadir, cfs$file)
-  cfs
-
-}
+}

R/extract.R

@@ -54,15 +54,7 @@ setOldClass("trip")
     pb$tick(0) ## ---------------------------------------------
     ## TODO, will have to figure out how to do this
     args <- list(...)
-    # nousexylim <- FALSE
-    # if ("xylim" %in% names(args)) {
-    #   if (inherits(args$xylim, "SpatRaster") ||  inherits(args$xylim, "BasicRaster")) {
-    #    nousexylim <- TRUE 
-    #   } else {
-    #     warning("xylim argument ignored (determined automatically from the input data)")
-    #     args$xylim <- NULL
-    #   }
-    # }
+
     if ("inputfiles" %in% names(args)) {
       warning("inputfiles argument ignored")
       args$inputfiles <- NULL
@@ -87,15 +79,7 @@ setOldClass("trip")
     dummy <- x(inputfiles = files, ...)
     yp <- reproj::reproj_xy(y1, projection(dummy), source = "+proj=longlat")
     pb$tick(0) ## ---------------------------------------------
-    # xylim <- extent(yp)
-    # ## expand out a bit for single-location queries
-    # if (xmax(xylim) == xmin(xylim) | ymax(xylim) == ymin(xylim)) {
-    #   xylim <- xylim + res(dummy) * 3
-    # }
-    # 
-    # ## never crop
-    # xylim <- NULL
-    # 
+
     if (notime) {
       ## assume we want topo/bathy values
       thisx1 <- x(...)
@@ -187,7 +171,7 @@ setOldClass("trip")
         if(resize) thisx <- aggregate(thisx, fact = fact, fun = "mean")
         asub <- windex == findex[i]
         ## no interpolation in time, controlled by "method" for xy
-        if (any(asub)) {result[asub] <- suppressWarnings(extract(thisx, y[asub, ], ...))}
+        if (any(asub)) {result[asub] <- suppressWarnings(v000 <- extract(thisx, y[asub, , drop = FALSE], ...))}
 
         if (interactive() & verbose) {
           pb$tick(tokens = list(what = time.resolution, ith = i, nn = length(date)))