Tif won't convert to grd without complaining about xres/yres

[cbworden-usgs]

I download a file from ScienceBase in .tif format, and grdinfo is happy with it, but as soon as I do anything with the file that turns it into a .grd (or .nc), grdinfo (and other GMT processes) throw an error: "grdinfo [ERROR]: (x_max-x_min) must equal (NX + eps) * x_inc), where NX is an integer and |eps| <= 0.0001. [etc.]". The file comes from this page on SB: https://www.sciencebase.gov/catalog/item/5a5fa029e4b06e28e9bfc43a

To get the .tif file do:

curl https://www.sciencebase.gov/catalog/file/get/5a5fa029e4b06e28e9bfc43a\?f=__disk__2b%2F4b%2Fac%2F2b4bacfcbe272e38b35637f7b9633401e1c2f191 > ca_map.tif

You may need to remove the "\" before the "?" depending on your shell.

Now:

gmt grdinfo ca_map.tif

is happy, and the numbers look good. But if we do:

gmt grdconvert ca_map.tif -Gca_map.grd

Now:

gmt grdinfo ca_map.grd

Gives:

grdinfo [ERROR]: (x_max-x_min) must equal (NX + eps) * x_inc), where NX is an integer and |eps| <= 0.0001.
grdinfo [ERROR]: (y_max-y_min) must equal (NY + eps) * y_inc), where NY is an integer and |eps| <= 0.0001.
grdinfo (gmtlib_read_grd_info): Use grdedit -A on your grid file to make region and increments compatible [ca_map.nc]
[...]

Despite the numbers still looking good. Following the instructions in the error message:

gmt grdedit -A ca_map.grd -Gca_map_fixed.grd

Results in the same errors when running

gmt grdinfo ca_map_fixed.grd

curl https://www.sciencebase.gov/catalog/file/get/5a5fa029e4b06e28e9bfc43a\?f=__disk__2b%2F4b%2Fac%2F2b4bacfcbe272e38b35637f7b9633401e1c2f191 > ca_map.tif

gmt grdconvert ca_map.tif -Gca_map.grd

gmt grdedit -A ca_map.grd -Gca_map_fixed.grd

gmt grdinfo ca_map_fixed.grd

curl https://www.sciencebase.gov/catalog/file/get/5a5fa029e4b06e28e9bfc43a\?f=__disk__2b%2F4b%2Fac%2F2b4bacfcbe272e38b35637f7b9633401e1c2f191 > ca_map.tif

Full error message

grdinfo [ERROR]: (x_max-x_min) must equal (NX + eps) * x_inc), where NX is an integer and |eps| <= 0.0001.
grdinfo [ERROR]: (y_max-y_min) must equal (NY + eps) * y_inc), where NY is an integer and |eps| <= 0.0001.

Actual outcome

bash-3.2$ curl https://www.sciencebase.gov/catalog/file/get/5a5fa029e4b06e28e9bfc43a\?f=__disk__2b%2F4b%2Fac%2F2b4bacfcbe272e38b35637f7b9633401e1c2f191 > ca_map.tif
  % Total    % Received % Xferd  Average Speed   Time    Time     Time  Current
                                 Dload  Upload   Total   Spent    Left  Speed
100 84.1M  100 84.1M    0     0   775k      0  0:01:51  0:01:51 --:--:--  822k
bash-3.2$ 
bash-3.2$ gmt grdconvert ca_map.tif -Gca_map.grd
bash-3.2$ gmt grdinfo ca_map.grd
grdinfo [ERROR]: (x_max-x_min) must equal (NX + eps) * x_inc), where NX is an integer and |eps| <= 0.0001.
grdinfo [ERROR]: (y_max-y_min) must equal (NY + eps) * y_inc), where NY is an integer and |eps| <= 0.0001.
grdinfo (gmtlib_read_grd_info): Use grdedit -A on your grid file to make region and increments compatible [ca_map.grd]
[...]
bash-3.2$ gmt grdedit -A ca_map.grd
bash-3.2$ gmt grdinfo ca_map.grd
grdinfo [ERROR]: (x_max-x_min) must equal (NX + eps) * x_inc), where NX is an integer and |eps| <= 0.0001.
grdinfo [ERROR]: (y_max-y_min) must equal (NY + eps) * y_inc), where NY is an integer and |eps| <= 0.0001.
[...]

Expected outcome

The grid spacing seems fine for the bounds. I don't think an error is warranted here, but in any event, I don't know how to fix it.

System information

• Operating system: macOS Sonoma 14.7.2
• GMT version (gmt --version): 6.5.0

[welcome]

👋 Thanks for opening your first issue here! Please make sure you filled out the template with as much detail as possible. We appreciate that you took the time to contribute!

Please make sure you read our Contributing Guide and abide by our Code of Conduct.

[joa-quim]

Hi,

The real issue here is that grdedit -A is not doing its job as advertised. The message itself is correct but IMO it should not show up as error but instead as a WARNING, and again IMO, it should only be shown under -V.

The x_inc|y_inc are not correct, but just by that tinny bit that triggers the message.

3 arc secs = 3*1/3600 = 0.0008333333333333334

not 0.000833333300008

[cbworden-usgs]

Hi, thanks for your comment.

I think it is an issue with x_inc/y_inc being treated as both single-precision and double-precision floating points and not being reset correctly by grdconvert or grdedit, but the code logic is too dense for me to follow. From using -Vd, it seems that grdconvert tries to reset x_inc/y_inc, but it doesn't seem to stick. And then grdedit -A tries again, but it just goes back to the original single-precision settings for some reason.

I may be barking up the wrong tree, but that's what seems to be going on.

[joa-quim]

Nope, all computations are carried in double precision. It's just the format printing that shows less decimals. I increased them and now -Vd says

grdedit [INFORMATION]: Reset grid-spacing in file ca_map.grd to 0.00083333330000801/0.00083333330000868

which equal to original (but it was recomputed from the -R limits). The issue is that the -A doesn't even try to fix the (NX + eps) part.

[joa-quim]

No, wrong conclusion above. -A computes the right inc but somehow when we run grdinfo the increments were slightly changed to 0.00083333333333333339 (the 3 arc sec) and that's why we have that warning that prints like an error
grdinfo [ERROR]: (x_max-x_min) must equal (NX + eps) * x_inc) ...

[joa-quim]

So here is what is happening (in grdinfo).

gmtgrdio_doctor_geo_increments(...)
/* Check for sloppy arc min/sec increments due to divisions by 60 or 3600 */

Changes the increments ...

grdinfo ca_map.grd -V
...
grdinfo [INFORMATION]: Round-off patrol changed geographic grid increment for longitude from 0.000833333300008011826 to 0.000833333333333333387

but not the limits so when latter gmt_minmaxinc_verify() compares the original limits with the modified increments we fall in a case where the message

grdinfo [ERROR]: (x_max-x_min) must equal (NX + eps) * x_inc), where NX is an integer and |eps| <= 0.0001.

is triggered.

Not sure what to do here other than move that annoying warning up to at least the -V level. @seisman ?

EDIT: But basically GMT is right the limits are slightly off such that the increments that were clearly intended to 3 arc sec, are in fact not. That grid also does not announce its nodata-value so ocean and non-California are filled with zeros instead of NaNs.

[seisman]

Here is my understanding of the issue.

The origianl GeoTIFF file is stored in float32. That's why the grid increment of 3 arc second is 0.000833333300008. In GMT, computations are done in double precision, so GMT think 0.000833333300008 is inaccurate and changes it to 0.000833333333333333387.

gdalinfo ca_map.tif
Driver: GTiff/GeoTIFF
Files: ca_map.tif
Size is 12483, 11523
Coordinate System is:
GEOGCRS["WGS 84",
    ENSEMBLE["World Geodetic System 1984 ensemble",
        MEMBER["World Geodetic System 1984 (Transit)"],
        MEMBER["World Geodetic System 1984 (G730)"],
        MEMBER["World Geodetic System 1984 (G873)"],
        MEMBER["World Geodetic System 1984 (G1150)"],
        MEMBER["World Geodetic System 1984 (G1674)"],
        MEMBER["World Geodetic System 1984 (G1762)"],
        MEMBER["World Geodetic System 1984 (G2139)"],
        MEMBER["World Geodetic System 1984 (G2296)"],
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]],
        ENSEMBLEACCURACY[2.0]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["Horizontal component of 3D system."],
        AREA["World."],
        BBOX[-90,-180,90,180]],
    ID["EPSG",4326]]
Data axis to CRS axis mapping: 2,1
Origin = (-124.501041666999996,42.101041666999997)
Pixel Size = (0.000833333300008,-0.000833333300009)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-124.5010417,  42.1010417) (124d30' 3.75"W, 42d 6' 3.75"N)
Lower Left  (-124.5010417,  32.4985421) (124d30' 3.75"W, 32d29'54.75"N)
Upper Right (-114.0985421,  42.1010417) (114d 5'54.75"W, 42d 6' 3.75"N)
Lower Right (-114.0985421,  32.4985421) (114d 5'54.75"W, 32d29'54.75"N)
Center      (-119.2997919,  37.2997919) (119d17'59.25"W, 37d17'59.25"N)
Band 1 Block=12483x1 Type=Float32, ColorInterp=Gray
  Min=-1.000 Max=1634.973 
  Minimum=-1.000, Maximum=1634.973, Mean=nan, StdDev=nan
  NoData Value=-1
  Metadata:
    STATISTICS_MAXIMUM=1634.9729003906
    STATISTICS_MEAN=nan
    STATISTICS_MINIMUM=-1
    STATISTICS_STDDEV=nan

Personally, I feel nothing is wrong here. Instead, you may want to use the following command to edit the grid. The new -R region is calculated by assuming x_min/y_min is correct, and x_max=x_min + NX * 0.000833333333333333387, y_max=y_min + NY * 0.000833333333333333387.

gmt grdedit -A -R-124.501041667/-114.09854166699999/32.498542051/42.101042051 ca_map.grd -Gca_map_fixed.gr

[joa-quim]

Yes, but that would require to do grdedit -A as well because by changing the limits the condition (x_max-x_min) must equal (NX + eps) * x_inc) will no longer hold. I tried to address this in #8660

[cbworden-usgs]

Thanks for the workaround. It seems a bit fragile, like if the file's author changes anything about the file by just a smidge it won't work anymore and I'll have to revise my script, but this will work for now. So thanks.

I do think this qualifies as a bug because the behavior is unexpected: 1) grdconvert takes a .tif file that grdinfo thinks is fine and turns it into a .grd file that is no longer compliant, then 2) the recommended solution, grdedit -A, does not fix it. So maybe it is two bugs. Or maybe joa-quim is right and the ERROR message should be downgraded to a WARNING since it doesn't prevent anything from actually working.

[joa-quim]

Bruce, the whole thing was fixed in #8660

The fix suggested by @seisman followed by grdedit -A is not fragile, but annoying I agree. The alternative is to ignore the warnings (yes, they are warnings though the ERROR tag)

[joa-quim]

Closing as fixed in #8660