GenericMappingTools / pygmt

A Python interface for the Generic Mapping Tools.
https://www.pygmt.org
BSD 3-Clause "New" or "Revised" License
757 stars 219 forks source link

Failing test with the GMT 6.5.0 #2511

Open seisman opened 1 year ago

seisman commented 1 year ago
__________________ test_earth_relief_03s_default_registration __________________

    def test_earth_relief_03s_default_registration():
        """
        Test that the grid returned by default for the 3 arc-second resolution has
        a "gridline" registration.
        """
        data = load_earth_relief(resolution="03s", region=[-10, -9.8, 4.9, 5])
        assert data.shape == (121, 241)
        assert data.gmt.registration == 0
        npt.assert_allclose(data.coords["lat"].data.min(), 4.9)
        npt.assert_allclose(data.coords["lat"].data.max(), 5)
        npt.assert_allclose(data.coords["lon"].data.min(), -10)
        npt.assert_allclose(data.coords["lon"].data.max(), -9.8)
>       npt.assert_allclose(data.min(), -2069.996)

../pygmt/tests/test_datasets_earth_relief.py:263: 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 

args = (<function assert_allclose.<locals>.compare at 0x7fb55d6cccc0>, array(-2131.9707, dtype=float32), array(-2069.996))
kwds = {'equal_nan': True, 'err_msg': '', 'header': 'Not equal to tolerance rtol=1e-07, atol=0', 'verbose': True}

    @wraps(func)
    def inner(*args, **kwds):
        with self._recreate_cm():
>           return func(*args, **kwds)
E           AssertionError: 
E           Not equal to tolerance rtol=1e-07, atol=0
E           
E           Mismatched elements: 1 / 1 (100%)
E           Max absolute difference: 61.9[747](https://github.com/GenericMappingTools/pygmt/actions/runs/4765800212/jobs/8471996968?pr=2435#step:18:748)0312
E           Max relative difference: 0.02993953
E            x: array(-2131.9707, dtype=float32)
E            y: array(-2069.996)
seisman commented 1 year ago

GMT 6.4 and 6.5 produce different grids:

# Using GMT 6.4.0
$ gmt grdcut @earth_relief_03s_g -R-10/-9.8/4.9/5.0 -Grelief-6.4.nc
gmt grdinfo relief-6.4.nc 
relief-6.4.nc: Title: 
relief-6.4.nc: Command: gmt grdcut @earth_relief_03s_g/ -R-10/-9.8/4.9/5.0 -Grelief-6.4.nc
relief-6.4.nc: Remark: 
relief-6.4.nc: Gridline node registration used [Geographic grid]
relief-6.4.nc: Grid file format: nf = GMT netCDF format (32-bit float), CF-1.7
relief-6.4.nc: x_min: -10 x_max: -9.8 x_inc: 0.000833333333333 (3 sec) name: longitude n_columns: 241
relief-6.4.nc: y_min: 4.9 y_max: 5 y_inc: 0.000833333333333 (3 sec) name: latitude n_rows: 121
relief-6.4.nc: v_min: -2069.99609375 v_max: -924.080078125 name: z
relief-6.4.nc: scale_factor: 1 add_offset: 0
relief-6.4.nc: format: netCDF-4 chunk_size: 241,121 shuffle: on deflation_level: 3
relief-6.4.nc: Default CPT: geo

# Using 6.5.0_3122cc8_2023.04.21
$ gmt grdcut @earth_relief_03s_g -R-10/-9.8/4.9/5.0 -Grelief-6.5.nc
$ gmt grdinfo relief-6.5.nc                                        
relief-6.5.nc: Title: 
relief-6.5.nc: Command: gmt grdcut @earth_relief_03s_g/ -R-10/-9.8/4.9/5.0 -Grelief-6.5.nc
relief-6.5.nc: Remark: 
relief-6.5.nc: Gridline node registration used [Geographic grid]
relief-6.5.nc: Grid file format: nf = GMT netCDF format (32-bit float), CF-1.7
relief-6.5.nc: x_min: -10 x_max: -9.8 x_inc: 0.000833333333333 (3 sec) name: longitude n_columns: 241
relief-6.5.nc: y_min: 4.9 y_max: 5 y_inc: 0.000833333333333 (3 sec) name: latitude n_rows: 121
relief-6.5.nc: v_min: -2131.97070312 v_max: -924.080078125 name: z
relief-6.5.nc: scale_factor: 1 add_offset: 0
relief-6.5.nc: format: netCDF-4 chunk_size: 241,121 shuffle: on deflation_level: 3
relief-6.5.nc: Default CPT: geo
seisman commented 1 year ago

Here is the difference between the two grids:

gmt grdmath relief-6.5.nc relief-6.4.nc SUB = diff.nc
gmt grdimage diff.nc -Cturbo -Baf -png diff

diff

Obviously, the two grids diff at the lower-left corner.

Then I tried to compare the two grids in a larger region:

# GMT 6.4
$ gmt grdcut @earth_relief_03s_g -R-10.5/-9.8/4.6/5.0 -Grelief-6.4.nc
# GMT 6.5
$ gmt grdcut @earth_relief_03s_g -R-10.5/-9.8/4.6/5.0 -Grelief-6.5.nc

$ gmt grdmath relief-6.5.nc relief-6.4.nc SUB = diff.nc
$ gmt grdimage diff.nc -Cturbo -Baf -png diff

Here is the diff: diff The two grids diff at (-10, 4.6).

@PaulWessel Do you have any ideas why GMT 6.4/6.5 produce different grids and which one is correct?

PaulWessel commented 1 year ago

You see the same thing, @Esteban82 ?

Esteban82 commented 1 year ago

I can't test it for now because I have another problem. I get this message when I run anything with a remote data set.

*** buffer overflow detected ***: terminated

PaulWessel commented 1 year ago

Please post the command so we can test

Esteban82 commented 1 year ago

Please post the command so we can test

Now it works. I thinks that is something wrong in gmt_data_server.txt in the TEST server. We could see it later. I got the error with this command for example this. gmt grdcut @earth_relief_03s_g -R-10.5/-9.8/4.6/5.0 -Grelief-6.4.nc

We could see it later.

Esteban82 commented 1 year ago

You see the same thing, @Esteban82 ?

Yes, I also get a diff at 10, 4.6.

$ gmt grdinfo diff.nc 
diff.nc: Title: Produced by grdmath
diff.nc: Command: gmt grdmath relief-6.5.nc relief-6.4.0.nc SUB = diff.nc
diff.nc: Remark: 
diff.nc: Gridline node registration used [Geographic grid]
diff.nc: Grid file format: nf = GMT netCDF format (32-bit float), CF-1.7
diff.nc: x_min: -10.5 x_max: -9.8 x_inc: 0.000833333333333 (3 sec) name: longitude n_columns: 841
diff.nc: y_min: 4.6 y_max: 5 y_inc: 0.000833333333333 (3 sec) name: latitude n_rows: 481
diff.nc: v_min: -96.9716796875 v_max: 659.251953125 name: z
diff.nc: scale_factor: 1 add_offset: 0
diff.nc: format: netCDF-4 chunk_size: 141,161 shuffle: on deflation_level: 3
diff.nc: Default CPT: geo

diff_Larger_Region