Unnecessary warnings when working with some netcdf files

[uleysky]

I get a warning when using grdinfo with some netcdf files: Guessing of registration in conflict between x and y, using gridline

Reading the code showed that the problem is in line 826 of the gmt_nc.c file:

if (fabs (fmod (dummy[0], dy)) > threshold) {	/* Most likely pixel registration since off by dy/2 */

This check is incorrect for two reasons, firstly, it does not work correctly for grids that are offset relative to zero. Secondly, due to the use of the fmod function, which may give incorrect results due to rounding errors (see https://godbolt.org/z/KbdK4Yaa9).

Also, I read lines 818 and 819 with great interest )

dummy[0] = xy[0], dummy[1] = xy[header->n_rows-1];
if ((fabs(dummy[1] - dummy[0]) / fabs(xy[header->n_rows-1] - xy[0]) - 1.0 > 0.5 / (header->n_rows - 1)) && header->registration == GMT_GRID_NODE_REG)

[joa-quim]

due to rounding errors (see https://godbolt.org/z/KbdK4Yaa9).

If I compile it with MSVC the result is correct.

[uleysky]

Different results with different compilers. The fmod function is more unpredictable than I thought.

[joa-quim]

Different results with different compilers. The fmod function is more unpredictable than I thought.

Or a GCC bug?

Also, I read lines 818 and 819 with great interest )

dummy[0] = xy[0], dummy[1] = xy[header->n_rows-1]; if ((fabs(dummy[1] - dummy[0]) / fabs(xy[header->n_rows-1] - xy[0]) - 1.0 > 0.5 / (header->n_rows - 1)) && header->registration == GMT_GRID_NODE_REG)

OK, it does nothing that line but is not the source of the problem either.

it does not work correctly for grids that are offset relative to zero.

Do you mind to elaborate a bit more (I don't understand it)?

And a grid to grid to reproduce the problem is ... you know.

[uleysky]

Or a GCC bug?

This is unlikely to be a bug of any compiler. More like undefined behavior.

OK, it does nothing that line but is not the source of the problem either.

I was too lazy to fill out another bug report )

it does not work correctly for grids that are offset relative to zero. Do you mind to elaborate a bit more (I don't understand it)?

Something like 0.5, 1.5, 2.5, etc.

And a grid to grid to reproduce the problem is ... you know.

ptemp-sal-pdens-u-v_HYCOM_2022-09-17-09_8.nc.tar.gz Here is example. Just call grdinfo on this nc file.

[joa-quim]

Thanks, I get the same warning with a MSVC built GMT.

[PaulWessel]

no time to look, how about @joa-quim ?

[joa-quim]

No time for deep debug either but maybe something is related to the files coordinates. Here dy should be 0.04 and not

[PaulWessel]

Not fmod. GMT is a bit fickle when it comes to xmin or ymin not being a multiple of xinc and yinc. Here, the x min is 1724.96051879 times xinc and that interferes with the dumb check. I will see how to avoid this but big dinner to prep for today so...

[PaulWessel]

But also note that the y-vector we read from the file is junk:

69, 69.040000915527344, 69.080001831054688, ..., 70

So whatever tool write this file goofed and said yinc = 0.040000915527344 which makes no sense since the y-rage is 1 and 1 / 0.040000915527344 = 24.9994278085, a tad bit short of 25, no?

[uleysky]

HYCOM data obtained via OpenDAP. ncdump -v lat https://tds.hycom.org/thredds/dodsC/GLBy0.08/expt_93.0

[PaulWessel]

So ncdump -h tells me that the x and y vectors are floats:

	float latitude(latitude) ;
		latitude:standard_name = "latitude" ;
		latitude:long_name = "Latitude" ;

and dumping just the vector to stdout gives

latitude = 69, 69.04, 69.08, 69.12, 69.16, 69.2, 69.24, 69.28, 69.32, 69.36, 
   69.4, 69.44, 69.48, 69.52, 69.56, 69.6, 69.64, 69.68, 69.72, 69.76, 69.8, 
   69.84, 69.88, 69.92, 69.96, 70 ;

which does step 0.04, possibly due to formatting. However, in GMT we read that vector in via

if ((has_vector = !nc_get_var_double (ncid, ids[HH->xy_dim[1]], xy)) && header->n_rows > 1)

where xy is a double precision array. netCDF is supposed to handle the conversion between float and double in such calls. However, we get 69, 69.0400009, 69.0800018,..., 69.9599991, 70

For fun I tried nc_get_var_float with a float xy instead and it returns exactly the same rotten steps. So the float-> double conversion works as advertised. I can only draw two conclusions:

1. The netcdf file was written with inaccurate vector node information.
2. It is not unreasonable for us to give some warnings on such files since they are bad.

Maybe @remkos has some information about these things.

Hence, nothing to really fix here, but you could send a note to the HYCOM folks and ask why.

[uleysky]

1. I agree, the grid in the source file is not perfect.
2. Sounds reasonable. But the warning should be about a bad step, not about incorrect registration.

[joa-quim]

We will improve that test, but not urgent.