Wrap fitcircle

doc/api/index.rst

@@ -129,6 +129,7 @@ Operations on tabular data
     blockmedian
     blockmode
     filter1d
+    fitcircle
     nearneighbor
     project
     select

pygmt/__init__.py

@@ -35,6 +35,7 @@
     config,
     dimfilter,
     filter1d,
+    fitcircle,
     grd2cpt,
     grd2xyz,
     grdclip,

pygmt/helpers/caching.py

@@ -116,6 +116,7 @@ def cache_data() -> None:
         "@RidgeTest.prj",
         "@RidgeTest.shp",
         "@RidgeTest.shx",
+        "@sat_03.txt",
         "@SOEST_block4.png",
         "@Table_5_11.txt",
         "@Table_5_11_mean.xyz",

pygmt/src/__init__.py

@@ -9,6 +9,7 @@
 from pygmt.src.config import config
 from pygmt.src.dimfilter import dimfilter
 from pygmt.src.filter1d import filter1d
+from pygmt.src.fitcircle import fitcircle
 from pygmt.src.grd2cpt import grd2cpt
 from pygmt.src.grd2xyz import grd2xyz
 from pygmt.src.grdclip import grdclip

pygmt/src/fitcircle.py

@@ -0,0 +1,135 @@
+"""
+fitcircle - Find mean position and great [or small] circle fit to points on
+sphere.
+"""
+
+from typing import Literal
+
+import numpy as np
+import pandas as pd
+from pygmt._typing import PathLike, TableLike
+from pygmt.alias import AliasSystem
+from pygmt.clib import Session
+from pygmt.exceptions import GMTParameterError, GMTValueError
+from pygmt.helpers import (
+    build_arg_list,
+    fmt_docstring,
+    use_alias,
+    validate_output_table_type,
+)
+
+
+@fmt_docstring
+@use_alias(
+    L="norm",
+    S="small_circle",
+)
+def fitcircle(
+    data: PathLike | TableLike,
+    output_type: Literal["pandas", "numpy", "file"] = "pandas",
+    outfile: PathLike | None = None,
+    verbose: Literal["quiet", "error", "warning", "timing", "info", "compat", "debug"]
+    | bool = False,
+    **kwargs,
+) -> pd.DataFrame | np.ndarray | None:
+    r"""
+    Find mean position and great [or small] circle fit to points on sphere.
+
+    **fitcircle** reads (longitude, latitude) or (latitude, longitude) values from the
+    first two columns of the input data. These are converted to Cartesian
+    three-vectors on the unit sphere. Then two locations are found: the mean
+    of the input positions, and the pole to the great circle which best fits
+    the input positions. The user may choose one or both of two possible
+    solutions to this problem. When the data are closely grouped along a
+    great circle both solutions are similar. If the data have large
+    dispersion, the pole to the great circle will be less well determined
+    than the mean. Compare both solutions as a qualitative check.
+
+    Setting ``norm`` to **1** approximates the minimization of the sum of
+    absolute values of cosines of angular distances. This solution finds the
+    mean position as the Fisher average of the data, and the pole position
+    as the Fisher average of the cross-products between the mean and the
+    data. Averaging cross-products gives weight to points in proportion to
+    their distance from the mean, analogous to the "leverage" of distant
+    points in linear regression in the plane.
+
+    Setting ``norm`` to **2** approximates the minimization of the sum of
+    squares of cosines of angular distances. It creates a 3 by 3 matrix of
+    sums of squares of components of the data vectors. The eigenvectors of
+    this matrix give the mean and pole locations. This method may be more
+    subject to roundoff errors when there are thousands of data. The pole is
+    given by the eigenvector corresponding to the smallest eigenvalue; it is
+    the least-well represented factor in the data and is not easily
+    estimated by either method.
+
+    Full GMT docs at :gmt-docs:`fitcircle.html`.
+
+    $aliases
+       - V = verbose
+
+    Parameters
+    ----------
+    data
+        Pass in (longitude, latitude) or (latitude, longitude) values by
+        providing a file name to an ASCII data table, a 2-D
+        $table_classes.
+    $output_type
+    $outfile
+    norm : int or bool
+        Specify the desired *norm* as **1** or **2**\ , or use ``True`` or
+        **3** to see both solutions. Note that ``output_type="pandas"`` is
+        not supported when ``norm`` is ``True`` or **3**; use
+        ``output_type="numpy"`` or ``output_type="file"`` instead.
+    small_circle : bool or float
+        Attempt to fit a small circle instead of a great circle. The pole
+        will be constrained to lie on the great circle connecting the pole
+        of the best-fit great circle and the mean location of the data.
+        Optionally append the desired fixed latitude of the small circle
+        [Default will determine the optimal latitude].
+    $verbose
+
+    Returns
+    -------
+    ret
+        Return type depends on ``outfile`` and ``output_type``:
+
+        - ``None`` if ``outfile`` is set (output will be stored in the file set by
+          ``outfile``)
+        - :class:`pandas.DataFrame` or :class:`numpy.ndarray` if ``outfile`` is not set
+          (depends on ``output_type``)
+    """
+    if kwargs.get("L") is None:
+        raise GMTParameterError(required="norm")
+
+    output_type = validate_output_table_type(output_type, outfile=outfile)
+    norm = kwargs.get("L")
+    if output_type == "pandas" and (norm is True or norm == 3):
+        raise GMTValueError(
+            norm,
+            description="value for parameter 'norm'",
+            reason=(
+                "Pandas output is not supported when 'norm' is set to True or 3 "
+                "since both L1 and L2 solutions are stacked in the same rows. "
+                "Use output_type='numpy' or output_type='file' instead."
+            ),
+        )
+
+    aliasdict = AliasSystem().add_common(
+        V=verbose,
+    )
+    aliasdict.merge(kwargs)
+
+    with Session() as lib:
+        with (
+            lib.virtualfile_in(check_kind="vector", data=data) as vintbl,
+            lib.virtualfile_out(kind="dataset", fname=outfile) as vouttbl,
+        ):
+            lib.call_module(
+                module="fitcircle",
+                args=build_arg_list(aliasdict, infile=vintbl, outfile=vouttbl),
+            )
+        return lib.virtualfile_to_dataset(
+            vfname=vouttbl,
+            output_type=output_type,
+            column_names=["longitude", "latitude", "method"],
+        )

pygmt/tests/test_fitcircle.py

@@ -0,0 +1,130 @@
+"""
+Test pygmt.fitcircle.
+"""
+
+from pathlib import Path
+
+import numpy as np
+import numpy.testing as npt
+import pandas as pd
+import pytest
+from pygmt import fitcircle
+from pygmt.exceptions import GMTParameterError, GMTValueError
+from pygmt.helpers import GMTTempFile
+from pygmt.src import which
+
+
+@pytest.fixture(scope="module", name="data")
+def fixture_data():
+    """
+    Load the sample data from the @sat_03 remote file.
+    """
+    fname = which("@sat_03.txt", download="c")
+    return pd.read_csv(
+        fname, header=None, skiprows=1, sep="\t", names=["longitude", "latitude", "z"]
+    )
+
+
+@pytest.mark.benchmark
+def test_fitcircle_no_outfile(data):
+    """
+    Test fitcircle with no set outfile.
+    """
+    result = fitcircle(data=data, norm=2)
+    assert isinstance(result, pd.DataFrame)
+    assert result.shape == (4, 3)
+    # Test longitude results
+    npt.assert_allclose(result.longitude.min(), 52.7449849947)
+    npt.assert_allclose(result.longitude.max(), 330.243649573)
+    # Test latitude results
+    npt.assert_allclose(result.latitude.min(), -21.2046833116)
+    npt.assert_allclose(result.latitude.max(), 21.2046833116)
+
+
+def test_fitcircle_file_output(data):
+    """
+    Test that fitcircle returns a file output when it is specified.
+    """
+    with GMTTempFile(suffix=".txt") as tmpfile:
+        result = fitcircle(
+            data=data, norm=True, outfile=tmpfile.name, output_type="file"
+        )
+        assert result is None  # return value is None
+        assert Path(tmpfile.name).stat().st_size > 0  # check that outfile exists
+
+
+def test_fitcircle_invalid_format(data):
+    """
+    Test that fitcircle fails with an incorrect format for output_type.
+    """
+    with pytest.raises(GMTValueError):
+        fitcircle(data=data, norm=True, output_type="a")
+
+
+def test_fitcircle_no_norm(data):
+    """
+    Test that fitcircle fails when the required "norm" parameter is missing.
+    """
+    with pytest.raises(GMTParameterError):
+        fitcircle(data=data)
+
+
+def test_fitcircle_no_outfile_specified(data):
+    """
+    Test that fitcircle fails when output_type is set to "file" but no outfile
+    is specified.
+    """
+    with pytest.raises(GMTParameterError):
+        fitcircle(data=data, norm=True, output_type="file")
+
+
+def test_fitcircle_outfile_incorrect_output_type(data):
+    """
+    Test that fitcircle raises a warning when an outfile filename is set but the
+    output_type is not set to "file".
+    """
+    with GMTTempFile(suffix=".txt") as tmpfile:
+        with pytest.warns(RuntimeWarning) as record:
+            result = fitcircle(
+                data=data, norm=True, outfile=tmpfile.name, output_type="numpy"
+            )
+        assert len(record) == 1  # check that only one warning was raised
+        assert result is None  # return value is None
+        assert Path(tmpfile.name).stat().st_size > 0  # check that outfile exists
+
+
+def test_fitcircle_format(data):
+    """
+    Test that correct formats are returned.
+    """
+    circle_default = fitcircle(data=data, norm=2)
+    assert isinstance(circle_default, pd.DataFrame)
+    circle_array = fitcircle(data=data, norm=2, output_type="numpy")
+    assert isinstance(circle_array, np.ndarray)
+    circle_df = fitcircle(data=data, norm=2, output_type="pandas")
+    assert isinstance(circle_df, pd.DataFrame)
+
+
+@pytest.mark.parametrize("norm", [True, 3])
+def test_fitcircle_pandas_unsupported_for_both_norms(data, norm):
+    """
+    Test that fitcircle raises an exception when output_type is "pandas" (the
+    default) and norm is True or 3, since the L1 and L2 solutions are stacked
+    in the same rows and can't be represented as a single pandas.DataFrame.
+    """
+    with pytest.raises(GMTValueError):
+        fitcircle(data=data, norm=norm)
+    with pytest.raises(GMTValueError):
+        fitcircle(data=data, norm=norm, output_type="pandas")
+    result = fitcircle(data=data, norm=norm, output_type="numpy")
+    assert isinstance(result, np.ndarray)
+
+
+def test_fitcircle_small_circle(data):
+    """
+    Test that fitcircle can fit a small circle instead of a great circle.
+    """
+    result = fitcircle(data=data, norm=2, small_circle=True)
+    assert isinstance(result, pd.DataFrame)
+    assert result.shape == (5, 3)
+    assert "Small Circle Pole" in result.method.iloc[-1]