diff --git a/README.rst b/README.rst index 35e0a82898d..e005bbce80a 100644 --- a/README.rst +++ b/README.rst @@ -45,7 +45,7 @@ Disclaimer All of the API (functions/classes/interfaces) is subject to change until we reach v1.0.0 as per the `semantic versioning specification `__. There may be non-backward compatible changes as we experiment with new design ideas and -implement new features. **This is not a finished product, use with caution** +implement new features. **This is not a finished product, use with caution.** We welcome any feedback and ideas! Let us know by submitting @@ -85,9 +85,6 @@ Contacting Us open issue or pull request. * We have a `Discourse forum `__ where you can ask questions and leave comments. -* This project is released with a `Contributor Code of Conduct - `__. - By participating in this project you agree to abide by its terms. Contributing diff --git a/examples/tutorials/first-figure.py b/examples/tutorials/first-figure.py index 44d48f86d2b..752f010ab5e 100644 --- a/examples/tutorials/first-figure.py +++ b/examples/tutorials/first-figure.py @@ -66,8 +66,7 @@ # # You’ll probably have noticed several things that are different from classic # command-line GMT. Many of these changes reflect the new GMT modern execution mode that -# will be part of the future 6.0 release. A few are PyGMT exclusive (like the -# ``savefig`` method). +# are part of GMT 6. A few are PyGMT exclusive (like the ``savefig`` method). # # 1. The name of method is ``coast`` instead of ``pscoast``. As a general rule, all # ``ps*`` modules had their ``ps`` prefix removed. The exceptions are: diff --git a/pygmt/modules.py b/pygmt/modules.py index 94512ec325d..1cacf67d6df 100644 --- a/pygmt/modules.py +++ b/pygmt/modules.py @@ -60,14 +60,17 @@ def info(table, **kwargs): """ Get information about data tables. - Reads from files and finds the extreme values in each of the columns. - It recognizes NaNs and will print warnings if the number of columns vary - from record to record. As an option, it will find the extent of the first - n columns rounded up and down to the nearest multiple of the supplied - increments. By default, this output will be in the form *-Rw/e/s/n*, - or the output will be in column form for as many columns as there are - increments provided. The *nearest_multiple* option will provide a - *-Tzmin/zmax/dz* string for makecpt. + Reads from files and finds the extreme values in each of the columns + reported as min/max pairs. It recognizes NaNs and will print warnings if + the number of columns vary from record to record. As an option, it will + find the extent of the first two columns rounded up and down to the nearest + multiple of the supplied increments given by *spacing*. Such output will be + in a numpy.ndarray form ``[w, e, s, n]``, which can be used directly as the + *region* argument for other modules (hence only dx and dy are needed). If + the *per_column* option is combined with *spacing*, then the numpy.ndarray + output will be rounded up/down for as many columns as there are increments + provided in *spacing*. A similar option *nearest_multiple* option will + provide a numpy.ndarray in the form of ``[zmin, zmax, dz]`` for makecpt. Full option list at :gmt-docs:`gmtinfo.html` @@ -83,12 +86,21 @@ def info(table, **kwargs): spacing : str ``'[b|p|f|s]dx[/dy[/dz...]]'``. Report the min/max of the first n columns to the nearest multiple of - the provided increments and output results in the form *-Rw/e/s/n* - (unless *per_column* is set). + the provided increments and output results in the form + ``[w, e, s, n]``. nearest_multiple : str ``'dz[+ccol]'`` Report the min/max of the first (0'th) column to the nearest multiple - of dz and output this as the string *-Tzmin/zmax/dz*. + of dz and output this in the form ``[zmin, zmax, dz]``. + + Returns + ------- + output : np.ndarray or str + Return type depends on whether any of the 'per_column', 'spacing', or + 'nearest_multiple' parameters are set. + + - np.ndarray if either of the above parameters are used. + - str if none of the above parameters are used. """ kind = data_kind(table) with Session() as lib: @@ -108,7 +120,16 @@ def info(table, **kwargs): [fname, build_arg_string(kwargs), "->" + tmpfile.name] ) lib.call_module("info", arg_str) - return tmpfile.read() + result = tmpfile.read() + + if any(arg in kwargs for arg in ["C", "I", "T"]): + # Converts certain output types into a numpy array + # instead of a raw string that is less useful. + if result.startswith(("-R", "-T")): # e.g. -R0/1/2/3 or -T0/9/1 + result = result[2:].replace("/", " ") + result = np.loadtxt(result.splitlines()) + + return result @fmt_docstring diff --git a/pygmt/tests/test_grdimage.py b/pygmt/tests/test_grdimage.py index f686af883a4..214de9def42 100644 --- a/pygmt/tests/test_grdimage.py +++ b/pygmt/tests/test_grdimage.py @@ -69,6 +69,27 @@ def test_grdimage_file(): return fig +@pytest.mark.xfail(reason="Upstream bug in GMT 6.1.1") +@check_figures_equal() +def test_grdimage_xarray_shading(grid, fig_ref, fig_test): + """ + Test that shading works well for xarray. + See https://github.com/GenericMappingTools/pygmt/issues/364 + """ + fig_ref, fig_test = Figure(), Figure() + kwargs = dict( + region=[-180, 180, -90, 90], + frame=True, + projection="Cyl_stere/6i", + cmap="geo", + shading=True, + ) + + fig_ref.grdimage("@earth_relief_01d_g", **kwargs) + fig_test.grdimage(grid, **kwargs) + return fig_ref, fig_test + + def test_grdimage_fails(): "Should fail for unrecognized input" fig = Figure() diff --git a/pygmt/tests/test_grdview.py b/pygmt/tests/test_grdview.py index 65af6e6eef3..f9166b57b90 100644 --- a/pygmt/tests/test_grdview.py +++ b/pygmt/tests/test_grdview.py @@ -3,53 +3,64 @@ """ import pytest -from .. import Figure, which -from ..datasets import load_earth_relief +from .. import Figure, grdcut, which from ..exceptions import GMTInvalidInput -from ..helpers import data_kind +from ..helpers import GMTTempFile, data_kind +from ..helpers.testing import check_figures_equal -@pytest.fixture(scope="module", name="grid") -def fixture_grid(): - "Load the grid data from the sample earth_relief file" - return load_earth_relief(registration="gridline").sel( - lat=slice(-49, -42), lon=slice(-118, -107) - ) +@pytest.fixture(scope="module", name="region") +def fixture_region(): + "Test region as lonmin, lonmax, latmin, latmax" + return (-116, -109, -47, -44) -@pytest.mark.xfail( - reason="Baseline image generated using Cartesian instead of Geographic coordinates" -) -@pytest.mark.mpl_image_compare -def test_grdview_grid_dataarray(grid): +@pytest.fixture(scope="module", name="gridfile") +def fixture_gridfile(region): + """ + Load the NetCDF grid file from the sample earth_relief file + """ + with GMTTempFile(suffix=".nc") as tmpfile: + grdcut(grid="@earth_relief_01d_g", region=region, outgrid=tmpfile.name) + yield tmpfile.name + + +@pytest.fixture(scope="module", name="xrgrid") +def fixture_xrgrid(region): + """ + Load the xarray.DataArray grid from the sample earth_relief file + """ + return grdcut(grid="@earth_relief_01d_g", region=region) + + +@check_figures_equal() +def test_grdview_grid_dataarray(gridfile, xrgrid): """ Run grdview by passing in a grid as an xarray.DataArray. """ - fig = Figure() - fig.grdview(grid=grid) - return fig + fig_ref, fig_test = Figure(), Figure() + fig_ref.grdview(grid=gridfile) + fig_test.grdview(grid=xrgrid) + return fig_ref, fig_test -@pytest.mark.xfail( - reason="Baseline image not updated to use earth relief grid in GMT 6.1.0", -) @pytest.mark.mpl_image_compare -def test_grdview_grid_file_with_region_subset(): +def test_grdview_grid_file_with_region_subset(region): """ Run grdview by passing in a grid filename, and cropping it to a region. """ gridfile = which("@earth_relief_01d_g", download="a") fig = Figure() - fig.grdview(grid=gridfile, region=[-116, -109, -47, -44]) + fig.grdview(grid=gridfile, region=region) return fig -def test_grdview_wrong_kind_of_grid(grid): +def test_grdview_wrong_kind_of_grid(xrgrid): """ Run grdview using grid input that is not an xarray.DataArray or file. """ - dataset = grid.to_dataset() # convert xarray.DataArray to xarray.Dataset + dataset = xrgrid.to_dataset() # convert xarray.DataArray to xarray.Dataset assert data_kind(dataset) == "matrix" fig = Figure() @@ -57,208 +68,193 @@ 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): +@check_figures_equal() +def test_grdview_with_perspective(gridfile, xrgrid): """ Run grdview by passing in a grid and setting a perspective viewpoint with an azimuth from the SouthEast and an elevation angle 15 degrees from the z-plane. """ - fig = Figure() - fig.grdview(grid=grid, perspective=[135, 15]) - return fig + fig_ref, fig_test = Figure(), Figure() + fig_ref.grdview(grid=gridfile, perspective=[135, 15]) + fig_test.grdview(grid=xrgrid, perspective=[135, 15]) + return fig_ref, fig_test -@pytest.mark.xfail( - reason="Baseline image not updated to use earth relief grid in GMT 6.1.0", -) -@pytest.mark.mpl_image_compare -def test_grdview_with_perspective_and_zscale(grid): +@check_figures_equal() +def test_grdview_with_perspective_and_zscale(gridfile, xrgrid): """ Run grdview by passing in a grid and setting a perspective viewpoint with an azimuth from the SouthWest and an elevation angle 30 degrees from the z-plane, plus a z-axis scaling factor of 0.005. """ - fig = Figure() - fig.grdview(grid=grid, perspective=[225, 30], zscale=0.005) - return fig + fig_ref, fig_test = Figure(), Figure() + kwargs = dict(perspective=[225, 30], zscale=0.005) + fig_ref.grdview(grid=gridfile, **kwargs) + fig_test.grdview(grid=xrgrid, **kwargs) + return fig_ref, fig_test -@pytest.mark.xfail( - reason="Baseline image not updated to use earth relief grid in GMT 6.1.0", -) -@pytest.mark.mpl_image_compare -def test_grdview_with_perspective_and_zsize(grid): +@check_figures_equal() +def test_grdview_with_perspective_and_zsize(gridfile, xrgrid): """ Run grdview by passing in a grid and setting a perspective viewpoint with an azimuth from the SouthWest and an elevation angle 30 degrees from the z-plane, plus a z-axis size of 10cm. """ - fig = Figure() - fig.grdview(grid=grid, perspective=[225, 30], zsize="10c") - return fig + fig_ref, fig_test = Figure(), Figure() + kwargs = dict(perspective=[225, 30], zsize="10c") + fig_ref.grdview(grid=gridfile, **kwargs) + fig_test.grdview(grid=xrgrid, **kwargs) + return fig_ref, fig_test -@pytest.mark.xfail( - reason="Baseline image not updated to use earth relief grid in GMT 6.1.0", -) -@pytest.mark.mpl_image_compare -def test_grdview_with_cmap_for_image_plot(grid): +@check_figures_equal() +def test_grdview_with_cmap_for_image_plot(gridfile, xrgrid): """ Run grdview by passing in a grid and setting a colormap for producing an image plot. """ - fig = Figure() - fig.grdview(grid=grid, cmap="oleron", surftype="i") - return fig + fig_ref, fig_test = Figure(), Figure() + kwargs = dict(cmap="oleron", surftype="i") + fig_ref.grdview(grid=gridfile, **kwargs) + fig_test.grdview(grid=xrgrid, **kwargs) + return fig_ref, fig_test -@pytest.mark.xfail( - reason="Baseline image not updated to use earth relief grid in GMT 6.1.0", -) -@pytest.mark.mpl_image_compare -def test_grdview_with_cmap_for_surface_monochrome_plot(grid): +@check_figures_equal() +def test_grdview_with_cmap_for_surface_monochrome_plot(gridfile, xrgrid): """ Run grdview by passing in a grid and setting a colormap for producing a surface monochrome plot. """ - fig = Figure() - fig.grdview(grid=grid, cmap="oleron", surftype="s+m") - return fig + fig_ref, fig_test = Figure(), Figure() + kwargs = dict(cmap="oleron", surftype="s+m") + fig_ref.grdview(grid=gridfile, **kwargs) + fig_test.grdview(grid=xrgrid, **kwargs) + return fig_ref, fig_test -@pytest.mark.xfail( - reason="Baseline image not updated to use earth relief grid in GMT 6.1.0", -) -@pytest.mark.mpl_image_compare -def test_grdview_with_cmap_for_perspective_surface_plot(grid): +@check_figures_equal() +def test_grdview_with_cmap_for_perspective_surface_plot(gridfile, xrgrid): """ Run grdview by passing in a grid and setting a colormap for producing a surface plot with a 3D perspective viewpoint. """ - fig = Figure() - fig.grdview( - grid=grid, cmap="oleron", surftype="s", perspective=[225, 30], zscale=0.005 - ) - return fig + fig_ref, fig_test = Figure(), Figure() + kwargs = dict(cmap="oleron", surftype="s", perspective=[225, 30], zscale=0.005) + fig_ref.grdview(grid=gridfile, **kwargs) + fig_test.grdview(grid=xrgrid, **kwargs) + return fig_ref, fig_test -@pytest.mark.xfail( - reason="Baseline image not updated to use earth relief grid in GMT 6.1.0", -) -@pytest.mark.mpl_image_compare -def test_grdview_on_a_plane(grid): +@check_figures_equal() +def test_grdview_on_a_plane(gridfile, xrgrid): """ Run grdview by passing in a grid and plotting it on a z-plane, while setting a 3D perspective viewpoint. """ - fig = Figure() - fig.grdview(grid=grid, plane=-4000, perspective=[225, 30], zscale=0.005) - return fig + fig_ref, fig_test = Figure(), Figure() + kwargs = dict(plane=-4000, perspective=[225, 30], zscale=0.005) + fig_ref.grdview(grid=gridfile, **kwargs) + fig_test.grdview(grid=xrgrid, **kwargs) + return fig_ref, fig_test -@pytest.mark.xfail( - reason="Baseline image not updated to use earth relief grid in GMT 6.1.0", -) -@pytest.mark.mpl_image_compare -def test_grdview_on_a_plane_with_colored_frontal_facade(grid): +@check_figures_equal() +def test_grdview_on_a_plane_with_colored_frontal_facade(gridfile, xrgrid): """ Run grdview by passing in a grid and plotting it on a z-plane whose frontal facade is colored gray, while setting a 3D perspective viewpoint. """ - fig = Figure() - fig.grdview(grid=grid, plane="-4000+ggray", perspective=[225, 30], zscale=0.005) - return fig + fig_ref, fig_test = Figure(), Figure() + kwargs = dict(plane="-4000+ggray", perspective=[225, 30], zscale=0.005) + fig_ref.grdview(grid=gridfile, **kwargs) + fig_test.grdview(grid=xrgrid, **kwargs) + return fig_ref, fig_test -@pytest.mark.xfail( - reason="Baseline image not updated to use earth relief grid in GMT 6.1.0", -) -@pytest.mark.mpl_image_compare -def test_grdview_with_perspective_and_zaxis_frame(grid): +@check_figures_equal() +def test_grdview_with_perspective_and_zaxis_frame(gridfile, xrgrid, region): """ Run grdview by passing in a grid and plotting an annotated vertical - z-axis frame. + z-axis frame on a Transverse Mercator (T) projection. """ - fig = Figure() - fig.grdview(grid=grid, perspective=[225, 30], zscale=0.005, frame="zaf") - return fig + fig_ref, fig_test = Figure(), Figure() + projection = f"T{(region[0]+region[1])/2}/{abs((region[2]+region[3])/2)}" + kwargs = dict( + projection=projection, + perspective=[225, 30], + zscale=0.005, + frame=["xaf", "yaf", "zaf"], + ) + fig_ref.grdview(grid=gridfile, **kwargs) + fig_test.grdview(grid=xrgrid, **kwargs) + return fig_ref, fig_test -@pytest.mark.xfail( - reason="Baseline image not updated to use earth relief grid in GMT 6.1.0", -) -@pytest.mark.mpl_image_compare -def test_grdview_surface_plot_styled_with_contourpen(grid): +@check_figures_equal() +def test_grdview_surface_plot_styled_with_contourpen(gridfile, xrgrid): """ Run grdview by passing in a grid with styled contour lines plotted on top of a surface plot. """ - fig = Figure() - fig.grdview(grid=grid, cmap="relief", surftype="s", contourpen="0.5p,black,dash") - return fig + fig_ref, fig_test = Figure(), Figure() + kwargs = dict(cmap="relief", surftype="s", contourpen="0.5p,black,dash") + fig_ref.grdview(grid=gridfile, **kwargs) + fig_test.grdview(grid=xrgrid, **kwargs) + return fig_ref, fig_test -@pytest.mark.xfail( - reason="Baseline image not updated to use earth relief grid in GMT 6.1.0", -) -@pytest.mark.mpl_image_compare -def test_grdview_surface_mesh_plot_styled_with_meshpen(grid): +@check_figures_equal() +def test_grdview_surface_mesh_plot_styled_with_meshpen(gridfile, xrgrid): """ Run grdview by passing in a grid with styled mesh lines plotted on top of a surface mesh plot. """ - fig = Figure() - fig.grdview(grid=grid, cmap="relief", surftype="sm", meshpen="0.5p,black,dash") - return fig + fig_ref, fig_test = Figure(), Figure() + kwargs = dict(cmap="relief", surftype="sm", meshpen="0.5p,black,dash") + fig_ref.grdview(grid=gridfile, **kwargs) + fig_test.grdview(grid=xrgrid, **kwargs) + return fig_ref, fig_test -@pytest.mark.xfail( - reason="Baseline image not updated to use earth relief grid in GMT 6.1.0", -) -@pytest.mark.mpl_image_compare -def test_grdview_on_a_plane_styled_with_facadepen(grid): +@check_figures_equal() +def test_grdview_on_a_plane_styled_with_facadepen(gridfile, xrgrid): """ Run grdview by passing in a grid and plotting it on a z-plane with styled lines for the frontal facade. """ - fig = Figure() - fig.grdview( - grid=grid, - plane=-4000, - perspective=[225, 30], - zscale=0.005, - facadepen="0.5p,blue,dash", + fig_ref, fig_test = Figure(), Figure() + kwargs = dict( + plane=-4000, perspective=[225, 30], zscale=0.005, facadepen="0.5p,blue,dash" ) - return fig + fig_ref.grdview(grid=gridfile, **kwargs) + fig_test.grdview(grid=xrgrid, **kwargs) + return fig_ref, fig_test -@pytest.mark.xfail( - reason="Baseline image not updated to use earth relief grid in GMT 6.1.0", -) -@pytest.mark.mpl_image_compare -def test_grdview_drapegrid_dataarray(grid): +@check_figures_equal() +def test_grdview_drapegrid_dataarray(gridfile, xrgrid): """ Run grdview by passing in both a grid and drapegrid as an xarray.DataArray, setting a colormap for producing an image plot. """ - drapegrid = 1.1 * grid + drapegrid = 1.1 * xrgrid - fig = Figure() - fig.grdview(grid=grid, drapegrid=drapegrid, cmap="oleron", surftype="c") - return fig + fig_ref, fig_test = Figure(), Figure() + fig_ref.grdview(grid=gridfile, drapegrid=drapegrid, cmap="oleron", surftype="c") + fig_test.grdview(grid=xrgrid, drapegrid=drapegrid, cmap="oleron", surftype="c") + return fig_ref, fig_test -def test_grdview_wrong_kind_of_drapegrid(grid): +def test_grdview_wrong_kind_of_drapegrid(xrgrid): """ Run grdview using drapegrid input that is not an xarray.DataArray or file. """ - dataset = grid.to_dataset() # convert xarray.DataArray to xarray.Dataset + dataset = xrgrid.to_dataset() # convert xarray.DataArray to xarray.Dataset assert data_kind(dataset) == "matrix" fig = Figure() with pytest.raises(GMTInvalidInput): - fig.grdview(grid=grid, drapegrid=dataset) + fig.grdview(grid=xrgrid, drapegrid=dataset) diff --git a/pygmt/tests/test_info.py b/pygmt/tests/test_info.py index b7eadc53649..142b56dce8e 100644 --- a/pygmt/tests/test_info.py +++ b/pygmt/tests/test_info.py @@ -4,6 +4,7 @@ import os import numpy as np +import numpy.testing as npt import pandas as pd import pytest import xarray as xr @@ -57,25 +58,42 @@ def test_info_1d_array(): def test_info_per_column(): "Make sure the per_column option works" output = info(table=POINTS_DATA, per_column=True) - assert output == "11.5309 61.7074 -2.9289 7.8648 0.1412 0.9338\n" + npt.assert_allclose( + actual=output, desired=[11.5309, 61.7074, -2.9289, 7.8648, 0.1412, 0.9338] + ) def test_info_spacing(): "Make sure the spacing option works" output = info(table=POINTS_DATA, spacing=0.1) - assert output == "-R11.5/61.8/-3/7.9\n" + npt.assert_allclose(actual=output, desired=[11.5, 61.8, -3, 7.9]) + + +def test_info_spacing_bounding_box(): + "Make sure the spacing option for writing a bounding box works" + output = info(table=POINTS_DATA, spacing="b") + npt.assert_allclose( + actual=output, + desired=[ + [11.5309, -2.9289], + [61.7074, -2.9289], + [61.7074, 7.8648], + [11.5309, 7.8648], + [11.5309, -2.9289], + ], + ) def test_info_per_column_spacing(): "Make sure the per_column and spacing options work together" output = info(table=POINTS_DATA, per_column=True, spacing=0.1) - assert output == "11.5 61.8 -3 7.9 0.1412 0.9338\n" + npt.assert_allclose(actual=output, desired=[11.5, 61.8, -3, 7.9, 0.1412, 0.9338]) def test_info_nearest_multiple(): "Make sure the nearest_multiple option works" output = info(table=POINTS_DATA, nearest_multiple=0.1) - assert output == "-T11.5/61.8/0.1\n" + npt.assert_allclose(actual=output, desired=[11.5, 61.8, 0.1]) def test_info_fails(): diff --git a/pygmt/tests/test_x2sys_cross.py b/pygmt/tests/test_x2sys_cross.py index 5d3ced83aa2..50577b738df 100644 --- a/pygmt/tests/test_x2sys_cross.py +++ b/pygmt/tests/test_x2sys_cross.py @@ -31,7 +31,8 @@ def fixture_tracks(): Load track data from the sample bathymetry file """ dataframe = load_sample_bathymetry() - return [dataframe.query(expr="bathymetry > -20")] # reduce size of dataset + dataframe.columns = ["x", "y", "z"] # longitude, latitude, bathymetry + return [dataframe.query(expr="z > -20")] # reduce size of dataset def test_x2sys_cross_input_file_output_file(mock_x2sys_home): @@ -76,25 +77,57 @@ def test_x2sys_cross_input_file_output_dataframe(mock_x2sys_home): def test_x2sys_cross_input_dataframe_output_dataframe(mock_x2sys_home, tracks): """ Run x2sys_cross by passing in one dataframe, and output external crossovers - to a pandas.DataFrame. Not actually implemented yet, wait for - https://github.com/GenericMappingTools/gmt/issues/3717 + to a pandas.DataFrame. """ with TemporaryDirectory(prefix="X2SYS", dir=os.getcwd()) as tmpdir: tag = os.path.basename(tmpdir) x2sys_init(tag=tag, fmtfile="xyz", force=True) - with pytest.raises(NotImplementedError): - _ = x2sys_cross(tracks=tracks, tag=tag, coe="i", verbose="i") + output = x2sys_cross(tracks=tracks, tag=tag, coe="i", verbose="i") - # assert isinstance(output, pd.DataFrame) - # assert output.shape == (4, 12) - # columns = list(output.columns) - # assert columns[:6] == ["x", "y", "t_1", "t_2", "dist_1", "dist_2"] - # assert columns[6:] == ["head_1","head_2","vel_1","vel_2","z_X","z_M"] - # assert output.dtypes["t_1"].type == np.datetime64 - # assert output.dtypes["t_2"].type == np.datetime64 + assert isinstance(output, pd.DataFrame) + assert output.shape == (14, 12) + columns = list(output.columns) + assert columns[:6] == ["x", "y", "i_1", "i_2", "dist_1", "dist_2"] + assert columns[6:] == ["head_1", "head_2", "vel_1", "vel_2", "z_X", "z_M"] + assert output.dtypes["i_1"].type == np.object_ + assert output.dtypes["i_2"].type == np.object_ + + return output - # return output + +def test_x2sys_cross_input_two_dataframes(mock_x2sys_home): + """ + Run x2sys_cross by passing in two pandas.DataFrame tables with a time + column, and output external crossovers to a pandas.DataFrame + """ + with TemporaryDirectory(prefix="X2SYS", dir=os.getcwd()) as tmpdir: + tag = os.path.basename(tmpdir) + x2sys_init( + tag=tag, fmtfile="xyz", suffix="xyzt", units=["de", "se"], force=True + ) + + # Add a time row to the x2sys fmtfile + with open(file=os.path.join(tmpdir, "xyz.fmt"), mode="a") as fmtfile: + fmtfile.write("time\ta\tN\t0\t1\t0\t%g\n") + + # Create pandas.DataFrame track tables + tracks = [] + for i in range(2): + np.random.seed(seed=i) + track = pd.DataFrame(data=np.random.rand(10, 3), columns=("x", "y", "z")) + track["time"] = pd.date_range(start=f"2020-{i}1-01", periods=10, freq="ms") + tracks.append(track) + + output = x2sys_cross(tracks=tracks, tag=tag, coe="e", verbose="i") + + assert isinstance(output, pd.DataFrame) + assert output.shape == (30, 12) + columns = list(output.columns) + assert columns[:6] == ["x", "y", "t_1", "t_2", "dist_1", "dist_2"] + assert columns[6:] == ["head_1", "head_2", "vel_1", "vel_2", "z_X", "z_M"] + assert output.dtypes["t_1"].type == np.datetime64 + assert output.dtypes["t_2"].type == np.datetime64 def test_x2sys_cross_input_two_filenames(mock_x2sys_home): @@ -131,7 +164,7 @@ def test_x2sys_cross_invalid_tracks_input_type(tracks): Run x2sys_cross using tracks input that is not a pandas.DataFrame (matrix) or str (file) type, which would raise a GMTInvalidInput error. """ - invalid_tracks = tracks[0].to_xarray().bathymetry + invalid_tracks = tracks[0].to_xarray().z assert data_kind(invalid_tracks) == "grid" with pytest.raises(GMTInvalidInput): x2sys_cross(tracks=[invalid_tracks]) diff --git a/pygmt/x2sys.py b/pygmt/x2sys.py index 4f55b23f3dc..22294c3e791 100644 --- a/pygmt/x2sys.py +++ b/pygmt/x2sys.py @@ -2,6 +2,8 @@ GMT supplementary X2SYS module for crossover analysis. """ import contextlib +import os +from pathlib import Path import pandas as pd @@ -14,10 +16,45 @@ dummy_context, fmt_docstring, kwargs_to_strings, + unique_name, use_alias, ) +@contextlib.contextmanager +def tempfile_from_dftrack(track, suffix): + """ + Saves pandas.DataFrame track table to a temporary tab-separated ASCII text + file with a unique name (to prevent clashes when running x2sys_cross), + adding a suffix extension to the end. + + Parameters + ---------- + track : pandas.DataFrame + A table holding track data with coordinate (x, y) or (lon, lat) values, + and (optionally) time (t). + suffix : str + File extension, e.g. xyz, tsv, etc. + + Yields + ------ + tmpfilename : str + A temporary tab-separated value file with a unique name holding the + track data. E.g. 'track-1a2b3c4.tsv'. + """ + try: + tmpfilename = f"track-{unique_name()[:7]}.{suffix}" + track.to_csv( + path_or_buf=tmpfilename, + sep="\t", + index=False, + date_format="%Y-%m-%dT%H:%M:%S.%fZ", + ) + yield tmpfilename + finally: + os.remove(tmpfilename) + + @fmt_docstring @use_alias( D="fmtfile", @@ -158,9 +195,10 @@ def x2sys_cross(tracks=None, outfile=None, **kwargs): Parameters ---------- - tracks : str or list + tracks : pandas.DataFrame or str or list A table or a list of tables with (x, y) or (lon, lat) values in the - first two columns. Supported formats are ASCII, native binary, or + first two columns. Track(s) can be provided as pandas DataFrame tables + or file names. Supported file formats are ASCII, native binary, or COARDS netCDF 1-D data. More columns may also be present. If the filenames are missing their file extension, we will append the @@ -263,8 +301,20 @@ def x2sys_cross(tracks=None, outfile=None, **kwargs): if kind == "file": file_contexts.append(dummy_context(track)) elif kind == "matrix": - raise NotImplementedError(f"{type(track)} inputs are not supported yet") - # file_contexts.append(lib.virtualfile_from_matrix(track.values)) + # find suffix (-E) of trackfiles used (e.g. xyz, csv, etc) from + # $X2SYS_HOME/TAGNAME/TAGNAME.tag file + lastline = ( + Path(os.environ["X2SYS_HOME"], kwargs["T"], f"{kwargs['T']}.tag") + .read_text() + .strip() + .split("\n")[-1] + ) # e.g. "-Dxyz -Etsv -I1/1" + for item in sorted(lastline.split()): # sort list alphabetically + if item.startswith(("-E", "-D")): # prefer -Etsv over -Dxyz + suffix = item[2:] # e.g. tsv (1st choice) or xyz (2nd choice) + + # Save pandas.DataFrame track data to temporary file + file_contexts.append(tempfile_from_dftrack(track=track, suffix=suffix)) else: raise GMTInvalidInput(f"Unrecognized data type: {type(track)}") @@ -287,8 +337,8 @@ def x2sys_cross(tracks=None, outfile=None, **kwargs): parse_dates=[2, 3], # Datetimes on 3rd and 4th column ) # Remove the "# " from "# x" in the first column - result = table.rename(columns={table.columns[0]: table.columns[0][2:]}) + table = table.rename(columns={table.columns[0]: table.columns[0][2:]}) elif outfile != tmpfile.name: # if outfile is set, output in outfile only - result = None + table = None - return result + return table