Stars: Data reading and writing in chunked NetCDF files

[adrfantini]

Just a heads up.
As you know NetCDF4 files are chunked: the data is stored in big data blocks which can be of any dimension. The chunks can be zipped. Different variables can have different chunking schemes and compression.

Say I have a XYT 100x100x50 dataset. If I want to access a single time slice, a chunking of 100x100x1 makes much more sense than a chunking of 10x10x50. The opposite of course is true if I want to read time series.

1. When writing NetCDF files (if this option will be provided), stars should provide an interface to choose the chunking. Zipping should be an option too.
2. When opening NetCDF files, stars should provide information about the chunking.
3. When reading data from NetCDF files, stars should use the most efficient access method given the chunk sizes, i.e. read as few chunks as possible and avoid reading them multiple times (yes, some programs do this, ignoring the chunking specs).

AFAIK this should all be doable through GDAL, but I might be wrong.

[edzer]

@adrfantini what makes you think this can be done through GDAL?

[adrfantini]

Compression is supported through the COMPRESS=DEFLATE option, and I assumed the chunking was available too (supporting compression and not chunking seems weird to me).

https://trac.osgeo.org/gdal/wiki/NetCDF reports that the new driver supports chunking via HDF5, but I am not sure of what kind of support that is.

[edzer]

The gdal netcdf driver document doesn't mention chunking. The options arg to st_stars can be used to specify opening options. Writing is not yet supported by stars.

I ran into a big 4D netcdf file with funny variable order that took ages to read using gdal, but seconds to read using raster (which uses the ncdf4 package for netcdf files).

[adrfantini]

Indeed it doesn't, but it must use some chunking options. Maybe they are not configurable, or is there any chance they are hidden?

Funnily enough I was just running some tests on the same file you used for the blog post and stars was 2x faster than raster. I guess it depends on the chunking, or maybe you were not reading all the layers with raster? I recalled that reading 4D files with raster wasn't possible?

[edzer]

If hidden, then hard-coded, see the source.

Good to hear stars is faster than raster, for once! Yes, the driver might be better optimized than the ncdf4 driver used (I assume) by raster. And indeed, raster will only read one 3D variable, stars potentially multiple 4D or higher-D ones.

[adrfantini]

Quite a lot of chunking is mentioned there, but with my limited knowledge of c++ and gdal I'm not able to identify where it sets the chunk dimensions.
Anyway, of course stars is set up to be so much more versatile than raster. There's no way I can translate all of my data analysis in raster, however stars looks much more promising. However performance will never be as good as NCO and CDO if GDAL pretty much ignores chunking when reading/writing weirdly shaped datasets. I'll do some tests maybe.

[mdsumner]

stars could align its access to files organised by tiles/chunks whether with GDAL or otherwise. I think this hinges a bit on possible moves by GDAL to multi-d support. But either way, using sf structures to map tiling/chunking and organize access based on that is a good stars topic. A related issue is raster being slow for tiled sources, because it doesn't leverage the tiles even though GDAL can provide them. Raster is always line-based, killing performance when tiles are present

[adrfantini]

With "tiled" you mean chunked in "horizontal" tiles?

I agree access patterns based on chunking is a must with NetCDF files (my point 3 above). It just destroys performance otherwise, I have seen programs taking 1s to read a file while others took 2 hours.

[mdsumner]

Yes I consider tile and chunk to be synonymous, there's no important difference AFAIK. It's not easy to design a general scheme, though, workflow specifics really matter

[mdsumner]

My only suggestion is to help collate examples of good and bad performance, with ncdf4, RNetCDF, raster and others, with real-world files and sources that are of interest. Having concrete examples is the best motivation - "this could be better", especially if the operation of existing tools is sub-optimal. Generally, using RNetCDF or ncdf4 is going to be fastest, but low-level, and specific - whereas raster will be general, but sometimes perform poorly, or not at all. These are ever-present trade-offs, so motivating case-studies really matter.

[adrfantini]

Yeah, that I can do, I'll try to find some time, maybe during the upcoming break.

[edzer]

Would be great!

[edzer]

@adrfantini : any news on this one?

[adrfantini]

Uh, I guess the break came and passed and I did not find some time for it. I'll really try to.
I am flooded with work for EGU, but maybe in that work I find some good usage examples.
I'll set this as permanent "to read" in my inbox, I swear.

[edzer]

Should this be kept open?

[adrfantini]

Oops, sorry, this passed by and I forgot. I do believe that this is now partially covered by stars_proxy objects both in reading and in writing.

As a sidenote, it's not clear to me what direction do we want to take with the read_stars/read_ncdf coexistence. read_ncdf is not able to create stars_proxy objects (iirc), instead it allows for manual reading of chunks. Do we want to reach feature parity between the two implementations?

As for writing, the chunk_size option only exists for proxy objects, or at least that's how it appears from the docs. This can be another issue.

[dblodgett-usgs]

The goal of the work I was doing in read_ncdf was to make this more possible. The difficulty will be that there is great diversity of NetCDF sources that will "work" with read_ncdf but not, for example, get assigned an accurate coordinate reference system. It is possible to do this with pure NetCDF, but it's not trivial for non-NetCDF-CF sources. That's just one example of many that would need to be handled to give all the needed functionality.

If there's a desire to go in the direction of read_ncdf having parity with read_stars, I'd be happy to help with details of the NetCDF implementation.

[adrfantini]

I think you only want to target NetCDF-CF sources. The files non compliant with the conventions can sometimes be read by read_stars or, well, not read at all using high-level packages.

[edzer]

Some observations:

• read_stars (GDAL) has huge problems with certain NetCDF files, taking essentially infinity, I suspect depending on axis ordering
• read_ncdf does not handle time properly (to POSIXct, PCICt, or even udunits), nor attribute units
• read_ncdf does everything in memory
• read_ncdf does not (?) deal with the coordinate reference systems properly
• handling units, time, and curvilinear grids with read_stars (GDAL) is clunky, handled in R, but mostly works
• write_stars does not write NetCDF's properly - it uses the GDAL model, everything becomes band, time gets lost

so until all this has been sorted out, or is considered fine, we still need both... Ideal would be an integration, a single function that kind of like always works.

[dblodgett-usgs]

Good observations. Some comments.

• Unfortunately, I don't think using GDAL for NetCDF is going to be workable for the breadth of functions desired in stars. I'd love to be proven wrong and there may be some hybrid approach that would primarily use GDAL but handle time independently.
• RNetCDF's time utilities seem quite good but my experience with them is limited.
• Minimal effort would be required to have read_ncdf work with proxy objects -- the key is getting axes of NetCDF sources mapped onto axes of a stars proxy. The underlying NetCDF machinery has the functionality to do the rest easily.
• I tried to contribute coordinate reference system handling but there were dependency issues. I know we can work around those and get CRS handling directly -- or we could possibly use a hybrid approach where GDAL is used to get the CRS information.
• Writing NetCDF is not very hard from a stars object if you are explicit and opinionated about how to write it out. I would enjoy contributing a simple NetCDF-CF serializer at some point if a stand alone writer were of interest.

[edzer]

I agree, and that would be great. I would be happy to sort out the time & units stuff, essentially did that the hard way (parsing METADATA tags) in the read_stars/GDAL path, and did the units package which forms the basis for netcdf (when properly used). Indeed, having GDAL sorting out the CRS would be a good option.

One thing we haven't done and stars can but GDAL can't do is rectilinear grids (netcdf may have dimension boundaries), and the information whether a cell refers to a point or an aggregate over an interval for each dimension.

[dblodgett-usgs]

Curious to get @mdsumner's input here. Some of what I contributed over in ncmeta is useful here. I wonder if working on this in a branch to deal with dependencies and some longer-term development might be the right way to go? I brought in CRS handling and some axis-order stuff that will be needed long term in: https://github.com/r-spatial/stars/pull/88 but it dropped out trying to get things released.

[mdsumner]

This is good discussion!

Can we take this discussion to #62