Autodetect grid type as Cartesian or Geographic coordinate system

Keep the default as Cartesian, but allow for either Cartesian (c) or Geographic (g) coordinate system grid types. Test case based on #375 (comment).

pygmt/base_plotting.py

@@ -228,7 +228,7 @@ def colorbar(self, **kwargs):
         W="pen",
     )
     @kwargs_to_strings(R="sequence", L="sequence", A="sequence_plus")
-    def grdcontour(self, grid, registration=None, **kwargs):
+    def grdcontour(self, grid, coord_sys=None, registration=None, **kwargs):
         """
         Convert grids or images to contours and plot them on maps
 
@@ -275,6 +275,7 @@ def grdcontour(self, grid, registration=None, **kwargs):
         {G}
         {U}
         {W}
+        {coord_sys}
         {r}
         """
         kwargs = self._preprocess(**kwargs)
@@ -283,7 +284,7 @@ def grdcontour(self, grid, registration=None, **kwargs):
             if kind == "file":
                 file_context = dummy_context(grid)
             elif kind == "grid":
-                file_context = lib.virtualfile_from_grid(grid, registration)
+                file_context = lib.virtualfile_from_grid(grid, coord_sys, registration)
             else:
                 raise GMTInvalidInput("Unrecognized data type: {}".format(type(grid)))
             with file_context as fname:
@@ -293,7 +294,7 @@ def grdcontour(self, grid, registration=None, **kwargs):
     @fmt_docstring
     @use_alias(R="region", J="projection", W="pen", B="frame", I="shading", C="cmap")
     @kwargs_to_strings(R="sequence")
-    def grdimage(self, grid, registration=None, **kwargs):
+    def grdimage(self, grid, coord_sys=None, registration=None, **kwargs):
         """
         Project grids or images and plot them on maps.
 
@@ -307,6 +308,7 @@ def grdimage(self, grid, registration=None, **kwargs):
         ----------
         grid : str or xarray.DataArray
             The file name of the input grid or the grid loaded as a DataArray.
+        {coord_sys}
         {r}
         """
         kwargs = self._preprocess(**kwargs)
@@ -315,7 +317,7 @@ def grdimage(self, grid, registration=None, **kwargs):
             if kind == "file":
                 file_context = dummy_context(grid)
             elif kind == "grid":
-                file_context = lib.virtualfile_from_grid(grid, registration)
+                file_context = lib.virtualfile_from_grid(grid, coord_sys, registration)
             else:
                 raise GMTInvalidInput("Unrecognized data type: {}".format(type(grid)))
             with file_context as fname:
@@ -339,7 +341,7 @@ def grdimage(self, grid, registration=None, **kwargs):
         p="perspective",
     )
     @kwargs_to_strings(R="sequence", p="sequence")
-    def grdview(self, grid, registration=None, **kwargs):
+    def grdview(self, grid, coord_sys=None, registration=None, **kwargs):
         """
         Create 3-D perspective image or surface mesh from a grid.
 
@@ -403,6 +405,7 @@ def grdview(self, grid, registration=None, **kwargs):
         perspective : list or str
             ``'[x|y|z]azim[/elev[/zlevel]][+wlon0/lat0[/z0]][+vx0/y0]'``.
             Select perspective view.
+        {coord_sys}
         {r}
         """
         kwargs = self._preprocess(**kwargs)
@@ -411,7 +414,7 @@ def grdview(self, grid, registration=None, **kwargs):
             if kind == "file":
                 file_context = dummy_context(grid)
             elif kind == "grid":
-                file_context = lib.virtualfile_from_grid(grid, registration)
+                file_context = lib.virtualfile_from_grid(grid, coord_sys, registration)
             else:
                 raise GMTInvalidInput(f"Unrecognized data type for grid: {type(grid)}")
 
@@ -422,7 +425,7 @@ def grdview(self, grid, registration=None, **kwargs):
                     if data_kind(drapegrid) in ("file", "grid"):
                         if data_kind(drapegrid) == "grid":
                             drape_context = lib.virtualfile_from_grid(
-                                drapegrid, registration
+                                drapegrid, coord_sys, registration
                             )
                             drapefile = stack.enter_context(drape_context)
                             kwargs["G"] = drapefile

pygmt/clib/session.py

@@ -1174,7 +1174,7 @@ def virtualfile_from_matrix(self, matrix):
             yield vfile
 
     @contextmanager
-    def virtualfile_from_grid(self, grid, registration=None):
+    def virtualfile_from_grid(self, grid, coord_sys=None, registration=None):
         """
         Store a grid in a virtual file.
 
@@ -1201,6 +1201,9 @@ def virtualfile_from_grid(self, grid, registration=None):
         ----------
         grid : :class:`xarray.DataArray`
             The grid that will be included in the virtual file.
+        coord_sys : str or None
+            Use a Cartesian (c) or Geographic (g) coordinate system. Default is
+            auto (None), with a fallback to Cartesian (c).
         registration : str or None
             ``[g|p]``
             Use gridline (g) or pixel (p) node registration to make the virtual
@@ -1237,37 +1240,56 @@ def virtualfile_from_grid(self, grid, registration=None):
         >>> # The output is: w e s n z0 z1 dx dy n_columns n_rows
 
         """
-        if registration is None:
+        try:
+            coord_sys_dict = {
+                None: "GMT_GRID_IS_CARTESIAN",  # Default to Cartesian
+                "c": "GMT_GRID_IS_CARTESIAN",
+                "g": "GMT_GRID_IS_GEO",
+            }
+            _coord_sys = coord_sys_dict[coord_sys]
+        except KeyError:
+            raise GMTInvalidInput(
+                "The coord_sys argument should be either 'c', 'g', or None (auto)"
+            )
+
+        try:
+            registration_dict = {
+                None: "GMT_GRID_NODE_REG",  # Default to gridline registration
+                "g": "GMT_GRID_NODE_REG",
+                "p": "GMT_GRID_PIXEL_REG",
+            }
+            _registration = registration_dict[registration]
+        except KeyError:
+            raise GMTInvalidInput(
+                "Invalid registration type, must be either 'g', 'p', or None (auto)"
+            )
+
+        if coord_sys is None or registration is None:
             # Automatically detect whether the NetCDF source of an
-            # xarray.DataArray grid uses gridline or pixel registration.
-            # Defaults to gridline if grdinfo cannot find any source file.
-            registration = "GMT_GRID_NODE_REG"  # default to gridline registration
+            # xarray.DataArray grid uses a Cartesian or Geographic coordinate
+            # system, and whether it is gridline or pixel registered.
             try:
                 gridfile = grid.encoding["source"]
                 with GMTTempFile() as gridinfotext:
                     arg_str = " ".join([gridfile, "->" + gridinfotext.name])
                     self.call_module("grdinfo", arg_str)
-                    if "Pixel node registration used" in gridinfotext.read():
-                        registration = "GMT_GRID_PIXEL_REG"
+                    if coord_sys is None and "[Geographic grid]" in gridinfotext.read():
+                        _coord_sys = "GMT_GRID_IS_GEO"
+                    if (
+                        registration is None
+                        and "Pixel node registration used" in gridinfotext.read()
+                    ):
+                        _registration = "GMT_GRID_PIXEL_REG"
             except KeyError:
                 pass
-        else:
-            if registration == "g":
-                registration = "GMT_GRID_NODE_REG"
-            elif registration == "p":
-                registration = "GMT_GRID_PIXEL_REG"
-            else:
-                raise GMTInvalidInput(
-                    "Invalid registration type, must be either 'g', 'p', or None (auto)"
-                )
 
         # Conversion to a C-contiguous array needs to be done here and not in
         # put_matrix because we need to maintain a reference to the copy while
         # it is being used by the C API. Otherwise, the array would be garbage
         # collected and the memory freed. Creating it in this context manager
         # guarantees that the copy will be around until the virtual file is
         # closed. The conversion is implicit in dataarray_to_matrix.
-        matrix, region, inc = dataarray_to_matrix(grid, registration)
+        matrix, region, inc = dataarray_to_matrix(grid, _registration)
 
         if Version(self.info["version"]) < Version("6.1.0"):
             # Avoid crashing PyGMT because GMT < 6.1.0's C API cannot handle
@@ -1280,10 +1302,10 @@ def virtualfile_from_grid(self, grid, registration=None):
         gmt_grid = self.create_data(
             family,
             geometry,
-            mode="GMT_CONTAINER_ONLY|GMT_GRID_IS_CARTESIAN",
+            mode=f"GMT_CONTAINER_ONLY|{_coord_sys}",
             ranges=region,
             inc=inc,
-            registration=registration,
+            registration=_registration,
         )
         self.put_matrix(gmt_grid, matrix)
         args = (family, geometry, "GMT_IN|GMT_IS_REFERENCE", gmt_grid)

pygmt/helpers/decorators.py

@@ -54,6 +54,11 @@
             ``[g|p]``
             Force gridline (g) or pixel (p) node registration. Default is
             gridline.""",
+    "coord_sys": """\
+        coord_sys : bool
+            Coordinate System of grid. Can be either Cartesian (c) or
+            Geographic (g). Default is auto-detect (None), with a fallback
+            to Cartesian (c).""",
 }
 
 

pygmt/modules.py

@@ -14,7 +14,7 @@
 
 
 @fmt_docstring
-def grdinfo(grid, registration=None, **kwargs):
+def grdinfo(grid, coord_sys=None, registration=None, **kwargs):
     """
     Get information about a grid.
 
@@ -27,6 +27,7 @@ def grdinfo(grid, registration=None, **kwargs):
     grid : str or xarray.DataArray
         The file name of the input grid or the grid loaded as a DataArray.
 
+    {coord_sys}
     {r}
 
     Returns
@@ -41,7 +42,7 @@ def grdinfo(grid, registration=None, **kwargs):
             if kind == "file":
                 file_context = dummy_context(grid)
             elif kind == "grid":
-                file_context = lib.virtualfile_from_grid(grid, registration)
+                file_context = lib.virtualfile_from_grid(grid, coord_sys, registration)
             else:
                 raise GMTInvalidInput("Unrecognized data type: {}".format(type(grid)))
             with file_context as infile:

pygmt/sampling.py

@@ -17,7 +17,15 @@
 
 @fmt_docstring
 @use_alias(n="interpolation")
-def grdtrack(points, grid, newcolname=None, outfile=None, registration=None, **kwargs):
+def grdtrack(
+    points,
+    grid,
+    newcolname=None,
+    outfile=None,
+    registration=None,
+    coord_sys=None,
+    **kwargs,
+):
     """
     Sample grids at specified (x,y) locations.
 
@@ -55,6 +63,7 @@ def grdtrack(points, grid, newcolname=None, outfile=None, registration=None, **k
         Required if 'points' is a file. The file name for the output ASCII
         file.
 
+    {coord_sys}
     {r}
     {n}
 
@@ -85,7 +94,7 @@ def grdtrack(points, grid, newcolname=None, outfile=None, registration=None, **k
 
             # Store the xarray.DataArray grid in virtualfile
             if data_kind(grid) == "grid":
-                grid_context = lib.virtualfile_from_grid(grid, registration)
+                grid_context = lib.virtualfile_from_grid(grid, coord_sys, registration)
             elif data_kind(grid) == "file":
                 grid_context = dummy_context(grid)
             else:

pygmt/tests/test_grdimage.py

@@ -2,13 +2,35 @@
 Test Figure.grdimage
 """
 import numpy as np
+import xarray as xr
 import pytest
 
 from .. import Figure
 from ..exceptions import GMTInvalidInput
 from ..datasets import load_earth_relief
 
 
+@pytest.fixture(scope="module")
+def xrgrid():
+    """
+    Create a sample xarray.DataArray grid for testing
+    """
+    longitude = np.arange(0, 360, 1)
+    latitude = np.arange(-89, 91, 1)
+    x = np.sin(np.deg2rad(longitude))
+    y = np.linspace(start=0, stop=1, num=180)
+    data = y[:, np.newaxis] * x
+
+    return xr.DataArray(
+        data,
+        coords=[
+            ("latitude", latitude, {"units": "degrees_north"}),
+            ("longitude", longitude, {"units": "degrees_east"}),
+        ],
+        attrs={"actual_range": [-1, 1]},
+    )
+
+
 @pytest.mark.mpl_image_compare
 def test_grdimage():
     "Plot an image using an xarray grid"
@@ -46,3 +68,17 @@ def test_grdimage_fails():
     fig = Figure()
     with pytest.raises(GMTInvalidInput):
         fig.grdimage(np.arange(20).reshape((4, 5)))
+
+
+@pytest.mark.mpl_image_compare
+def test_gridimage_over_dateline(xrgrid):
+    """
+    Ensure no gaps are plotted over the 180 degree international dateline.
+    Specifically checking that coord_sys="g" sets `GMT_GRID_IS_GEO`, and that
+    registration="p" sets `GMT_GRID_PIXEL_REG`.
+    """
+    fig = Figure()
+    fig.grdimage(
+        grid=xrgrid, region="g", projection="A0/0/1i", coord_sys="g", registration="p"
+    )
+    return fig