diff --git a/pygmt/base_plotting.py b/pygmt/base_plotting.py index 85aae6425fc..427c0855bee 100644 --- a/pygmt/base_plotting.py +++ b/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 diff --git a/pygmt/clib/session.py b/pygmt/clib/session.py index 4eef5dc707e..aee21e7eff4 100644 --- a/pygmt/clib/session.py +++ b/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,29 +1240,48 @@ 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 @@ -1267,7 +1289,7 @@ def virtualfile_from_grid(self, grid, registration=None): # 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) diff --git a/pygmt/helpers/decorators.py b/pygmt/helpers/decorators.py index 15b13aa65f4..c5ecb200949 100644 --- a/pygmt/helpers/decorators.py +++ b/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).""", } diff --git a/pygmt/modules.py b/pygmt/modules.py index a63f4b492ba..52468d0edd5 100644 --- a/pygmt/modules.py +++ b/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: diff --git a/pygmt/sampling.py b/pygmt/sampling.py index 4cc27fe6651..26503c5e07d 100644 --- a/pygmt/sampling.py +++ b/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: diff --git a/pygmt/tests/baseline/test_gridimage_over_dateline.png b/pygmt/tests/baseline/test_gridimage_over_dateline.png new file mode 100644 index 00000000000..94e57898801 Binary files /dev/null and b/pygmt/tests/baseline/test_gridimage_over_dateline.png differ diff --git a/pygmt/tests/test_grdimage.py b/pygmt/tests/test_grdimage.py index f74499c54d1..e970c398026 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 @@ -9,6 +10,27 @@ 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