diff --git a/Makefile b/Makefile index 3f643a7d2ad..cb27e1d9a3d 100644 --- a/Makefile +++ b/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) diff --git a/pygmt/clib/conversion.py b/pygmt/clib/conversion.py index 4334eb14daa..8dd8adbeb93 100644 --- a/pygmt/clib/conversion.py +++ b/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 diff --git a/pygmt/clib/session.py b/pygmt/clib/session.py index d6ba299090c..a2b8ef1e355 100644 --- a/pygmt/clib/session.py +++ b/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,10 +1245,13 @@ 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 @@ -1254,10 +1259,16 @@ def virtualfile_from_grid(self, grid): # 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) diff --git a/pygmt/helpers/decorators.py b/pygmt/helpers/decorators.py index 0e7933ea36c..518717851cf 100644 --- a/pygmt/helpers/decorators.py +++ b/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).""", } diff --git a/pygmt/tests/baseline/test_grdimage_over_dateline.png b/pygmt/tests/baseline/test_grdimage_over_dateline.png new file mode 100644 index 00000000000..34b909c2beb Binary files /dev/null and b/pygmt/tests/baseline/test_grdimage_over_dateline.png differ diff --git a/pygmt/tests/test_grdimage.py b/pygmt/tests/test_grdimage.py index 1db1428666c..45be76f6d01 100644 --- a/pygmt/tests/test_grdimage.py +++ b/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 diff --git a/pygmt/tests/test_grdinfo.py b/pygmt/tests/test_grdinfo.py index f27768e4797..6991000b668 100644 --- a/pygmt/tests/test_grdinfo.py +++ b/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(): diff --git a/pygmt/tests/test_grdview.py b/pygmt/tests/test_grdview.py index 125f200e467..65af6e6eef3 100644 --- a/pygmt/tests/test_grdview.py +++ b/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): """ diff --git a/setup.cfg b/setup.cfg index d607b9e76fd..812f2e81962 100644 --- a/setup.cfg +++ b/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