Properly allow for either pixel or gridline registered grids

pygmt/clib/conversion.py

@@ -49,7 +49,7 @@ def dataarray_to_matrix(grid):
     >>> grid = load_earth_relief(resolution='01d')
     >>> matrix, region, inc = dataarray_to_matrix(grid)
     >>> print(region)
-    [-179.5, 179.5, -89.5, 89.5]
+    [-180.0, 180.0, -90.0, 90.0]
     >>> print(inc)
     [1.0, 1.0]
     >>> type(matrix)
@@ -66,7 +66,7 @@ def dataarray_to_matrix(grid):
     >>> print(matrix.shape)
     (31, 71)
     >>> print(region)
-    [-149.5, -79.5, -79.5, -49.5]
+    [-150.0, -79.0, -80.0, -49.0]
     >>> print(inc)
     [1.0, 1.0]
     >>> # but not if only taking every other grid point.
@@ -76,7 +76,7 @@ def dataarray_to_matrix(grid):
     >>> print(matrix.shape)
     (16, 36)
     >>> print(region)
-    [-149.5, -79.5, -79.5, -49.5]
+    [-150.5, -78.5, -80.5, -48.5]
     >>> print(inc)
     [2.0, 2.0]
 
@@ -101,14 +101,19 @@ def dataarray_to_matrix(grid):
                     dim
                 )
             )
-        region.extend([coord.min(), coord.max()])
+        region.extend(
+            [
+                coord.min() - coord_inc / 2 * grid.gmt.registration,
+                coord.max() + coord_inc / 2 * grid.gmt.registration,
+            ]
+        )
         inc.append(coord_inc)
 
-    if any([i < 0 for i in inc]):  # Sort grid when there are negative increments
+    if any(i < 0 for i in inc):  # Sort grid when there are negative increments
         inc = [abs(i) for i in inc]
         grid = grid.sortby(variables=list(grid.dims), ascending=True)
 
-    matrix = as_c_contiguous(grid.values[::-1])
+    matrix = as_c_contiguous(grid[::-1].values)
     return matrix, region, inc
 
 

pygmt/clib/session.py

@@ -43,7 +43,7 @@
     "GMT_IS_SURFACE",
 ]
 
-MODES = ["GMT_CONTAINER_ONLY", "GMT_OUTPUT"]
+MODES = ["GMT_CONTAINER_ONLY", "GMT_IS_OUTPUT"]
 
 REGISTRATIONS = ["GMT_GRID_PIXEL_REG", "GMT_GRID_NODE_REG"]
 
@@ -111,7 +111,7 @@ class Session:
     ...             )
     ...             # Read the contents of the temp file before it's deleted.
     ...             print(fout.read().strip())
-    -179.5 179.5 -89.5 89.5 -8182 5651.5 1 1 360 180 0 0
+    -180 180 -90 90 -8182 5651.5 1 1 360 180 1 1
     """
 
     # The minimum version of GMT required
@@ -509,13 +509,13 @@ def create_data(self, family, geometry, mode, **kwargs):
         ----------
         family : str
             A valid GMT data family name (e.g., ``'GMT_IS_DATASET'``). See the
-            ``data_families`` attribute for valid names.
+            ``FAMILIES`` attribute for valid names.
         geometry : str
             A valid GMT data geometry name (e.g., ``'GMT_IS_POINT'``). See the
-            ``data_geometries`` attribute for valid names.
+            ``GEOMETRIES`` attribute for valid names.
         mode : str
-            A valid GMT data mode (e.g., ``'GMT_OUTPUT'``). See the
-            ``data_modes`` attribute for valid names.
+            A valid GMT data mode (e.g., ``'GMT_IS_OUTPUT'``). See the
+            ``MODES`` attribute for valid names.
         dim : list of 4 integers
             The dimensions of the dataset. See the documentation for the GMT C
             API function ``GMT_Create_Data`` (``src/gmt_api.c``) for the full
@@ -528,7 +528,7 @@ def create_data(self, family, geometry, mode, **kwargs):
         inc : list of 2 floats
             The increments between points of the dataset. See the C function
             documentation.
-        registration : int
+        registration : str
             The node registration (what the coordinates mean). Can be
             ``'GMT_GRID_PIXEL_REG'`` or ``'GMT_GRID_NODE_REG'``. Defaults to
             ``'GMT_GRID_NODE_REG'``.
@@ -561,7 +561,9 @@ def create_data(self, family, geometry, mode, **kwargs):
 
         family_int = self._parse_constant(family, valid=FAMILIES, valid_modifiers=VIAS)
         mode_int = self._parse_constant(
-            mode, valid=MODES, valid_modifiers=["GMT_GRID_IS_GEO"]
+            mode,
+            valid=MODES,
+            valid_modifiers=["GMT_GRID_IS_CARTESIAN", "GMT_GRID_IS_GEO"],
         )
         geometry_int = self._parse_constant(geometry, valid=GEOMETRIES)
         registration_int = self._parse_constant(
@@ -841,12 +843,12 @@ def write_data(self, family, geometry, mode, wesn, output, data):
         ----------
         family : str
             A valid GMT data family name (e.g., ``'GMT_IS_DATASET'``). See the
-            ``data_families`` attribute for valid names. Don't use the
+            ``FAMILIES`` attribute for valid names. Don't use the
             ``GMT_VIA_VECTOR`` or ``GMT_VIA_MATRIX`` constructs for this. Use
             ``GMT_IS_VECTOR`` and ``GMT_IS_MATRIX`` instead.
         geometry : str
             A valid GMT data geometry name (e.g., ``'GMT_IS_POINT'``). See the
-            ``data_geometries`` attribute for valid names.
+            ``GEOMETRIES`` attribute for valid names.
         mode : str
             How the data is to be written to the file. This option varies
             depending on the given family. See the GMT API documentation for
@@ -1228,21 +1230,30 @@ def virtualfile_from_grid(self, grid):
         ...             args = '{} -L0 -Cn ->{}'.format(fin, fout.name)
         ...             ses.call_module('grdinfo', args)
         ...             print(fout.read().strip())
-        -179.5 179.5 -89.5 89.5 -8182 5651.5 1 1 360 180 0 0
+        -180 180 -90 90 -8182 5651.5 1 1 360 180 1 1
         >>> # The output is: w e s n z0 z1 dx dy n_columns n_rows reg gtype
 
         """
+        _gtype = {0: "GMT_GRID_IS_CARTESIAN", 1: "GMT_GRID_IS_GEO"}[grid.gmt.gtype]
+        _reg = {0: "GMT_GRID_NODE_REG", 1: "GMT_GRID_PIXEL_REG"}[grid.gmt.registration]
+
         # 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)
+
         family = "GMT_IS_GRID|GMT_VIA_MATRIX"
         geometry = "GMT_IS_SURFACE"
         gmt_grid = self.create_data(
-            family, geometry, mode="GMT_CONTAINER_ONLY", ranges=region, inc=inc
+            family,
+            geometry,
+            mode=f"GMT_CONTAINER_ONLY|{_gtype}",
+            ranges=region,
+            inc=inc,
+            registration=_reg,
         )
         self.put_matrix(gmt_grid, matrix)
         args = (family, geometry, "GMT_IN|GMT_IS_REFERENCE", gmt_grid)

pygmt/helpers/decorators.py

@@ -49,6 +49,11 @@
             - 'c' for bicubic [Default]
             - 'l' for bilinear
             - 'n' for nearest-neighbor""",
+    "registration": """\
+        registration : str
+            ``[g|p]``
+            Force output grid to be gridline (g) or pixel (p) node registered.
+            Default is gridline (g).""",
 }
 
 

pygmt/tests/test_grdimage.py

@@ -2,6 +2,7 @@
 Test Figure.grdimage
 """
 import numpy as np
+import xarray as xr
 import pytest
 
 from .. import Figure
@@ -15,6 +16,27 @@ def fixture_grid():
     return load_earth_relief(registration="gridline")
 
 
+@pytest.fixture(scope="module", name="xrgrid")
+def fixture_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(grid):
     "Plot an image using an xarray grid"
@@ -51,3 +73,18 @@ 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_grdimage_over_dateline(xrgrid):
+    """
+    Ensure no gaps are plotted over the 180 degree international dateline.
+    Specifically checking that `xrgrid.gmt.gtype = 1` sets `GMT_GRID_IS_GEO`,
+    and that `xrgrid.gmt.registration = 0` sets `GMT_GRID_NODE_REG`. Note that
+    there would be a gap over the dateline if a pixel registered grid is used.
+    """
+    fig = Figure()
+    assert xrgrid.gmt.registration == 0  # gridline registration
+    xrgrid.gmt.gtype = 1  # geographic coordinate system
+    fig.grdimage(grid=xrgrid, region="g", projection="A0/0/1i", V="d")
+    return fig

pygmt/tests/test_grdinfo.py

@@ -13,7 +13,7 @@ def test_grdinfo():
     "Make sure grd info works as expected"
     grid = load_earth_relief(registration="gridline")
     result = grdinfo(grid, L=0, C="n")
-    assert result.strip() == "-180 180 -90 90 -8592.5 5559 1 1 361 181 0 0"
+    assert result.strip() == "-180 180 -90 90 -8592.5 5559 1 1 361 181 0 1"
 
 
 def test_grdinfo_file():

pygmt/tests/test_grdview.py

@@ -17,6 +17,9 @@ def fixture_grid():
     )
 
 
+@pytest.mark.xfail(
+    reason="Baseline image generated using Cartesian instead of Geographic coordinates"
+)
 @pytest.mark.mpl_image_compare
 def test_grdview_grid_dataarray(grid):
     """
@@ -54,6 +57,9 @@ def test_grdview_wrong_kind_of_grid(grid):
         fig.grdview(grid=dataset)
 
 
+@pytest.mark.xfail(
+    reason="Baseline image generated using Cartesian instead of Geographic coordinates"
+)
 @pytest.mark.mpl_image_compare
 def test_grdview_with_perspective(grid):
     """