diff --git a/doc/api/index.rst b/doc/api/index.rst index 75dedbf737d..3aa176f1f2e 100644 --- a/doc/api/index.rst +++ b/doc/api/index.rst @@ -129,6 +129,7 @@ Operations on tabular data blockmedian blockmode filter1d + fitcircle nearneighbor project select diff --git a/pygmt/__init__.py b/pygmt/__init__.py index 67c9c583d5a..b8c87737138 100644 --- a/pygmt/__init__.py +++ b/pygmt/__init__.py @@ -34,6 +34,7 @@ config, dimfilter, filter1d, + fitcircle, grd2cpt, grd2xyz, grdclip, diff --git a/pygmt/helpers/caching.py b/pygmt/helpers/caching.py index e315a0eddcc..51c886f4fee 100644 --- a/pygmt/helpers/caching.py +++ b/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", diff --git a/pygmt/src/__init__.py b/pygmt/src/__init__.py index 7f77088edb9..4e500b5b7be 100644 --- a/pygmt/src/__init__.py +++ b/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 diff --git a/pygmt/src/fitcircle.py b/pygmt/src/fitcircle.py new file mode 100644 index 00000000000..3cdf8b29159 --- /dev/null +++ b/pygmt/src/fitcircle.py @@ -0,0 +1,142 @@ +""" +fitcircle - Find mean position and great or small circle fit to points on sphere. +""" + +from typing import Literal + +from pygmt._typing import PathLike, TableLike +from pygmt.alias import Alias, AliasSystem +from pygmt.clib import Session +from pygmt.helpers import build_arg_list, fmt_docstring +from pygmt.helpers.utils import is_given + + +@fmt_docstring +def fitcircle( + data: PathLike | TableLike | None = None, + x=None, + y=None, + norm: Literal["absolutes", "squares"] = "squares", + small_circle: bool | float = False, + verbose: Literal["quiet", "error", "warning", "timing", "info", "compat", "debug"] + | bool = False, + **kwargs, +) -> dict[str, tuple[float, float] | float]: + """ + Find mean position and great or small circle fit to points on sphere. + + This method takes (longitude, latitude) values and converts them 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. + + Setting ``norm`` to ``"absolutes"`` 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 ``"squares"`` 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. + + 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 by calling :func:`pygmt.fitcircle` + twice, once for each ``norm``. + + Takes a matrix, (x, y) pairs, or a file name as input. + + Must provide either ``data`` or ``x`` and ``y``. + + Full GMT docs at :gmt-docs:`fitcircle.html`. + + **Aliases:** + + .. hlist:: + :columns: 3 + + - L = norm + - S = small_circle + - V = verbose + + Parameters + ---------- + data + Pass in (longitude, latitude) values by providing a file name to an ASCII data + table, a 2-D $table_classes. + x/y : 1-D arrays + Arrays of x and y coordinates of the data points. + norm + Specify the desired norm, either ``"absolutes"`` or ``"squares"`` + [Default is ``"squares"``]. + small_circle + 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 set the desired fixed + latitude of the small circle [Default will determine the optimal latitude]. + $verbose + + Returns + ------- + ret + A dictionary with the following keys, each mapping to a + ``(longitude, latitude)`` tuple: + + - ``"flat_mean"``: the flat Earth mean position + - ``"mean"``: the mean position (Fisher or eigenvalue method, + depending on ``norm``) + - ``"north_pole"``: the north hemisphere great circle pole + - ``"south_pole"``: the south hemisphere great circle pole + + If ``small_circle`` is set, two more keys are added: + + - ``"small_circle_pole"``: the small circle pole + - ``"small_circle_distance"``: the colatitude/distance in degrees + from the small circle pole to the small circle (a ``float``, not a + tuple) + """ + aliasdict = AliasSystem( + L=Alias(norm, name="norm", mapping={"absolutes": 1, "squares": 2}), + S=Alias(small_circle, name="small_circle"), + ).add_common( + V=verbose, + ) + aliasdict.merge(kwargs) + + # "c" (small-circle pole and colatitude) is only valid with -S; GMT errors + # ("Cannot select c without setting -S") if "c" is requested without it. + aliasdict["F"] = "fmnsc" if is_given(small_circle) else "fmns" + + with Session() as lib: + with ( + lib.virtualfile_in( + check_kind="vector", data=data, x=x, y=y, mincols=2 + ) as vintbl, + lib.virtualfile_out(kind="dataset") as vouttbl, + ): + lib.call_module( + module="fitcircle", + args=build_arg_list(aliasdict, infile=vintbl, outfile=vouttbl), + ) + row = lib.virtualfile_to_dataset(vfname=vouttbl, output_type="numpy")[0] + values = [float(value) for value in row] + + solution: dict[str, tuple[float, float] | float] = { + "flat_mean": (values[0], values[1]), + "mean": (values[2], values[3]), + "north_pole": (values[4], values[5]), + "south_pole": (values[6], values[7]), + } + if is_given(small_circle): + solution["small_circle_pole"] = (values[8], values[9]) + solution["small_circle_distance"] = values[10] + return solution diff --git a/pygmt/tests/test_fitcircle.py b/pygmt/tests/test_fitcircle.py new file mode 100644 index 00000000000..b9ab3f64c81 --- /dev/null +++ b/pygmt/tests/test_fitcircle.py @@ -0,0 +1,75 @@ +""" +Test pygmt.fitcircle. +""" + +import numpy.testing as npt +import pandas as pd +import pytest +from pygmt import fitcircle +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_absolutes(data): + """ + Test fitcircle with norm="absolutes". + """ + result = fitcircle(data=data, norm="absolutes") + assert isinstance(result, dict) + assert set(result.keys()) == {"flat_mean", "mean", "north_pole", "south_pole"} + npt.assert_allclose(result["flat_mean"], (330.243649573, -18.3910128205)) + npt.assert_allclose(result["mean"], (330.16313328, -18.4067771888)) + npt.assert_allclose(result["north_pole"], (52.7434273422, 21.2085369093)) + npt.assert_allclose(result["south_pole"], (232.743427342, -21.2085369093)) + + +def test_fitcircle_squares(data): + """ + Test fitcircle with norm="squares", which is also the default. + """ + result = fitcircle(data=data, norm="squares") + assert isinstance(result, dict) + assert set(result.keys()) == {"flat_mean", "mean", "north_pole", "south_pole"} + npt.assert_allclose(result["flat_mean"], (330.243649573, -18.3910128205)) + npt.assert_allclose(result["mean"], (330.163207808, -18.4067882988)) + npt.assert_allclose(result["north_pole"], (52.7449849947, 21.2046833116)) + npt.assert_allclose(result["south_pole"], (232.744984995, -21.2046833116)) + assert fitcircle(data=data) == result # norm="squares" is the default + + +def test_fitcircle_small_circle(data): + """ + Test that fitcircle can fit a small circle instead of a great circle, and + that the returned dict includes the small-circle keys. + """ + result = fitcircle(data=data, norm="squares", small_circle=True) + assert isinstance(result, dict) + assert set(result.keys()) == { + "flat_mean", + "mean", + "north_pole", + "south_pole", + "small_circle_pole", + "small_circle_distance", + } + npt.assert_allclose(result["small_circle_distance"], 87.6072781238) + + +def test_fitcircle_input_xy(data): + """ + Run fitcircle by passing in x/y as input. + """ + result = fitcircle(x=data.longitude, y=data.latitude, norm="absolutes") + assert isinstance(result, dict) + npt.assert_allclose(result["flat_mean"], (330.243649573, -18.3910128205))