Repository: GenericMappingTools/pygmt

Issue:
PyGMT is too strict when checking for irregular grids
**Description of the problem**

<!-- Please be as detailed as you can when describing an issue. The more information we have, the easier it will be for us to track this down. -->

Since GMT only supports regular grids, PyGMT has a check to make sure all point spacings are the same ([here](https://github.com/GenericMappingTools/pygmt/blob/8a6db70f1a9786007705a3888300028af4044de1/pygmt/clib/conversion.py#L97)):

```python
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):
    raise GMTInvalidInput(
        "Grid appears to have irregular spacing in the '{}' dimension.".format(
            dim
        )
    )
```

The problem is that the `np.allclose` function has variable precision arguments to determine whether it fails or passes. So the degree to which we check for uniform spacing depends on coordinate units. This is not ideal as this check can fail in cases where the grid is evenly spaced within an acceptable accuracy for coordinates but maybe not for `allclose`. 

I'm not sure there is a good solution here since:

1. This will vary with units so it's hard to generalize.
2. This part of the code is deep in the internals so getting the `rtol` and `atol` arguments from the user into this part of the code will be awkward and bug prone.
3. Should this really be an exception if GMT may be OK with the coordinates?

I would propose turning this into a warning instead and let GMT fail if it really thinks it can't handle the grid. What do others think? I'm not sure that would solve the problem or if GMT will be equally strict.

**Full code that generated the error**

```python
grid = xr.open_dataarray("topography-earth-10arcmin.nc")
fig = pygmt.Figure()
fig.grdimage(grid)
fig.show()
```

Using this grid file: [topography-earth-10arcmin.nc.zip](https://github.com/GenericMappingTools/pygmt/files/7065984/topography-earth-10arcmin.nc.zip) (GitHub required a zip archive for this)

For context, this grid is coming from the [ICGEM](http://icgem.gfz-potsdam.de/home) calculation service. I tracked down the error to them using the equivalent of `np.arange` to generate the grid coordinates instead of `np.linspace` (which is the correct way to do it). 

**Full error message**

<details>

```
---------------------------------------------------------------------------
GMTInvalidInput                           Traceback (most recent call last)
/tmp/ipykernel_207888/3701624484.py in <module>
      1 fig = pygmt.Figure()
----> 2 fig.grdimage(grid)
      3 fig.show()

~/bin/conda/envs/rockhound-data/lib/python3.9/site-packages/pygmt/helpers/decorators.py in new_module(*args, **kwargs)
    580                     )
    581                     warnings.warn(msg, category=SyntaxWarning, stacklevel=2)
--> 582             return module_func(*args, **kwargs)
    583 
    584         new_module.aliases = aliases

~/bin/conda/envs/rockhound-data/lib/python3.9/site-packages/pygmt/helpers/decorators.py in new_module(*args, **kwargs)
    723                         kwargs[arg] = separators[fmt].join(f"{item}" for item in value)
    724             # Execute the original function and return its output
--> 725             return module_func(*args, **kwargs)
    726 
    727         return new_module

~/bin/conda/envs/rockhound-data/lib/python3.9/site-packages/pygmt/src/grdimage.py in grdimage(self, grid, **kwargs)
    171                 kwargs["I"] = stack.enter_context(shading_context)
    172 
--> 173             fname = stack.enter_context(file_context)
    174             arg_str = " ".join([fname, build_arg_string(kwargs)])
    175             lib.call_module("grdimage", arg_str)

~/bin/conda/envs/rockhound-data/lib/python3.9/contextlib.py in enter_context(self, cm)
    427         _cm_type = type(cm)
    428         _exit = _cm_type.__exit__
--> 429         result = _cm_type.__enter__(cm)
    430         self._push_cm_exit(cm, _exit)
    431         return result

~/bin/conda/envs/rockhound-data/lib/python3.9/contextlib.py in __enter__(self)
    115         del self.args, self.kwds, self.func
    116         try:
--> 117             return next(self.gen)
    118         except StopIteration:
    119             raise RuntimeError("generator didn't yield") from None

~/bin/conda/envs/rockhound-data/lib/python3.9/site-packages/pygmt/clib/session.py in virtualfile_from_grid(self, grid)
   1344         # guarantees that the copy will be around until the virtual file is
   1345         # closed. The conversion is implicit in dataarray_to_matrix.
-> 1346         matrix, region, inc = dataarray_to_matrix(grid)
   1347 
   1348         family = "GMT_IS_GRID|GMT_VIA_MATRIX"

~/bin/conda/envs/rockhound-data/lib/python3.9/site-packages/pygmt/clib/conversion.py in dataarray_to_matrix(grid)
     96         coord_inc = coord_incs[0]
     97         if not np.allclose(coord_incs, coord_inc):
---> 98             raise GMTInvalidInput(
     99                 "Grid appears to have irregular spacing in the '{}' dimension.".format(
    100                     dim

GMTInvalidInput: Grid appears to have irregular spacing in the 'longitude' dimension.
```

</details>

**System information**

Please paste the output of `python -c "import pygmt; pygmt.show_versions()"`:

```
PyGMT information:
  version: v0.4.1
System information:
  python: 3.9.0 | packaged by conda-forge | (default, Nov 26 2020, 07:57:39)  [GCC 9.3.0]
  executable: /home/leo/bin/conda/envs/rockhound-data/bin/python
  machine: Linux-5.13.11-1-MANJARO-x86_64-with-glibc2.33
Dependency information:
  numpy: 1.21.2
  pandas: 1.3.2
  xarray: 0.19.0
  netCDF4: 1.5.7
  packaging: 21.0
  ghostscript: 9.54.0
  gmt: 6.2.0
GMT library information:
  binary dir: /home/leo/bin/conda/envs/rockhound-data/bin
  cores: 8
  grid layout: rows
  library path: /home/leo/bin/conda/envs/rockhound-data/lib/libgmt.so
  padding: 2
  plugin dir: /home/leo/bin/conda/envs/rockhound-data/lib/gmt/plugins
  share dir: /home/leo/bin/conda/envs/rockhound-data/share/gmt
  version: 6.2.0
```


Discussion:

Patch:
diff --git a/pygmt/clib/conversion.py b/pygmt/clib/conversion.py
index 57b2f92c..7c94b788 100644
--- a/pygmt/clib/conversion.py
+++ b/pygmt/clib/conversion.py
@@ -1,6 +1,8 @@
 """
 Functions to convert data types into ctypes friendly formats.
 """
+import warnings
+
 import numpy as np
 import pandas as pd
 from pygmt.exceptions import GMTInvalidInput
@@ -95,9 +97,14 @@ def dataarray_to_matrix(grid):
         coord_incs = coord[1:] - coord[0:-1]
         coord_inc = coord_incs[0]
         if not np.allclose(coord_incs, coord_inc):
-            raise GMTInvalidInput(
-                f"Grid appears to have irregular spacing in the '{dim}' dimension."
+            # 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."


Test patch:
diff --git a/pygmt/tests/test_clib.py b/pygmt/tests/test_clib.py
index e88c04a8..b240015b 100644
--- a/pygmt/tests/test_clib.py
+++ b/pygmt/tests/test_clib.py
@@ -805,16 +805,18 @@ def test_dataarray_to_matrix_dims_fails():
         dataarray_to_matrix(grid)
 
 
-def test_dataarray_to_matrix_inc_fails():
+def test_dataarray_to_matrix_irregular_inc_warning():
     """
-    Check that it fails for variable increments.
+    Check that it warns for variable increments, see also
+    https://github.com/GenericMappingTools/pygmt/issues/1468.
     """
     data = np.ones((4, 5), dtype="float64")
     x = np.linspace(0, 1, 5)
     y = np.logspace(2, 3, 4)
     grid = xr.DataArray(data, coords=[("y", y), ("x", x)])
-    with pytest.raises(GMTInvalidInput):
+    with pytest.warns(expected_warning=RuntimeWarning) as record:
         dataarray_to_matrix(grid)
+        assert len(record) == 1
 
 
 def test_dataarray_to_matrix_zero_inc_fails():