Enable writing of 3-D netCDF data cubes from greenspline and grdinterpolate

[PaulWessel]

Description of proposed changes

This PR completes the gmt_write_nc_cube coding to enable the writing of 3-D layers (with layers an unlimited dimension). It is WIP because I want to complete what GMT_Read_Data with CONTAINER_ONLY should do for cubes, and update gspline_4.sh to no longer write lots of 2-D grids.

Related question whose solution can be separate PR is:

• Do we add -Q for "cube" options to grd2xyz and xyz2grd so they can handle 4-D data or do we clone and write cube2xyzw/xyzw2cube? My vote would be for -Q since adding two new modules for this seems overkill. But we do want the ability to dump and load cubes.
• What about grdinfo for 3-D cubes? We currently have no way to determine the equivalent info for cubes that we create. I see -Q is available there too.

[PaulWessel]

Alternatively, no -Q is needed for grdinfo and grd2xyz since we could presumably check if a given netCDF file is a grid or cube, and then take action. As for xyz2grd one has to specify -Rw/e/s/n -Idx/dy, but if a cube this needs to be -Rw/e/s/n/z0/z1 -Idx/dxy/dz so perhaps there is no need for -Q there either. What do you think, @joa-quim ?

[joa-quim]

grd2xyz and xyz2grd so they can handle 4-D data

What is the interest of this? Why care to convert a 4-D to xyzt and back?
A more useful tool would a stacker, which would read a pack of grids and stack them into a 3-D nc file.

[PaulWessel]

grdinterpolate grid1 grid2 ... -Zlevels -Gcube.nc will do the stacking (in next commit).
We cared enough to convert grd <--> xyz so I guess it is not that different to do cube <--> xyzw ?

[joa-quim]

We cared enough to convert grd <--> xyz

Yes, > 20 years ago when people still wrote grids manually. I fail to see why spending time coding a 4D analog

[PaulWessel]

OK, we can wait to see if there is any demand. What about grdinfo for cubes?

[joa-quim]

Please break line. Lines like these are an horror to debug.

[PaulWessel]

OK, done in my workspace - not ready to commit yet.

[joa-quim]

grdinfo for cubes seems useful.

[PaulWessel]

Pinging @remkos to have a look if he has time. In addition to the email I sent, also see the issue with lack of assigning variable names and units a few posts above.

[remkos]

I was looking at your proposal for units. As far as I remember (but I may have done that in other programs), if one adds "[m]" to the x-name then "m" is used as the unit. This should also hold for y-name, z-name and any other specified units, of course.

But yes, this would be rather "hidden" way of doing things and probably a more direct way of specifying the units, as you suggested above would be better. Or maybe support both.

Indeed, there was no way to specify the variable name. It was always x and y, unless the projection was geographical and lon and lat where used. This is not required by COARDS, but I have seen applications react to it appropriately, even without the units like degrees_east, degrees_north or the appropriate standard_name. Still, what is required by COARDS is that the dimensions and the coordinate axes are called the same. So the dimensions and the coordinate arrays are called lon and lat, consistently.

Now checking the code:

• Yes, gmtnc_put_units does split off the [units] part from the name_units variable which can be either x_units or y_units. So if in the present form you would specify "x-unit" as "temperature [K]" it would write long_name = "temperature" and units = "K".
• There are also some provisions for time axes. But it does not look like that I made any provisions to put the time references in correctly, e.g. "seconds since 2000-01-01 00:00:00".

[PaulWessel]

It is great that name [unit] is parsed like that - I guess I thought it was a special thing for [degrees_longitude] etc. and not general. So we can improve the docs on that for sure.

[PaulWessel]

Since (re-)discovering the ability to specify the name of a dimension or data and add the unit in brackets, the proposal above has been simplified, and implemented. We now have, via the -Dinfo setting, these new things:

1. +ddname sets the name and unit of the data value (formerly just +zname for grids - this is still backwards compatible) [z]
2. +zname sets the 3rd dimension name and unit for cubes [z or time]
3. +vvarname names the netCDF data variable [z or cube]

With that I think we are good. I updated gspline_4.sh to use -D (which I added to greenspline) via the command

gmt greenspline -R12/32/0/6/5/10 -I0.25 -G3D_cube.grdTable_5_23.txt -Sr0.85 -Z5 -D+x"x-distance [km]"+y"y-distance [km]"+z"z-distance [km]"+d"Uranium Oxide [%]"+vuoxide"

and now ncdump -h gives

variables:
	double x(x) ;
		x:long_name = "x-distance" ;
		x:units = "km" ;
		x:actual_range = 12., 32. ;
	double y(y) ;
		y:long_name = "y-distance" ;
		y:units = "km" ;
		y:actual_range = 0., 6. ;
	double z(z) ;
		z:long_name = "z-distance" ;
		z:units = "km" ;
		z:actual_range = 5., 10. ;
	float uoxide(z, y, x) ;
		uoxide:long_name = "Uranium Oxide" ;
		uoxide:units = "%" ;
		uoxide:_FillValue = NaNf ;
		uoxide:actual_range = 0., 23.2000007629395 ;

So, the -Dinfo was added to both greenspline and grdinterpolate. Since greenspline already used -D I moved that option to -Z but we handle this with backwards compatibility. All tests work fine - I think this PR is now good to be approved and merged.

[joa-quim]

Tests pass for me too.