Stars: If netcdf files with identical timestamps are read in, only the first one is actually read

[adrfantini]

I can attach files if needed, but it's trivial to reproduce by taking one file and copying it with a different name.
This holds true for files with any number of timesteps.

read_stars(c("file.nc", "file_copy.nc"), sub="pr")
pr, 
pr, 
pr, 
stars object with 3 dimensions and 1 attribute
attribute(s):
 NETCDF:"file.nc":pr [kg/m^2/s]
 Min.   :0.000e+00                          
 1st Qu.:7.000e-08                          
 Median :4.038e-06                          
 Mean   :1.205e-05                          
 3rd Qu.:1.793e-05                          
 Max.   :4.492e-03                          
dimension(s):
     from  to              offset  delta                       refsys point
x       1 527            -3162000  12000 +proj=lcc +lat_1=30 +lat_...    NA
y       1 527             3162000 -12000 +proj=lcc +lat_1=30 +lat_...    NA
time    1   1 1970-03-31 01:30:00     NA                      POSIXct    NA
     values
x      NULL
y      NULL
time   NULL

[edzer]

I think this is a feature. Why would you want to have layers with identical time stamps in a single object?

[adrfantini]

I guess that's a matter of opinion, but a warning should at least be issued.

One example where this can happen to me using netCDF files, CDO and NCO: take 100 different simulations covering the same period which I want to compare with something. I take the time average of each of the simulations using NCO or CDO : this creates files with a single timestep (the average time). All the files will have the same time.
This might be a bad design decision on part of CDO and NCO, but there are reasons why this might be useful (e.g. being able to set standard time_bnds attributes and values).

This means that before being able to import them into stars as a single object I have to reopen them, remove the time, and add a stacking or ensemble dimension with a different value for each file.

Compare this with raster, where I can simply stack any number of files I want:

obs <- raster(observations)
s <- stack(filenames)
bias <- s - obs
#... do data analysis on the bias
#... or transform now to stars: st_as_stars(bias) will have a "band" dimension

This is a limitation in stars by which dimensions must always be tied to some (usually physical) quantity present in the files themselves beforehand. While I do agree that this is the most correct way of doing things, it is much more tedious, which is annoying especially for taking a quick and dirty look at the data (which is often most of a scientist's work, isn't it).

I understand that stars is this way by design; I'm not saying it shouldn't be so by default. However I think that either you should add a big warning when something like this happens, and/or that you should provide an option to layer multiple input files on an additional "virtual" variable.

[edzer]

There's not so much design yet behind it, we just have to make it work. The point is, all layers should have equal time, but should vary along a different dimension (like ensemble member, or simulation).

Right now I think there's only checking for regular dimensions (where we have offset and delta), not for rectilinear ones, where the dimension values are in the "values" field of the dimension. I don't think there is checking going on that these values should be unique, so the array model of having a single record for each unique combination of dimensions is not held: dimension index yes, but dimension values are not uniqued, right now.

[adrfantini]

The more I think about this the more I see how common of a use-case this is: opening the outputs from several different simulations and compare them on the fly. This is very easy to do with raster, which stacks the file regardless of the layer dimension.
Maybe the most convenient approach in stars would be an option to forcefully add a faux dimension in which to record from which file the data was taken. (it's not very elegant, is it)

[mdsumner]

Isn't it a new attribute? (what am I missing?)

[adrfantini]

I guess you could see it as a dimension but also as an attribute, yes. However I do not think you can perform actions (e.g. subsetting or using dplyr verbs) over attributes. My "solution" feels a bit hackish to me, i'm sure something more appropriate is possible.

[edzer]

@mdsumner if it is not clear (as with nc files of consecutive times), the set of arrays can be bound as a collection of attributes, or ordered as subarrays along a new dimension in a single new target array. I believe there's no way ahead of time to say what makes more sense. You can now specify this with the along arg to read_stars. along=NA will create a set of attributes, along="ensemble" will create a new dim called ensemble. The latter is easier if you want to reduce them using st_apply (which you typically want). Also, split and merge (or c) methods let you move stuff back and forth between attributes and dimension.

[adrfantini]

@edzer awesome! Tested, works. Feel free to close if you think this is definitive.

[edzer]

Thanks!

[mdsumner]

Thanks from me too, to both of you. A few interesting aspects here