Failing doctest on pygmt.clib.conversion.dataarray_to_matrix with NumPy 2.0

[weiji14]

Description of the problem

The following doctest started failing on the GMT Dev Tests on 28 Jul 2023, previously passing on 26 Jul 2023. See https://github.com/GenericMappingTools/pygmt/actions/runs/5686923585/job/15414542417#step:16:723.

Doesn't seem to be an issue with upstream GMT, at least I don't see anything obvious in these 7 commits GenericMappingTools/gmt@f497711...c94da14.

The region variable is created here -

pygmt/pygmt/clib/conversion.py

Lines 88 to 117 in 855f6fa

	# Extract region and inc from the grid
	region = []
	inc = []
	# Reverse the dims because it is rows, columns ordered. In geographic
	# grids, this would be North-South, East-West. GMT's region and inc are
	# East-West, North-South.
	for dim in grid.dims[::-1]:
	coord = grid.coords[dim].values
	coord_incs = coord[1:] - coord[0:-1]
	coord_inc = coord_incs[0]
	if not np.allclose(coord_incs, coord_inc):
	# calculate the increment if irregular spacing is found
	coord_inc = (coord[-1] - coord[0]) / (coord.size - 1)
	msg = (
	f"Grid may have irregular spacing in the '{dim}' dimension, "
	"but GMT only supports regular spacing. Calculated regular spacing "
	f"{coord_inc} is assumed in the '{dim}' dimension."
	)
	warnings.warn(msg, category=RuntimeWarning)
	if coord_inc == 0:
	raise GMTInvalidInput(
	f"Grid has a zero increment in the '{dim}' dimension."
	)
	region.extend(
	[
	coord.min() - coord_inc / 2 * grid.gmt.registration,
	coord.max() + coord_inc / 2 * grid.gmt.registration,
	]
	)
	inc.append(coord_inc)

Might be related to changes in the dev version of numpy/pandas/xarray?

Minimal Complete Verifiable Example

from pygmt.datasets import load_earth_relief
# Use the global Earth relief grid with 1 degree spacing
grid = load_earth_relief(resolution="01d", registration="pixel")
matrix, region, inc = dataarray_to_matrix(grid)
print(region)

Full error message

=================================== FAILURES ===================================
_____________ [doctest] pygmt.clib.conversion.dataarray_to_matrix ______________
043         If the grid has more than two dimensions or variable grid spacing.
044 
045     Examples
046     --------
047 
048     >>> from pygmt.datasets import load_earth_relief
049     >>> # Use the global Earth relief grid with 1 degree spacing
050     >>> grid = load_earth_relief(resolution="01d", registration="pixel")
051     >>> matrix, region, inc = dataarray_to_matrix(grid)
052     >>> print(region)
Expected:
    [-180.0, 180.0, -90.0, 90.0]
Got:
    [np.float64(-180.0), np.float64(180.0), np.float64(-90.0), np.float64(90.0)]

/home/runner/work/pygmt/pygmt/pygmt/clib/conversion.py:52: DocTestFailure

System information

Failing on:

PyGMT information:
  version: v0.9.1.dev110+gd0ca06a0
System information:
  python: 3.11.4 | packaged by conda-forge | (main, Jun 10 2023, 18:08:17) [GCC 12.2.0]
  executable: /home/runner/micromamba/envs/pygmt/bin/python
  machine: Linux-5.15.0-1041-azure-x86_64-with-glibc2.35
Dependency information:
  numpy: 2.0.0.dev0+646.g3dd9dba09
  pandas: 2.1.0.dev0+1303.g577bb7239c
  xarray: 2023.7.1.dev18+g52f5cf1f
  netCDF4: 1.6.4
  packaging: 23.1
  contextily: 1.3.0
  geopandas: 0.13.2
  IPython: 8.15.0.dev
  rioxarray: 0.14.1
  ghostscript: 9.54.0
GMT library information:
  binary version: 6.5.0_c94da14_2023.07.27
  cores: 2
  grid layout: rows
  image layout: 
  library path: /home/runner/work/pygmt/pygmt/gmt-install-dir/lib/libgmt.so
  padding: 2
  plugin dir: /home/runner/work/pygmt/pygmt/gmt-install-dir/lib/gmt/plugins
  share dir: /home/runner/work/pygmt/pygmt/gmt-install-dir/share
  version: 6.5.0

Passing version:

PyGMT information:
  version: v0.9.1.dev110+gd0ca06a0
System information:
  python: 3.11.4 | packaged by conda-forge | (main, Jun 10 2023, 18:08:17) [GCC 12.2.0]
  executable: /home/runner/micromamba/envs/pygmt/bin/python
  machine: Linux-5.15.0-1041-azure-x86_64-with-glibc2.35
Dependency information:
  numpy: 2.0.0.dev0+554.g0aaa5d397
  pandas: 2.1.0.dev0+1288.g1a4ac0ecb8
  xarray: 2023.7.1.dev16+gbb501ba5
  netCDF4: 1.6.4
  packaging: 23.1
  contextily: 1.3.0
  geopandas: 0.13.2
  IPython: 8.15.0.dev
  rioxarray: 0.14.1
  ghostscript: 9.54.0
GMT library information:
  binary version: 6.5.0_f497711_2023.07.25
  cores: 2
  grid layout: rows
  image layout: 
  library path: /home/runner/work/pygmt/pygmt/gmt-install-dir/lib/libgmt.so
  padding: 2
  plugin dir: /home/runner/work/pygmt/pygmt/gmt-install-dir/lib/gmt/plugins
  share dir: /home/runner/work/pygmt/pygmt/gmt-install-dir/share
  version: 6.5.0

Diff

Dependency information:
-  numpy: 2.0.0.dev0+554.g0aaa5d397
-  pandas: 2.1.0.dev0+1288.g1a4ac0ecb8
-  xarray: 2023.7.1.dev16+gbb501ba5
+  numpy: 2.0.0.dev0+646.g3dd9dba09
+  pandas: 2.1.0.dev0+1303.g577bb7239c
+  xarray: 2023.7.1.dev18+g52f5cf1f
GMT library information:
-  binary version: 6.5.0_f497711_2023.07.25
+  binary version: 6.5.0_c94da14_2023.07.27

[seisman]

It's caused by upstream changes in NumPy 2.0.

References:

• https://numpy.org/devdocs/release/2.0.0-notes.html#representation-of-numpy-scalars-changed
• https://numpy.org/neps/nep-0051-scalar-representation.html#nep51

[weiji14]

Nice job tracking this down! What do you think we should do here? Seems like this only affects the repr, and the other tests work fine. Don't think we can xfail a doctest though, and we need to keep compatibility with NumPy 1.22-1.26 (whose beta just came out https://github.com/numpy/numpy/releases/tag/v1.26.0b1) still.

Since dataarray_to_matrix is somewhat internal, we could maybe just use print([float(r) for r in region])? Or see what pytest-doctestplus is doing at scientific-python/pytest-doctestplus#210

[seisman]

Or see what pytest-doctestplus is doing at scientific-python/pytest-doctestplus#210

Let's wait for more days or weeks.

[bsipocz]

Honestly, I'm not sure when we have a proper way to handle this in doctestplus. In the meantime, you can do the workaround in your conftest:

https://github.com/astropy/astroquery/blob/main/astroquery/conftest.py#L17

[seisman]

Currently, the returned region is a list of np.float64 object. I think it makes more sense to return a list of float values. Then, we no longer have to worry about the changes in the numpy scalar representation. See #3226 for the fix.

[seisman]

We still have the issue in the doctests workflow (see https://github.com/GenericMappingTools/pygmt/actions/runs/9629554292/job/26559112184):

_____________________ [doctest] pygmt.src.grdclip.grdclip ______________________
101 Example
102 -------
103 >>> import pygmt
104 >>> # Load a grid of @earth_relief_30m data, with a longitude range of
105 >>> # 10° E to 30° E, and a latitude range of 15° N to 25° N
106 >>> grid = pygmt.datasets.load_earth_relief(
107 ...     resolution="30m", region=[10, 30, 15, 25]
108 ... )
109 >>> # Report the minimum and maximum data values
110 >>> [grid.data.min(), grid.data.max()]
Expected:
    [183.5, 1807.0]
Got:
    [np.float64(183.5), np.float64(1807.0)]