Improvements and fixes to NetCDFWriter

[tomchor]

This PR contains:

• A few improvements that will make Support ImmersedBoundaryGrid reconstruction using NetCDF #4848 easier.
• A couple of fixes for edge cases.
• Include dimension_type to NetCDFWriter which controls the precision type for dimensions, with a default value of Float64.
• Perhaps most importantly, it brings the treatment of arrays by the NetCDFWriter closer to what Oceananigans does. In particular singleton dimensions are kept if the location has a well-defined value (i.e. the value is not nothing, which NetCDF can't handle). This means will remove singleton dimensions only when:
  • the location on that dimension is nothing (e.g. a reduced field)
  • the topology in that dimension is Flat

To- do list:

• Should we enforce that spatial dimension is always Float64 to avoid alignment issues? We already do this with time for the same reason (see NetCDFOutputWriter sometimes outputs twice at approximately the same time step #3056).
• Special treatment for free surfaces
• Do not squeeze singleton dimensions in slices
• Add test for when the grid passed to NetCDFWriter is different from model.grid
• Add test to make sure that singleton dimensions aren't improperly squeezed/removed

[glwagner]

Suggested change

	indices = Any[:,:,:]
	indices = Any[:, : ,:]

what does this mean exactly?

[tomchor]

The comment is meant for the whole algorithm, not just that line. Basically it gets Fields that have Nothing in any of the locations and it gets rid of that dimension. (I'm calling these reduced_dims but there's probably a better name.)

So for example, if data comes from an array that has loc = Center, Center, Nothing, it'll reduce its parent data by indexing it as data[:, :, 1]; i.e. a 2D array in this case.

I'm testing this out (hence a draft PR) after some discussions with @ali-ramadhan. Basically we want to keep all NetCDF arrays 3D (to mimic Oceananigans' behavior) as much as possible, but in NetCDF there's no concept of location, so we can't easily pass LZ = Nothing to NetCDF. In these cases it probably makes sense to squeeze that dimension (or whatever the verb for this is) so that we know that dimension is reduced.

[glwagner]

Can you add attributes that indicate the location? We can add custom attributes to data even if they are not built in, right?

[tomchor]

I can (and in fact I probably will do it regardless), but that doesn't affect calculations with any of the software (that I'm aware of) which reads NetCDF. So broadcasting and other useful things wouldn't work. Without the singleton dimensions (i.e. with the change that I'm exploring here), those things will work out of the box for most software.

Also if we include the singleton reduced dimension, we still would need to pick a value for that length-1 coordinate. In Oceananigans, that value is nothing I think, but there's no analog to that for NetCDF, so we'd need to pick an arbitrary value like 1 or 0, which I'd like to avoid.

[tomchor]

Also to be clear, this PR would move the NetCDF behavior closer to Oceananigans' behavior. Currently the NetCDF writer also drops singleton dimensions whenever we slice a field, even though that sliced dimension (or length 1) has a proper value. This PR would change the behavior and keep that sliced length-1 dimension in the NetCDF, just like Oceananigans does.

[tomchor]

I think this PR is ready to be reviewed (the only failure is some keyring docs thing). I updated the top post to include a description of what was done. The most important thing is that slices with a singleton dimension now keep the singleton dimension, which is what Oceanigans does.

I'm only removing singleton dimensions in the data when there isn't a well-defined coordinate value associated with it. Namely when the dimension is Flat and when that direction is Reduced. In both cases the coordinate value in the singleton dimension is nothing, which NetCDF doesn't accept. To be clear, this was already done before this PR; I'm just keeping this behavior.

For the record, if I try to replace nothing with a placeholder value (like a missing flag which gets filled or a NaN) then broadcasting will not work properly (at least with all of the tools that I am aware of to read NetCDF). So I do think that's a necessary departure from the Oceananigans behavior to ensure a smooth user experience with the NetCDF files.

[tomchor]

@glwagner @ali-ramadhan one thing I'd like to discuss is whether we should make the spatial dimensions always Float64. I've faced alignment problems before when working with two separate NetCDF files from the same simulation: one in Float64 for calculations, and another in Float32 for plotting. We already enforce Float64 with the time dimension, so I think it'd be natural to do it with xyz as well.

[tomchor]

Btw, after this PR it'll be straightforward to merge #4848 and then it should be straightforward to tackle #3935 in one PR!

[glwagner]

@glwagner @ali-ramadhan one thing I'd like to discuss is whether we should make the spatial dimensions always Float64. I've faced alignment problems before when working with two separate NetCDF files from the same simulation: one in Float64 for calculations, and another in Float32 for plotting. We already enforce Float64 with the time dimension, so I think it'd be natural to do it with xyz as well.

happy for that. You could also make it user choice, with default Float64? That would be the Oceananigans Way.

[tomchor]

@glwagner @ali-ramadhan one thing I'd like to discuss is whether we should make the spatial dimensions always Float64. I've faced alignment problems before when working with two separate NetCDF files from the same simulation: one in Float64 for calculations, and another in Float32 for plotting. We already enforce Float64 with the time dimension, so I think it'd be natural to do it with xyz as well.

happy for that. You could also make it user choice, with default Float64? That would be the Oceananigans Way.

Sounds good! Just added that option with a Float64 default, plus a test :)