Access_* functions using THREDDS servers

[ErikKusch]

ClimHub contains functions called Access_KlimaiNorge2100 and Access_NORA3 which pull data from a thredds server. It should be investigated whether this one can be improved or made more robust potentially with ncdfcf.

[pvanlaake]

From looking at the code it is easy to see that multiple topics are included in a single function:

• Domain knowledge: where to find the data, layout of the THREDDS server, how file names are structured, what they contain (scenarios, models, variables, time periods).
• Catalog discovery: What data is already available locally.
• Data access: download the data from the THREDDS server and store it somewhere locally.

I would suggest that we break these things into separate functions. Comparing the KSS and NORA3 data access functions right now reveals that there is substantial overlap in the functions. That will only get worse as more domains are added. Breaking up the functions in these blocks would also make them better aligned with the 4 building blocks that are in the graphic on the landing page of the package.

Domains

There is an infinite number of domains (for all practical intents and purposes) that may be of interest to users of ClimHub. I would suggest that we make a catalog of them and include them with the package as an Rds file. We should see about all the details of file names and contents and try to automate creating URIs from the information that can plug into the next steps.

This should also include data sets other than netCDF on THREDDS servers. Topics such as topography, hydrology networks, soils, land cover, etc should all be included.

Extending the included domains should then be relatively straightforward and may build on community support.

Catalog

We need an arrangement that is user-led, non-intrusive, not limiting. I would suggest we treat this in the same way as the domains above: facilitate but not require it. As above, it should produce an URI, most likely a simple filename in this case, to be used in the next step.

Data access

When we are talking netCDF resources on THREDDS servers then the solution is there. The ncdfCF package uses a lazy approach: first it reads metadata which is used to subset data sets by variable and then over any of the dimensions of the data (spatial, temporal, others) before downloading just the data that the user needs.

For other data sets, like those mentioned above, we'd need a different solution.

I would not store the data locally immediately after downloading. Users are probably going to use the downloaded data in some data processing flow and may be interested in keeping only the result of that processing. I'd say we leave persisting data to a direct user command and in a separate function. That would also allow us to use other back-ends for persistence, such as Zarr (I am working on a package that should be ready for CRAN before the end of the year).

In short, I think we should retire these two functions and define the interface between the building blocks. If you agree, I can work on the data access part, which would also cover #7 and #6.

[ErikKusch]

Fully agree with more compartmentalisation and happy for you to take point on it!

We are already planning to source relevant data needs from the community via #20 but will need to build our own overview and functionality until then (@RichardDavy84 is compiling an overview for whats needed via #12).

Currently, @thisistaimur is planning to look into improvements to #6 so there may be overlap in your work.

[pvanlaake]

The below script is an example of what the ACCESS* functions could look like. It loops over all the files in the NORA3 data set within a defined date range, does a possible spatial subset, and then downloads the relevant data as CFVariable instances that are merged into a single CFVariable for each requested variable.

There is limited testing of argument validity: if the arguments are bad, the caller simply gets no data.

This is pretty much bare-bones code, many things can and should be improved. But it gives you a good idea of what the code could look like and what outputs are available for further processing.

The CFDataset can be written to file - so with multiple data variables - with ds$save().

You need the latest dev install from ncdfCF.

#' Extract NORA3 time series
#'
#' @param .vars Character vector with data variables to extract.
#' @param .dates Vector of a date range to extract, either as character strings,
#'   `Dates` or `POSIXct`.
#' @param .lead_time Optional. Integer value with the lead time of the forecast.
#'   6 (default) and 9 give the most complete data.
#' @param .verbose Logical that indicates if progress messages should be displayed.
#' @param ... Arguments of the form `axis = c(minimum, maximum)` to subset the
#'   spatial extent.
#' @returns A `CFDataset` instance with the extracted data variables.
#' @export
#' @examples
#' ds <- nora3(.vars = c("T2M", "Q2M"), .dates = c("2025-06-01", "2025-06-30"),
#'             longitude = 5:8, latitude = 55:57)
nora3 <- function(.vars, .dates, .lead_time = 6L, .verbose = FALSE, ...) {
  dates <- suppressWarnings(range(as.Date(.dates), na.rm = TRUE))
  if (is.infinite(dates[1L]) || is.na(dates[1L]))
    stop("Bad `.dates` argument")
  dates <- seq(from = dates[1L], to = dates[2L], by = 1L)
  dates_pack <- format(dates, "%Y%m%d")
  dates_slash <- format(dates, "%Y/%m/%d")

  .vars <- unique(.vars)
  dv <- vector("list", length(.vars))
  names(dv) <- .vars

  if (.verbose) cat("Downloading NORA3 data...")
  hr <- c("00", "06", "12", "18")
  base <- "https://thredds.met.no/thredds/dodsC/nora3/%s/%s/fc%s%2$s_%03d_sfx.nc"
  for (d in seq_along(dates)) {
    for (h in seq_along(hr)) {
      if (.verbose) cat("\n", dates_slash[d], hr[h], ":")
      ds <- open_ncdf(sprintf(base, dates_slash[d], hr[h], dates_pack[d], .lead_time))
      if (inherits(ds, "CFDataset")) {
        for (v in seq_along(.vars)) {
          if (.verbose) cat("", .vars[v])
          var <- ds[[.vars[v]]]
          if (inherits(var, "CFVariable")) {
            subs <- var$subset(...)
            subs$detach() # Make sure we grab the data
            if (is.null(dv[[v]])) dv[v] <- list(subs)
            else dv[[v]]$append(subs, "time")
          }
        }
      } else cat("(no file)")
    }
  }

  ds <- create_ncdf()
  ds$set_attribute("source", "NC_CHAR", paste("NORA3 (DOI:10.5194/wes-6-1501-2021) data provided by the The Norwegian Meteorological institute obtained on", format(Sys.time(), usetz = TRUE, digits = 0L)))
  ds$set_attribute("comment", "NC_CHAR", "Testing NORA3 download, subsetting, merging and saving")
  for (v in seq_along(dv))
    if (!is.null(dv[[v]]))
      ds$add_variable(dv[[v]])
  ds
}

[ErikKusch]

Just made some commits/pushes to the Access function branch.

A few heads-up and requests for some input:

1. Our solution for saving CFVariable objects that I have ontroduced in 4640235 is extremely rudimentary. Please investigate how this may be improved. CFVariable objects written to disk this way don't actually end up on disk with the additional attributes I am setting just above in the Access_NORA3() function - @pvanlaake, I am probably missing something basic, no? I imagine the rework of the file saving I implemented there may feed directly into NC Write Operations #6.
2. I have started to introduce a provenance attribute to the files we create to address Provenance String #25 (but this attribute is currently not saved properly)
3. I have reworked the metadata json file for NORA3 and the branch currently makes Discovery_* functionality target this development branch to test these changes to metadata. We will likely move away from json implementation in the future, but I figured some quality of life changes might be good right now.

[pvanlaake]

On point 1:

Commit R-CF/ncdfCF@5b18ea7 fixes the writing of attributes to new netCDF files, as well as writing multiple CFVariables to a single CFDataset, reading terra:SpatRaster objects and other improvements.

[ErikKusch]

With e50f14d, Access_* functions for NORA3 and Klima i Norge 2100 are now reworked to produce CFVariable products and can subset both for time and space (big thanks to @pvanlaake ) for the quick rework of some ncdfCF to facilitate this.

Only thing left here is to look into an alternative NC_writefor CFVariable objects so a case for #6

[pvanlaake]

I was looking at the Access_* functions already and I have several improvements to make. I'll work on the latest commit and will report back in a few days.

I am about to submit a new release of ncdfCF ro CRAN so will hold off until that is generally available on CRAN. Should be no more than a few days.

[pvanlaake]

ncdfCF is now out in version 0.8.0 and that supports all currently necessary features for ClimHub. It may be a few days before the Windows download is available so I have not changed the version specification in the DESCRIPTION file for now.

I have pushed a commit to the 18-* branch with a number of significant changes. In summary:

• There is now a clear separation of functions between Access*, Helper_AccessCF.R and NCRead.R. Access* functions have domain knowledge of the data sources so that is where file name and subset parameter composition takes place. Helper_AccessCF.R performs some common tasks for the Access* functions, such as looping over URLs and appending data from different URLs. NCRead.R does the low-level netCDF work, including subsetting with parameters prepared in the Access* functions.
• Saving to file is optional. In either case, the Access* functions return a CFDataset instance.
• I have commented out the lat/lon extent check. Like this it is also possible to subset on the planar coordinates. All subsetting takes place in one command, which is far more efficient than first lon/lat and then time after.
• Various other minor improvements.

The biggest drag, apart from downloading the netCDF data, is the QuickFacts download. These files should be in an Rds file and distributed with the package.

The ClimHub package builds with an error (in HelperAttrOverride.R which can probably be deleted) but that does not impact use of the AccessKlimaiNorge2100 and AccessNORA3 functions.

Have a look, let me know.

[ErikKusch]

Thanks for the contribution, @pvanlaake. Really nice efficiency improvements in there! Will look into the documentation a bit more and then merge the branches.

To add: The AttrOverride is obosolete with the move to ncdfCF as the main object driver of the package. I will remove it with #33

[ErikKusch]

Afraid I have run into some issues with your latest commit, @pvanlaake. This happened as I attempted to add some examples to all functions (something I would recommend we do going forward for all functions). The issues I have run into are as follows:

1. Data gaps in downloaded data (this looks like it only affects NORA3 data)
2. Failure to write data which seems to happen when I do not specify an extent

Any suggestions for how to resolve these would be great 😅 Happy to hop on a call for this!

Let me demonstrate them 1-by-1:

Data gaps

When downloading data like so:

library(ggplot2)
library(ncdfCF)
NORA3 <- Access_NORA3(variable = "TS", dateStart = "2001-08-01 00", dateStop = "2001-08-01 00", leadTimeHour = 3)
ts <- NORA3[["TS"]]
ggplot() +
    geom_ncdf(data = ts$subset(time = "2001-08-01T03:00:00")) +
    coord_equal() +
    scale_fill_viridis_c()

I get this plot:

There is a lot of data missing here. Any idea why?

For reference, the way we have been treating data so far (simplified):

library(terra)
URL <- "https://thredds.met.no/thredds/fileServer/nora3/2001/08/01/00/fc2001080100_003_sfx.nc"
download.file(
    url = URL,
    destfile = "fc2001080100_003_sfx.nc",
    method = "wget"
)
terra_obj <- rast("fc2001080100_003_sfx.nc")
plot(terra_obj[["T2M_grib1_vLevel1053=2"]])

No such gaps can be seen.

Failure to Write Data

Seemingly spontaneously, the Access_* functions can fail to write the data:

NORA3 <- Access_NORA3(variable = "TS", dateStart = "2001-08-01 00", dateStop = "2001-08-01 00", leadTimeHour = 3, fileName = "NORA3_NoExtent.nc")

can produce the error:

Error in RNetCDF::att.delete.nc(h, nm, names[a]) : 
  NetCDF: Write to read only

This still results in the file being created but seemingly not data being added to it:

r$> NC_Read("NORA3_NoExtent.nc")
<Dataset> NORA3_NoExtent 
Resource   : NORA3_NoExtent.nc 
Format     : netcdf4 
Collection : Generic netCDF data 
Conventions: CF-1.13 
Has groups : FALSE 

Attributes:
 name        type    length value
 Conventions NC_CHAR  7     CF-1.13
 history     NC_CHAR 67     Created with R package ncdfCF 0.8.0.9000 on 202...

At first, I figured this might be an issue of writing the file to a OneDrive-enabled folder, but the issue persists even when OneDrive synching is paused.

The same is true for the behaviour of Access_KlimaiNorge2100:

## fails with same error
KiN_rast <- Access_KlimaiNorge2100(
    variable = "mean_air_temperature",
    dateStart = "1995-01-01",
    dateStop = "1995-01-10",
    model = "noresm-r1i1p1f1-hclim",
    method = "eqm",
    scenario = "ssp370",
    fileName = file.path(getwd(), "inst/extdata", "KiN_rast.nc")
)

## executes just fine
KiN_rast <- Access_KlimaiNorge2100(
    variable = "mean_air_temperature",
    dateStart = "1995-01-01",
    dateStop = "1995-01-10",
    model = "noresm-r1i1p1f1-hclim",
    method = "eqm",
    scenario = "ssp370",
    extent = c(7, 12, 58, 62),
    fileName = file.path(getwd(), "inst/extdata", "KiN_rast2.nc")
)

[pvanlaake]

The writing failure gives me the shivers. I know what causes it (axes that are not subset on auxiliary lat-long grids still link to the original source, which is not writable and shouldn't be written to) but solving that is going to require some deep thinking and hard work...

The data gaps gives me the shivers about terra. I downloaded the exact same file from the fileServer instance and then I do this:

library(RNetCDF)
nc <- open.nc("~/Downloads/fc2001080100_003_sfx.nc")
# Get the area from NE coast of Ireland to SW corner of Norway
ts <- var.get.nc(nc, "TS", start = c(400, 300, 1, 1), count = c(200, 200, 1, 1), unpack = T)
# Detailed look of Ireland data: Lots of NAs!
> ts[1:10, 1:10]
        [,1]   [,2]   [,3]   [,4]   [,5]   [,6]   [,7]   [,8]   [,9]  [,10]
 [1,] 285.91 285.87 285.83 285.80     NA 282.86 283.46     NA 284.03 283.99
 [2,] 285.85 285.26 285.46 285.74     NA     NA 283.08 283.72 283.99     NA
 [3,] 285.79 285.61 285.25 285.69 285.25 283.00 282.78 283.22     NA     NA
 [4,] 284.45 284.94     NA 285.36 285.09     NA 282.44 282.98 283.51 283.88
 [5,] 283.51 283.72     NA 283.98     NA 282.99 282.60 283.50     NA     NA
 [6,] 283.68     NA 283.57 283.42     NA 282.98 283.21 284.55 284.08 284.02
 [7,] 284.87     NA 284.05 284.71 285.04     NA 283.90 284.50 283.32 283.17
 [8,] 285.56 285.42     NA 285.21 285.31     NA 284.91 284.30 283.44 283.41
 [9,] 285.51 285.49 285.46 285.39 284.80 285.14 285.32 284.43     NA     NA
[10,] 285.47 285.44 285.42 285.39 285.33 285.34 285.31 284.90 284.49 283.53
image(ts)

I have no idea where terra gets its data from but it is not in the file.

As I understand it, NORA3 is focused on open sea data so are land surfaces supposed to be fully covered?