Abbreviated "time" specification, documentation update

Author: pvanlaake

DESCRIPTION

@@ -63,6 +63,7 @@ Suggests:
     ggplot2,
     knitr,
     rmarkdown,
+    stars,
     terra,
     testthat (>= 3.0.0),
     units 

NEWS.md

@@ -2,11 +2,11 @@
 
 #### API
 - New `geom_ncdf()` function to create a `geom` for use in map composition using package `ggplot2`.
-- In `CFVariable$subset()` and `$profile()` subsetting a "time" axis can now use abbreviated specification such as `time = "2025-12"` to select all data for December 2025, or `time = c("2020-S1", "2025-S4")` for all meteorological season data from 2019-12-01 to 2025-11-30 (inclusive). Abbreviation can be by month, meteorological season (S1 to S4), quarter (Q1 to Q4) or dekad (D01 to D36). 
+- In `CFVariable$subset()` subsetting a "time" axis can now use abbreviated specification such as `time = "2025-12"` to select all data for December 2025, or `time = c("2020-S1", "2025-S4")` for all meteorological season data from 2019-12-01 to 2025-11-30 (inclusive). Abbreviation can be by year, month, meteorological season (S1 to S4), quarter (Q1 to Q4) or dekad (D01 to D36). 
 
 #### Code
 - Fixed writing attributes for groups.
-- Fixed reading of scalar variables following saving of `CFDataset`.
+- Fixed reading of scalar variables following subsetting of a scalar axis.
 - Documentation updates.
 
 # ncdfCF 0.8.0

R/CFAxisTime.R

@@ -63,7 +63,8 @@ CFAxisTime <- R6::R6Class("CFAxisTime",
         if (ym$month[2L] == 12L) {
           ym$year[2L] <- ym$year[2L] + 1L
           ym$month[2L] <- 1L
-        }
+        } else
+          ym$month[2L] <- ym$month[2L] + 1L
         return(sprintf("%04d-%02d-01", ym$year, ym$month))
       }
 
@@ -111,7 +112,7 @@ CFAxisTime <- R6::R6Class("CFAxisTime",
       if (ts[1L] == ts[2L]) {
         ymd <- private$.tm$calendar$parse(ts[1L])
         if (is.na(ymd$year))
-          stop("Bad format for timestamps.", call. = FALSE)
+          stop("Bad format for timestamps: Date not valid calendar.", call. = FALSE)
         ymd <- rbind(ymd, private$.tm$calendar$add_day(ymd))
         return(sprintf("%04d-%02d-%02d", ymd$year, ymd$month, ymd$day))
       }
@@ -360,12 +361,13 @@ CFAxisTime <- R6::R6Class("CFAxisTime",
     #' @param rightmost.closed Whether or not to include the upper limit.
     #'   Default is `FALSE`.
     #' @return An integer vector giving the indices in the time axis between
-    #'   values in `x`, or `integer(0)` if none of the values are valid.
+    #'   values in `x`, or `NULL` if none of the values are valid.
     slice = function(x, rightmost.closed = FALSE) {
       x <- private$expand_timestamps(x)
       time <- private$.tm
-      idx <- time$slice(x, rightmost.closed)
-      range((1L:length(time))[idx])
+      idx <- suppressWarnings(time$slice(x, rightmost.closed))
+      if (all(!idx)) NULL
+      else range((1L:length(time))[idx])
     },
 
     #' @description Return an axis spanning a smaller coordinate range. This

README.Rmd

@@ -35,141 +35,19 @@ fn <- system.file("extdata", "ERA5land_Rwanda_20160101.nc", package = "ncdfCF")
 
 # Open the file, all metadata is read
 (ds <- open_ncdf(fn))
-
-# ...or very brief details
-ds$var_names
-ds$axis_names
-
-# Variables and axes can be accessed through standard list-type extraction syntax
-(t2m <- ds[["t2m"]])
-
-ds[["longitude"]]
-
-# Regular base R operations simplify life further
-dimnames(ds[["pev"]]) # A variable: list of axis names
-dimnames(ds[["longitude"]]) # An axis: vector of axis coordinate values
-
-# Access attributes
-ds[["pev"]]$attribute("long_name")
-```
-
-If you just want to inspect what data is included in the netCDF resource, use the `peek_ncdf()` function: 
-
-```{r peek}
-peek_ncdf(fn)
-```
-
-## Extracting data
-
-There are various ways to read data for a data variable from the resource: 
-
- * **`[]`:** The usual R array operator gives you access to the raw, non-interpreted data in the netCDF resource. This uses index values into the dimensions and requires you to know the order in which the dimensions are specified for the variable. With a bit of tinkering and some helper functions in `ncdfCF` this is still very easy to do.
- * **`raw()`:** This also gets the data in the layout of the file (or the data set) but with dimnames set. Importantly, you can call this after calling `subset()` and you will get the raw data for the specific spatial and temporal domain that you are interested in.
- * **`array()`:** Like `raw()`, this extracts all the (subsetted) data, but now the data will be oriented in the standard R way of column-major order. Y coordinates will run from the top to the bottom (so latitude values, for instance, will be decreasing).
- * **`subset()`:** The `subset()` method lets you specify what you want to extract from each dimension in real-world coordinates and timestamps, in whichever order. This can also rectify non-Cartesian grids to regular longitude-latitude grids. Subsetting is lazy: data is not loaded so long as a direct relationship to the data in the netCDF resource is maintained.
- * **`profile()`:** Extract "profiles" from the data variable. This can take different forms, such as a temporal or depth profile for a single location, but it could also be a zonal field (such as a transect in latitude - atmospheric depth for a given longitude) or some other profile in the physical space of the data variable.
- 
-```{r extract}
-# Extract a timeseries for a specific location - see also the `profile()` method
-ts <- t2m[5, 4, ]
-str(ts)
-
-# Extract the full spatial extent for one time step
-ts <- t2m[, , 12]
-str(ts)
 ```
 
-Note that the results contain degenerate dimensions (of length 1). This by design when using basic `[]` data access because it allows attributes to be attached in a consistent manner. When using the `subset()` method, the data is returned as an instance of `CFVariable`, including axes and attributes:
-
-```{r subset}
-# Extract a specific region, full time dimension
-(ts <- t2m$subset(list(X = 29:30, Y = -1:-2)))
-
-# Extract specific time slices for a specific region
-# Note that the dimensions are specified out of order and using alternative
-# specifications: only the extreme values are used.
-(ts <- t2m$subset(list(T = c("2016-01-01 09:00", "2016-01-01 15:00"),
-                       X = c(29.6, 28.8),
-                       Y = seq(-2, -1, by = 0.05))))
-```
-
-Data loading is lazy. In the examples above, you can see that data did not yet get loaded. This is intentional: you can subset your data in multiple ways before actually reading the data from the resource. This is particularly important when getting data from an online location, such as a remote THREDDS server. Use `raw()` or `array()` to get the arrays.
-
-##### Make a profile of data
-
-It is often useful to extract a "profile" of data for a given location or zone, such as a timeseries of data. The `profile()` method has some flexible options to support this:
-
- - Profile specific locations, with multiple locations specified per call, returning the data as a (set of) `CFVariable` instance(s) or as a single `data.table`.
- - Profile zones, such as a latitude band or an atmospheric level. Data is returned as a new `CFVariable` instance(s).
- 
-In all cases, you can profile over any of the axes and over any number of axes.
-
-Note that the `profile()` method returns data for the grid cells closest to the specified location. That is different from the `subset()` method, which will return data as it is recorded in the netCDF resource.
-
-```{r profile}
-rwa <- t2m$profile(longitude = c(30.07, 30.07, 29.74), latitude = c(-1.94, -1.58, -2.60), 
-                   .names = c("Kigali", "Byumba", "Butare"), .as_table = TRUE)
-head(rwa)
-attr(rwa, "value")
-```
-
-Some critical metadata is recorded in the "value" attribute: original long name and the physical unit.
-
-When you provide coordinates for all axes but one, you get a profile of values along the remaining axis, as shown above. If you provide fewer axis coordinates you get progressively higher-order results. To get a latitudinal transect, for instance, provide only a longitude coordinate:
-
-```{r transect}
-(trans29_74 <- t2m$profile(longitude = 29.74, .names = "lon_29_74"))
-```
-
-Note that there is only a single longitude coordinate left, at exactly the specified longitude.
-
-## Summarising data over time
-
-With the `summarise()` method you can apply a function over the data to generate summaries. You could, for instance, summarise daily data to monthly means. These methods use the specific calendar of the "time" axis. The return value is a new `CFVariable` object.
-
-```{r summarise}
-# Summarising hourly temperature data to calculate the daily maximum temperature
-t2m$summarise("tmax", max, "day")
-```
-
-A function may also return a vector of multiple values, in which case a list is returned with a new `CFVariable` object for each return value of the function. This allows you to calculate multiple results with a single call. You could write your own function to tailor the calculations to your needs. Rather than just calculating the daily maximum, you could get the daily maximum, minimum and diurnal range in one go:
-
-```{r summarise_multi}
-# Function to calculate multiple daily stats
-# It is good practice to include a `na.rm` argument in all your functions
-daily_stats <- function(x, na.rm = TRUE) {
-  # x is the vector of values for one day
-  minmax <- range(x, na.rm = na.rm)
-  diurnal <- minmax[2L] - minmax[1L]
-  c(minmax, diurnal)
-}
-
-# Call summarise() with your own function
-# The `name` argument should have as many names as the function returns results
-(stats <- t2m$summarise(c("tmin", "tmax", "diurnal_range"), daily_stats, "day"))
-```
-
-Note that you may have to update some attributes after calling `summarise()`. You can use the `set_attribute()` method on the `CFVariable` objects to do that.
-
-## Create new netCDF objects
-
 You can convert a suitable R object into a `CFVariable` instance quite easily. R objects that are supported include arrays, matrices and vectors of type logical, integer, numeric or logical.
 
-```{r create_basic}
-arr <- array(rnorm(120), dim = c(6, 5, 4))
-as_CF("my_first_CF_object", arr)
-```
-
-Usable but not very impressive. The axes have dull names without any meaning and the coordinates are just a sequence along the axis.
- 
 If the R object has `dimnames` set, these will be used to create more informed axes. More interestingly, if your array represents some spatial data you can give your `dimnames` appropriate names ("lat", "lon", "latitude", "longitude", case-insensitive) and the corresponding axis will be created (if the coordinate values in the `dimnames` are within the domain of the axis type). For "time" coordinates, these are automatically detected irrespective of the name.
 
 ```{r create_dimnames}
-# Note the use of named dimnames here - these will become the names of the axes
+# Note the use of named dimnames: these will become the names of the axes
+arr <- array(rnorm(120), dim = c(6, 5, 4))
 dimnames(arr) <- list(lat = c(45, 44, 43, 42, 41, 40), lon = c(0, 1, 2, 3, 4), 
                       time = c("2025-07-01", "2025-07-02", "2025-07-03", "2025-07-04"))
 
-(obj <- as_CF("a_better_CF_object", arr))
+(obj <- as_CF("a_new_CF_object", arr))
 
 # Axes are of a specific type and have basic attributes set
 obj$axes[["lat"]]
@@ -179,34 +57,7 @@ obj$axes[["time"]]
 
 You can further modify the resulting `CFVariable` by setting other properties, such as attributes or a coordinate reference system. Once the object is complete, you can export or save it.
 
-## Exporting and saving data
-
-A `CFVariable` object can be exported to a `data.table` or to a `terra::SpatRaster` (3D) or `terra::SpatRasterDataset` (4D) for further processing. Obviously, these packages need to be installed to utilise these methods.
-
-```{r export}
-# install.packages("data.table")
-library(data.table)
-head(dt <- ts$data.table())
-
-#install.packages("terra")
-suppressMessages(library(terra))
-(r <- stats[["diurnal_range"]]$terra())
-terra::plot(r)
-```
-
-A `stars` object can be created from a `CFVariable` or a `CFDataset` with multiple variables with the function `stars::st_as_stars()`.
-
-```{r stars}
-library(stars)
-(st <- st_as_stars(ts))
-```
-
-A `CFVariable` object can also be written back to a netCDF file. The object will have all its relevant attributes and properties written together with the actual data: axes, bounds, attributes, CRS. The netCDF file is of version "netcdf4" and will have the axes oriented in such a way that the file has maximum portability (specifically, data will be stored in row-major order with increasing Y values).
-
-```{r save, eval=FALSE}
-# Save a CFVariable instance to a netCDF file on disk
-stats[["diurnal_range"]]$save("~/path/file.nc")
-```
+More detailed operations and options are given on the [package web site](https://r-cf.github.io/ncdfCF/).
 
 ## Development plan