We are currently solving the following issue within our repository. Here is the issue text:
--- BEGIN ISSUE ---
ENH: New geoaccessor to calculate two coordinates distance on earth
<!--
Thanks for contributing a pull request!

Please follow these standard acronyms to start the commit message:

- ENH: enhancement
- BUG: bug fix
- DOC: documentation
- TYP: type annotations
- TST: addition or modification of tests
- MAINT: maintenance commit (refactoring, typos, etc.)
- BLD: change related to building
- REL: related to releasing
- API: an (incompatible) API change
- DEP: deprecate something, or remove a deprecated object
- DEV: development tool or utility
- REV: revert an earlier commit
- PERF: performance improvement
- BOT: always commit via a bot
- CI: related to CI or CD
- CLN: Code cleanup
-->

- [ ] closes #xxxx
- [ ] whatsnew entry


--- END ISSUE ---

Below are some code segments, each from a relevant file. One or more of these files may contain bugs.
--- BEGIN FILES ---
Path: `dtoolkit/geoaccessor/geoseries/geodistance.py`
Content:
```
1 from __future__ import annotations
2 
3 from warnings import warn
4 
5 import geopandas as gpd
6 import numpy as np
7 import pandas as pd
8 from shapely.geometry.base import BaseGeometry
9 
10 from dtoolkit.geoaccessor.register import register_geoseries_method
11 from dtoolkit.util._exception import find_stack_level
12 
13 
14 @register_geoseries_method
15 def geodistance(
16     s: gpd.GeoSeries,
17     /,
18     other: BaseGeometry | gpd.GeoSeries | gpd.GeoDataFrame,
19     align: bool = True,
20     radius: float = 6371008.7714150598,
21 ) -> pd.Series:
22     """
23     Returns a ``Series`` containing the great-circle distance to aligned other.
24 
25     The algorithm uses the Vincenty formula which is more accurate than the Haversine
26     formula.
27 
28     .. math::
29 
30         D(x, y) = \\arctan[
31             \\frac{
32                 \\sqrt{
33                     (
34                         \\cos(y_1) \\sin(y_2)
35                         - \\sin(y_1) \\cos(y_2) \\cos(x_2 - x_1)
36                     )^2
37                     + (\\cos(y_2) \\sin(x_2 - x_1))^2
38                 }
39             }{
40                 \\sin(y_1) \\sin(y_2)
41                 + \\cos(y_1) \\cos(y_2) \\cos(x_2 - x_1)
42             }
43         ]
44 
45     Parameters
46     ----------
47     other : BaseGeometry, GeoSeries, or GeoDataFrame
48 
49     align : bool, default True
50         If True, automatically aligns GeoSeries based on their indices. If False,
51         the order of elements is preserved.
52 
53     radius : float, default 6371008.7714150598
54         Great-circle distance uses a spherical model of the earth, using the mean earth
55         radius as defined by the International Union of Geodesy and Geophysics,
56         (2\\ *a* + *b*)/3 = 6371008.7714150598 meters for WGS-84.
57 
58     Returns
59     -------
60     Series
61         The values are the great-circle distances and its unit is meters.
62 
63     Raises
64     ------
65     ValueError
66         If the CRS is not ``ESGP:4326``.
67 
68     TypeError
69         If the other is not a ``BaseGeometry``, ``GeoSeries``, or ``GeoDataFrame``.
70 
71     See Also
72     --------
73     geopandas.GeoSeries.distance
74     dtoolkit.geoaccessor.geoseries.geodistance
75     dtoolkit.geoaccessor.geoseries.geodistance_matrix
76     dtoolkit.geoaccessor.geodataframe.geodistance
77     dtoolkit.geoaccessor.geodataframe.geodistance_matrix
78 
79     Notes
80     -----
81     - Currently, only supports Point geometry.
82     - The geodesic distance is the shortest distance on the surface of an ellipsoidal
83       model of the earth. Resulting in an error of up to about 0.5%.
84 
85     Examples
86     --------
87     >>> import dtoolkit.geoaccessor
88     >>> import pandas as pd
89     >>> import geopandas as gpd
90     >>> from shapely.geometry import Point
91     >>> df = (
92     ...     pd.DataFrame({"x": [122, 100], "y":[55, 1]})
93     ...     .from_xy("x", "y", crs=4326)
94     ... )
95     >>> df
96          x   y                    geometry
97     0  122  55  POINT (122.00000 55.00000)
98     1  100   1   POINT (100.00000 1.00000)
99     >>> df.geodistance(Point(120, 30)) / 1e6
100     0    2.784435
101     1    3.855604
102     dtype: float64
103 
104     Calculate the great-circle distance of corresponding points.
105 
106     >>> s = gpd.GeoSeries([Point(120, 30), Point(120, 50)], index=[1, 2], crs=4326)
107     >>> s
108     1    POINT (120.00000 30.00000)
109     2    POINT (120.00000 50.00000)
110     dtype: geometry
111     >>> df.geodistance(s) / 1e6
112     0         NaN
113     1    3.855604
114     2         NaN
115     dtype: float64
116     >>> df.geodistance(s, align=False) / 1e6
117     0    2.784435
118     1    5.768885
119     dtype: float64
120     """
121 
122     if s.crs != 4326:
123         raise ValueError(f"Only support 'EPSG:4326' CRS, but got {s.crs!r}.")
124     if not isinstance(other, (BaseGeometry, gpd.base.GeoPandasBase)):
125         raise TypeError(f"Unknown type: {type(other).__name__!r}.")
126 
127     if isinstance(other, gpd.base.GeoPandasBase):
128         if other.crs != 4326:
129             raise ValueError(f"Only support 'EPSG:4326' CRS, but got {other.crs!r}.")
130 
131         s = s.geometry
132         if align and not s.index.equals(other.index):
133             warn("The indices are different.", stacklevel=find_stack_level())
134             s, other = s.align(other)
135         else:
136             other = other.geometry
137 
138     return pd.Series(
139         distance(
140             s.geometry.x.to_numpy(),
141             s.geometry.y.to_numpy(),
142             other.x if isinstance(other, BaseGeometry) else other.x.to_numpy(),
143             other.y if isinstance(other, BaseGeometry) else other.y.to_numpy(),
144             radius=radius,
145         ),
146         index=s.index,
147     )
148 
149 
150 # based on https://github.com/geopy/geopy geopy/distance.py::great_circle.measure
151 def distance(
152     lng1: np.ndarray | float,
153     lat1: np.ndarray | float,
154     lng2: np.ndarray | float,
155     lat2: np.ndarray | float,
156     radius: float,
157 ) -> np.ndarray:
158     lng1, lat1, lng2, lat2 = map(np.radians, (lng1, lat1, lng2, lat2))
159     sin_lat1, cos_lat1 = np.sin(lat1), np.cos(lat1)
160     sin_lat2, cos_lat2 = np.sin(lat2), np.cos(lat2)
161 
162     delta_lng = lng2 - lng1
163     cos_delta_lng, sin_delta_lng = np.cos(delta_lng), np.sin(delta_lng)
164 
165     return radius * np.arctan2(
166         np.sqrt(
167             (cos_lat1 * sin_lat2 - sin_lat1 * cos_lat2 * cos_delta_lng) ** 2
168             + (cos_lat2 * sin_delta_lng) ** 2,
169         ),
170         sin_lat1 * sin_lat2 + cos_lat1 * cos_lat2 * cos_delta_lng,
171     )
172
```

Path: `dtoolkit/geoaccessor/geoseries/geodistance_matrix.py`
Content:
```
1 from __future__ import annotations
2 
3 import geopandas as gpd
4 import numpy as np
5 import pandas as pd
6 
7 from dtoolkit.geoaccessor.register import register_geoseries_method
8 
9 
10 @register_geoseries_method
11 def geodistance_matrix(
12     s: gpd.GeoSeries,
13     /,
14     other: gpd.GeoSeries | gpd.GeoDataFrame | None = None,
15     radius: float = 6371008.7714150598,
16 ) -> pd.DataFrame:
17     """
18     Returns a ``DataFrame`` containing the great-circle distances matrix between in
19     ``s`` and ``other`` via haversine formula.
20 
21     .. math::
22 
23         D(x, y) = 2 \\arcsin [
24             \\sqrt{
25                 \\sin^2 ((x_1 - y_1) / 2)
26                 + \\cos(x_1) \\cos(y_1) \\sin^2 ((x_2 - y_2) / 2)
27             }
28         ]
29 
30     Parameters
31     ----------
32     other : GeoSeries, or GeoDataFrame, default None
33         If None, uses ``other=s``.
34 
35     radius : float, default 6371008.7714150598
36         Great-circle distance uses a spherical model of the earth, using the mean earth
37         radius as defined by the International Union of Geodesy and Geophysics,
38         (2\\ *a* + *b*)/3 = 6371008.7714150598 meters for WGS-84.
39 
40     Returns
41     -------
42     DataFrame
43         - The index and columns are the same as the index of ``s`` and ``other``.
44         - The values are the great-circle distances and its unit is meters.
45 
46     Raises
47     ------
48     ModuleNotFoundError
49         If don't have module named 'sklearn'.
50 
51     ValueError
52         If the CRS is not ``ESGP:4326``.
53 
54     See Also
55     --------
56     sklearn.metrics.pairwise.haversine_distances
57     dtoolkit.geoaccessor.geoseries.geodistance
58     dtoolkit.geoaccessor.geoseries.geodistance_matrix
59     dtoolkit.geoaccessor.geodataframe.geodistance
60     dtoolkit.geoaccessor.geodataframe.geodistance_matrix
61 
62     Notes
63     -----
64     - Currently, only supports Point geometry.
65     - The great-circle distance is the angular distance between two points on the
66       surface of a sphere. As the Earth is nearly spherical, the haversine formula
67       provides a good approximation of the distance between two points of the Earth
68       surface, with a less than 1% error on average.
69 
70     Examples
71     --------
72     >>> import dtoolkit.geoaccessor
73     >>> df = pd.DataFrame(
74     ...     {
75     ...         "x": [120, 122, 100],
76     ...         "y":[30, 55, 1],
77     ...     },
78     ... ).from_xy("x", "y", crs=4326)
79     >>> df
80          x   y                    geometry
81     0  120  30  POINT (120.00000 30.00000)
82     1  122  55  POINT (122.00000 55.00000)
83     2  100   1   POINT (100.00000 1.00000)
84     >>> other = pd.DataFrame(
85     ...     {
86     ...         "x": [120, 110],
87     ...         "y":[30, 40],
88     ...     },
89     ... ).from_xy("x", "y", crs=4326)
90     >>> other
91          x   y                    geometry
92     0  120  30  POINT (120.00000 30.00000)
93     1  110  40  POINT (110.00000 40.00000)
94     >>> df.geodistance_matrix(other)
95                   0             1
96     0  0.000000e+00  1.203540e+06
97     1  1.439971e+06  1.511958e+06
98     2  2.418544e+06  1.522752e+06
99     """
100     from sklearn.metrics.pairwise import haversine_distances
101 
102     if s.crs != 4326:
103         raise ValueError(f"Only support 'EPSG:4326' CRS, but got {s.crs!r}.")
104     if other is not None and other.crs != 4326:
105         raise ValueError(f"Only support 'EPSG:4326' CRS, but got {other.crs!r}.")
106 
107     X = np.radians(np.stack((s.x, s.y), axis=1))
108     Y = np.radians(np.stack((other.x, other.y), axis=1)) if other is not None else other
109     return pd.DataFrame(
110         radius * haversine_distances(X, Y),
111         index=s.index,
112         columns=other.index,
113     )
114
```
--- END FILES ---

Please first localize the bug based on the issue statement, and then generate a patch according to the `git diff` format fenced by three backticks.

Here is an example:

```diff
diff --git a/examples/server_async.py b/examples/server_async.py
--- a/examples/server_async.py
+++ b/examples/server_async.py
@@ -313,4 +313,4 @@
 
 
 if __name__ == "__main__":
-    asyncio.run(run_async_server("."), debug=True)
+    asyncio.run(run_async_server(), debug=True)
diff --git a/examples/server_sync.py b/examples/server_sync.py
--- a/examples/server_sync.py
+++ b/examples/server_sync.py
@@ -313,5 +313,5 @@
 
 
 if __name__ == "__main__":
-    server = run_sync_server(".")
+    server = run_sync_server()
     server.shutdown()

```

diff --git a/dtoolkit/geoaccessor/geoseries/geodistance.py b/dtoolkit/geoaccessor/geoseries/geodistance.py
--- a/dtoolkit/geoaccessor/geoseries/geodistance.py
+++ b/dtoolkit/geoaccessor/geoseries/geodistance.py
@@ -20,7 +20,9 @@
     radius: float = 6371008.7714150598,
 ) -> pd.Series:
     """
-    Returns a ``Series`` containing the great-circle distance to aligned other.
+    Returns a ``Series`` containing the `great-circle`__ distance to aligned other.
+
+    __ https://en.wikipedia.org/wiki/Great-circle_distance
 
     The algorithm uses the Vincenty formula which is more accurate than the Haversine
     formula.
diff --git a/dtoolkit/geoaccessor/geoseries/geodistance_matrix.py b/dtoolkit/geoaccessor/geoseries/geodistance_matrix.py
--- a/dtoolkit/geoaccessor/geoseries/geodistance_matrix.py
+++ b/dtoolkit/geoaccessor/geoseries/geodistance_matrix.py
@@ -15,9 +15,11 @@
     radius: float = 6371008.7714150598,
 ) -> pd.DataFrame:
     """
-    Returns a ``DataFrame`` containing the great-circle distances matrix between in
+    Returns a ``DataFrame`` containing the `great-circle`__ distances matrix between in
     ``s`` and ``other`` via haversine formula.
 
+    __ https://en.wikipedia.org/wiki/Great-circle_distance
+
     .. math::
 
         D(x, y) = 2 \\arcsin [

{"golden_diff": "diff --git a/dtoolkit/geoaccessor/geoseries/geodistance.py b/dtoolkit/geoaccessor/geoseries/geodistance.py\n--- a/dtoolkit/geoaccessor/geoseries/geodistance.py\n+++ b/dtoolkit/geoaccessor/geoseries/geodistance.py\n@@ -20,7 +20,9 @@\n     radius: float = 6371008.7714150598,\n ) -> pd.Series:\n     \"\"\"\n-    Returns a ``Series`` containing the great-circle distance to aligned other.\n+    Returns a ``Series`` containing the `great-circle`__ distance to aligned other.\n+\n+    __ https://en.wikipedia.org/wiki/Great-circle_distance\n \n     The algorithm uses the Vincenty formula which is more accurate than the Haversine\n     formula.\ndiff --git a/dtoolkit/geoaccessor/geoseries/geodistance_matrix.py b/dtoolkit/geoaccessor/geoseries/geodistance_matrix.py\n--- a/dtoolkit/geoaccessor/geoseries/geodistance_matrix.py\n+++ b/dtoolkit/geoaccessor/geoseries/geodistance_matrix.py\n@@ -15,9 +15,11 @@\n     radius: float = 6371008.7714150598,\n ) -> pd.DataFrame:\n     \"\"\"\n-    Returns a ``DataFrame`` containing the great-circle distances matrix between in\n+    Returns a ``DataFrame`` containing the `great-circle`__ distances matrix between in\n     ``s`` and ``other`` via haversine formula.\n \n+    __ https://en.wikipedia.org/wiki/Great-circle_distance\n+\n     .. math::\n \n         D(x, y) = 2 \\\\arcsin [\n", "issue": "ENH: New geoaccessor to calculate two coordinates distance on earth\n<!--\r\nThanks for contributing a pull request!\r\n\r\nPlease follow these standard acronyms to start the commit message:\r\n\r\n- ENH: enhancement\r\n- BUG: bug fix\r\n- DOC: documentation\r\n- TYP: type annotations\r\n- TST: addition or modification of tests\r\n- MAINT: maintenance commit (refactoring, typos, etc.)\r\n- BLD: change related to building\r\n- REL: related to releasing\r\n- API: an (incompatible) API change\r\n- DEP: deprecate something, or remove a deprecated object\r\n- DEV: development tool or utility\r\n- REV: revert an earlier commit\r\n- PERF: performance improvement\r\n- BOT: always commit via a bot\r\n- CI: related to CI or CD\r\n- CLN: Code cleanup\r\n-->\r\n\r\n- [ ] closes #xxxx\r\n- [ ] whatsnew entry\r\n\n", "before_files": [{"content": "from __future__ import annotations\n\nfrom warnings import warn\n\nimport geopandas as gpd\nimport numpy as np\nimport pandas as pd\nfrom shapely.geometry.base import BaseGeometry\n\nfrom dtoolkit.geoaccessor.register import register_geoseries_method\nfrom dtoolkit.util._exception import find_stack_level\n\n\n@register_geoseries_method\ndef geodistance(\n    s: gpd.GeoSeries,\n    /,\n    other: BaseGeometry | gpd.GeoSeries | gpd.GeoDataFrame,\n    align: bool = True,\n    radius: float = 6371008.7714150598,\n) -> pd.Series:\n    \"\"\"\n    Returns a ``Series`` containing the great-circle distance to aligned other.\n\n    The algorithm uses the Vincenty formula which is more accurate than the Haversine\n    formula.\n\n    .. math::\n\n        D(x, y) = \\\\arctan[\n            \\\\frac{\n                \\\\sqrt{\n                    (\n                        \\\\cos(y_1) \\\\sin(y_2)\n                        - \\\\sin(y_1) \\\\cos(y_2) \\\\cos(x_2 - x_1)\n                    )^2\n                    + (\\\\cos(y_2) \\\\sin(x_2 - x_1))^2\n                }\n            }{\n                \\\\sin(y_1) \\\\sin(y_2)\n                + \\\\cos(y_1) \\\\cos(y_2) \\\\cos(x_2 - x_1)\n            }\n        ]\n\n    Parameters\n    ----------\n    other : BaseGeometry, GeoSeries, or GeoDataFrame\n\n    align : bool, default True\n        If True, automatically aligns GeoSeries based on their indices. If False,\n        the order of elements is preserved.\n\n    radius : float, default 6371008.7714150598\n        Great-circle distance uses a spherical model of the earth, using the mean earth\n        radius as defined by the International Union of Geodesy and Geophysics,\n        (2\\\\ *a* + *b*)/3 = 6371008.7714150598 meters for WGS-84.\n\n    Returns\n    -------\n    Series\n        The values are the great-circle distances and its unit is meters.\n\n    Raises\n    ------\n    ValueError\n        If the CRS is not ``ESGP:4326``.\n\n    TypeError\n        If the other is not a ``BaseGeometry``, ``GeoSeries``, or ``GeoDataFrame``.\n\n    See Also\n    --------\n    geopandas.GeoSeries.distance\n    dtoolkit.geoaccessor.geoseries.geodistance\n    dtoolkit.geoaccessor.geoseries.geodistance_matrix\n    dtoolkit.geoaccessor.geodataframe.geodistance\n    dtoolkit.geoaccessor.geodataframe.geodistance_matrix\n\n    Notes\n    -----\n    - Currently, only supports Point geometry.\n    - The geodesic distance is the shortest distance on the surface of an ellipsoidal\n      model of the earth. Resulting in an error of up to about 0.5%.\n\n    Examples\n    --------\n    >>> import dtoolkit.geoaccessor\n    >>> import pandas as pd\n    >>> import geopandas as gpd\n    >>> from shapely.geometry import Point\n    >>> df = (\n    ...     pd.DataFrame({\"x\": [122, 100], \"y\":[55, 1]})\n    ...     .from_xy(\"x\", \"y\", crs=4326)\n    ... )\n    >>> df\n         x   y                    geometry\n    0  122  55  POINT (122.00000 55.00000)\n    1  100   1   POINT (100.00000 1.00000)\n    >>> df.geodistance(Point(120, 30)) / 1e6\n    0    2.784435\n    1    3.855604\n    dtype: float64\n\n    Calculate the great-circle distance of corresponding points.\n\n    >>> s = gpd.GeoSeries([Point(120, 30), Point(120, 50)], index=[1, 2], crs=4326)\n    >>> s\n    1    POINT (120.00000 30.00000)\n    2    POINT (120.00000 50.00000)\n    dtype: geometry\n    >>> df.geodistance(s) / 1e6\n    0         NaN\n    1    3.855604\n    2         NaN\n    dtype: float64\n    >>> df.geodistance(s, align=False) / 1e6\n    0    2.784435\n    1    5.768885\n    dtype: float64\n    \"\"\"\n\n    if s.crs != 4326:\n        raise ValueError(f\"Only support 'EPSG:4326' CRS, but got {s.crs!r}.\")\n    if not isinstance(other, (BaseGeometry, gpd.base.GeoPandasBase)):\n        raise TypeError(f\"Unknown type: {type(other).__name__!r}.\")\n\n    if isinstance(other, gpd.base.GeoPandasBase):\n        if other.crs != 4326:\n            raise ValueError(f\"Only support 'EPSG:4326' CRS, but got {other.crs!r}.\")\n\n        s = s.geometry\n        if align and not s.index.equals(other.index):\n            warn(\"The indices are different.\", stacklevel=find_stack_level())\n            s, other = s.align(other)\n        else:\n            other = other.geometry\n\n    return pd.Series(\n        distance(\n            s.geometry.x.to_numpy(),\n            s.geometry.y.to_numpy(),\n            other.x if isinstance(other, BaseGeometry) else other.x.to_numpy(),\n            other.y if isinstance(other, BaseGeometry) else other.y.to_numpy(),\n            radius=radius,\n        ),\n        index=s.index,\n    )\n\n\n# based on https://github.com/geopy/geopy geopy/distance.py::great_circle.measure\ndef distance(\n    lng1: np.ndarray | float,\n    lat1: np.ndarray | float,\n    lng2: np.ndarray | float,\n    lat2: np.ndarray | float,\n    radius: float,\n) -> np.ndarray:\n    lng1, lat1, lng2, lat2 = map(np.radians, (lng1, lat1, lng2, lat2))\n    sin_lat1, cos_lat1 = np.sin(lat1), np.cos(lat1)\n    sin_lat2, cos_lat2 = np.sin(lat2), np.cos(lat2)\n\n    delta_lng = lng2 - lng1\n    cos_delta_lng, sin_delta_lng = np.cos(delta_lng), np.sin(delta_lng)\n\n    return radius * np.arctan2(\n        np.sqrt(\n            (cos_lat1 * sin_lat2 - sin_lat1 * cos_lat2 * cos_delta_lng) ** 2\n            + (cos_lat2 * sin_delta_lng) ** 2,\n        ),\n        sin_lat1 * sin_lat2 + cos_lat1 * cos_lat2 * cos_delta_lng,\n    )\n", "path": "dtoolkit/geoaccessor/geoseries/geodistance.py"}, {"content": "from __future__ import annotations\n\nimport geopandas as gpd\nimport numpy as np\nimport pandas as pd\n\nfrom dtoolkit.geoaccessor.register import register_geoseries_method\n\n\n@register_geoseries_method\ndef geodistance_matrix(\n    s: gpd.GeoSeries,\n    /,\n    other: gpd.GeoSeries | gpd.GeoDataFrame | None = None,\n    radius: float = 6371008.7714150598,\n) -> pd.DataFrame:\n    \"\"\"\n    Returns a ``DataFrame`` containing the great-circle distances matrix between in\n    ``s`` and ``other`` via haversine formula.\n\n    .. math::\n\n        D(x, y) = 2 \\\\arcsin [\n            \\\\sqrt{\n                \\\\sin^2 ((x_1 - y_1) / 2)\n                + \\\\cos(x_1) \\\\cos(y_1) \\\\sin^2 ((x_2 - y_2) / 2)\n            }\n        ]\n\n    Parameters\n    ----------\n    other : GeoSeries, or GeoDataFrame, default None\n        If None, uses ``other=s``.\n\n    radius : float, default 6371008.7714150598\n        Great-circle distance uses a spherical model of the earth, using the mean earth\n        radius as defined by the International Union of Geodesy and Geophysics,\n        (2\\\\ *a* + *b*)/3 = 6371008.7714150598 meters for WGS-84.\n\n    Returns\n    -------\n    DataFrame\n        - The index and columns are the same as the index of ``s`` and ``other``.\n        - The values are the great-circle distances and its unit is meters.\n\n    Raises\n    ------\n    ModuleNotFoundError\n        If don't have module named 'sklearn'.\n\n    ValueError\n        If the CRS is not ``ESGP:4326``.\n\n    See Also\n    --------\n    sklearn.metrics.pairwise.haversine_distances\n    dtoolkit.geoaccessor.geoseries.geodistance\n    dtoolkit.geoaccessor.geoseries.geodistance_matrix\n    dtoolkit.geoaccessor.geodataframe.geodistance\n    dtoolkit.geoaccessor.geodataframe.geodistance_matrix\n\n    Notes\n    -----\n    - Currently, only supports Point geometry.\n    - The great-circle distance is the angular distance between two points on the\n      surface of a sphere. As the Earth is nearly spherical, the haversine formula\n      provides a good approximation of the distance between two points of the Earth\n      surface, with a less than 1% error on average.\n\n    Examples\n    --------\n    >>> import dtoolkit.geoaccessor\n    >>> df = pd.DataFrame(\n    ...     {\n    ...         \"x\": [120, 122, 100],\n    ...         \"y\":[30, 55, 1],\n    ...     },\n    ... ).from_xy(\"x\", \"y\", crs=4326)\n    >>> df\n         x   y                    geometry\n    0  120  30  POINT (120.00000 30.00000)\n    1  122  55  POINT (122.00000 55.00000)\n    2  100   1   POINT (100.00000 1.00000)\n    >>> other = pd.DataFrame(\n    ...     {\n    ...         \"x\": [120, 110],\n    ...         \"y\":[30, 40],\n    ...     },\n    ... ).from_xy(\"x\", \"y\", crs=4326)\n    >>> other\n         x   y                    geometry\n    0  120  30  POINT (120.00000 30.00000)\n    1  110  40  POINT (110.00000 40.00000)\n    >>> df.geodistance_matrix(other)\n                  0             1\n    0  0.000000e+00  1.203540e+06\n    1  1.439971e+06  1.511958e+06\n    2  2.418544e+06  1.522752e+06\n    \"\"\"\n    from sklearn.metrics.pairwise import haversine_distances\n\n    if s.crs != 4326:\n        raise ValueError(f\"Only support 'EPSG:4326' CRS, but got {s.crs!r}.\")\n    if other is not None and other.crs != 4326:\n        raise ValueError(f\"Only support 'EPSG:4326' CRS, but got {other.crs!r}.\")\n\n    X = np.radians(np.stack((s.x, s.y), axis=1))\n    Y = np.radians(np.stack((other.x, other.y), axis=1)) if other is not None else other\n    return pd.DataFrame(\n        radius * haversine_distances(X, Y),\n        index=s.index,\n        columns=other.index,\n    )\n", "path": "dtoolkit/geoaccessor/geoseries/geodistance_matrix.py"}], "after_files": [{"content": "from __future__ import annotations\n\nfrom warnings import warn\n\nimport geopandas as gpd\nimport numpy as np\nimport pandas as pd\nfrom shapely.geometry.base import BaseGeometry\n\nfrom dtoolkit.geoaccessor.register import register_geoseries_method\nfrom dtoolkit.util._exception import find_stack_level\n\n\n@register_geoseries_method\ndef geodistance(\n    s: gpd.GeoSeries,\n    /,\n    other: BaseGeometry | gpd.GeoSeries | gpd.GeoDataFrame,\n    align: bool = True,\n    radius: float = 6371008.7714150598,\n) -> pd.Series:\n    \"\"\"\n    Returns a ``Series`` containing the `great-circle`__ distance to aligned other.\n\n    __ https://en.wikipedia.org/wiki/Great-circle_distance\n\n    The algorithm uses the Vincenty formula which is more accurate than the Haversine\n    formula.\n\n    .. math::\n\n        D(x, y) = \\\\arctan[\n            \\\\frac{\n                \\\\sqrt{\n                    (\n                        \\\\cos(y_1) \\\\sin(y_2)\n                        - \\\\sin(y_1) \\\\cos(y_2) \\\\cos(x_2 - x_1)\n                    )^2\n                    + (\\\\cos(y_2) \\\\sin(x_2 - x_1))^2\n                }\n            }{\n                \\\\sin(y_1) \\\\sin(y_2)\n                + \\\\cos(y_1) \\\\cos(y_2) \\\\cos(x_2 - x_1)\n            }\n        ]\n\n    Parameters\n    ----------\n    other : BaseGeometry, GeoSeries, or GeoDataFrame\n\n    align : bool, default True\n        If True, automatically aligns GeoSeries based on their indices. If False,\n        the order of elements is preserved.\n\n    radius : float, default 6371008.7714150598\n        Great-circle distance uses a spherical model of the earth, using the mean earth\n        radius as defined by the International Union of Geodesy and Geophysics,\n        (2\\\\ *a* + *b*)/3 = 6371008.7714150598 meters for WGS-84.\n\n    Returns\n    -------\n    Series\n        The values are the great-circle distances and its unit is meters.\n\n    Raises\n    ------\n    ValueError\n        If the CRS is not ``ESGP:4326``.\n\n    TypeError\n        If the other is not a ``BaseGeometry``, ``GeoSeries``, or ``GeoDataFrame``.\n\n    See Also\n    --------\n    geopandas.GeoSeries.distance\n    dtoolkit.geoaccessor.geoseries.geodistance\n    dtoolkit.geoaccessor.geoseries.geodistance_matrix\n    dtoolkit.geoaccessor.geodataframe.geodistance\n    dtoolkit.geoaccessor.geodataframe.geodistance_matrix\n\n    Notes\n    -----\n    - Currently, only supports Point geometry.\n    - The geodesic distance is the shortest distance on the surface of an ellipsoidal\n      model of the earth. Resulting in an error of up to about 0.5%.\n\n    Examples\n    --------\n    >>> import dtoolkit.geoaccessor\n    >>> import pandas as pd\n    >>> import geopandas as gpd\n    >>> from shapely.geometry import Point\n    >>> df = (\n    ...     pd.DataFrame({\"x\": [122, 100], \"y\":[55, 1]})\n    ...     .from_xy(\"x\", \"y\", crs=4326)\n    ... )\n    >>> df\n         x   y                    geometry\n    0  122  55  POINT (122.00000 55.00000)\n    1  100   1   POINT (100.00000 1.00000)\n    >>> df.geodistance(Point(120, 30)) / 1e6\n    0    2.784435\n    1    3.855604\n    dtype: float64\n\n    Calculate the great-circle distance of corresponding points.\n\n    >>> s = gpd.GeoSeries([Point(120, 30), Point(120, 50)], index=[1, 2], crs=4326)\n    >>> s\n    1    POINT (120.00000 30.00000)\n    2    POINT (120.00000 50.00000)\n    dtype: geometry\n    >>> df.geodistance(s) / 1e6\n    0         NaN\n    1    3.855604\n    2         NaN\n    dtype: float64\n    >>> df.geodistance(s, align=False) / 1e6\n    0    2.784435\n    1    5.768885\n    dtype: float64\n    \"\"\"\n\n    if s.crs != 4326:\n        raise ValueError(f\"Only support 'EPSG:4326' CRS, but got {s.crs!r}.\")\n    if not isinstance(other, (BaseGeometry, gpd.base.GeoPandasBase)):\n        raise TypeError(f\"Unknown type: {type(other).__name__!r}.\")\n\n    if isinstance(other, gpd.base.GeoPandasBase):\n        if other.crs != 4326:\n            raise ValueError(f\"Only support 'EPSG:4326' CRS, but got {other.crs!r}.\")\n\n        s = s.geometry\n        if align and not s.index.equals(other.index):\n            warn(\"The indices are different.\", stacklevel=find_stack_level())\n            s, other = s.align(other)\n        else:\n            other = other.geometry\n\n    return pd.Series(\n        distance(\n            s.geometry.x.to_numpy(),\n            s.geometry.y.to_numpy(),\n            other.x if isinstance(other, BaseGeometry) else other.x.to_numpy(),\n            other.y if isinstance(other, BaseGeometry) else other.y.to_numpy(),\n            radius=radius,\n        ),\n        index=s.index,\n    )\n\n\n# based on https://github.com/geopy/geopy geopy/distance.py::great_circle.measure\ndef distance(\n    lng1: np.ndarray | float,\n    lat1: np.ndarray | float,\n    lng2: np.ndarray | float,\n    lat2: np.ndarray | float,\n    radius: float,\n) -> np.ndarray:\n    lng1, lat1, lng2, lat2 = map(np.radians, (lng1, lat1, lng2, lat2))\n    sin_lat1, cos_lat1 = np.sin(lat1), np.cos(lat1)\n    sin_lat2, cos_lat2 = np.sin(lat2), np.cos(lat2)\n\n    delta_lng = lng2 - lng1\n    cos_delta_lng, sin_delta_lng = np.cos(delta_lng), np.sin(delta_lng)\n\n    return radius * np.arctan2(\n        np.sqrt(\n            (cos_lat1 * sin_lat2 - sin_lat1 * cos_lat2 * cos_delta_lng) ** 2\n            + (cos_lat2 * sin_delta_lng) ** 2,\n        ),\n        sin_lat1 * sin_lat2 + cos_lat1 * cos_lat2 * cos_delta_lng,\n    )\n", "path": "dtoolkit/geoaccessor/geoseries/geodistance.py"}, {"content": "from __future__ import annotations\n\nimport geopandas as gpd\nimport numpy as np\nimport pandas as pd\n\nfrom dtoolkit.geoaccessor.register import register_geoseries_method\n\n\n@register_geoseries_method\ndef geodistance_matrix(\n    s: gpd.GeoSeries,\n    /,\n    other: gpd.GeoSeries | gpd.GeoDataFrame | None = None,\n    radius: float = 6371008.7714150598,\n) -> pd.DataFrame:\n    \"\"\"\n    Returns a ``DataFrame`` containing the `great-circle`__ distances matrix between in\n    ``s`` and ``other`` via haversine formula.\n\n    __ https://en.wikipedia.org/wiki/Great-circle_distance\n\n    .. math::\n\n        D(x, y) = 2 \\\\arcsin [\n            \\\\sqrt{\n                \\\\sin^2 ((x_1 - y_1) / 2)\n                + \\\\cos(x_1) \\\\cos(y_1) \\\\sin^2 ((x_2 - y_2) / 2)\n            }\n        ]\n\n    Parameters\n    ----------\n    other : GeoSeries, or GeoDataFrame, default None\n        If None, uses ``other=s``.\n\n    radius : float, default 6371008.7714150598\n        Great-circle distance uses a spherical model of the earth, using the mean earth\n        radius as defined by the International Union of Geodesy and Geophysics,\n        (2\\\\ *a* + *b*)/3 = 6371008.7714150598 meters for WGS-84.\n\n    Returns\n    -------\n    DataFrame\n        - The index and columns are the same as the index of ``s`` and ``other``.\n        - The values are the great-circle distances and its unit is meters.\n\n    Raises\n    ------\n    ModuleNotFoundError\n        If don't have module named 'sklearn'.\n\n    ValueError\n        If the CRS is not ``ESGP:4326``.\n\n    See Also\n    --------\n    sklearn.metrics.pairwise.haversine_distances\n    dtoolkit.geoaccessor.geoseries.geodistance\n    dtoolkit.geoaccessor.geoseries.geodistance_matrix\n    dtoolkit.geoaccessor.geodataframe.geodistance\n    dtoolkit.geoaccessor.geodataframe.geodistance_matrix\n\n    Notes\n    -----\n    - Currently, only supports Point geometry.\n    - The great-circle distance is the angular distance between two points on the\n      surface of a sphere. As the Earth is nearly spherical, the haversine formula\n      provides a good approximation of the distance between two points of the Earth\n      surface, with a less than 1% error on average.\n\n    Examples\n    --------\n    >>> import dtoolkit.geoaccessor\n    >>> df = pd.DataFrame(\n    ...     {\n    ...         \"x\": [120, 122, 100],\n    ...         \"y\":[30, 55, 1],\n    ...     },\n    ... ).from_xy(\"x\", \"y\", crs=4326)\n    >>> df\n         x   y                    geometry\n    0  120  30  POINT (120.00000 30.00000)\n    1  122  55  POINT (122.00000 55.00000)\n    2  100   1   POINT (100.00000 1.00000)\n    >>> other = pd.DataFrame(\n    ...     {\n    ...         \"x\": [120, 110],\n    ...         \"y\":[30, 40],\n    ...     },\n    ... ).from_xy(\"x\", \"y\", crs=4326)\n    >>> other\n         x   y                    geometry\n    0  120  30  POINT (120.00000 30.00000)\n    1  110  40  POINT (110.00000 40.00000)\n    >>> df.geodistance_matrix(other)\n                  0             1\n    0  0.000000e+00  1.203540e+06\n    1  1.439971e+06  1.511958e+06\n    2  2.418544e+06  1.522752e+06\n    \"\"\"\n    from sklearn.metrics.pairwise import haversine_distances\n\n    if s.crs != 4326:\n        raise ValueError(f\"Only support 'EPSG:4326' CRS, but got {s.crs!r}.\")\n    if other is not None and other.crs != 4326:\n        raise ValueError(f\"Only support 'EPSG:4326' CRS, but got {other.crs!r}.\")\n\n    X = np.radians(np.stack((s.x, s.y), axis=1))\n    Y = np.radians(np.stack((other.x, other.y), axis=1)) if other is not None else other\n    return pd.DataFrame(\n        radius * haversine_distances(X, Y),\n        index=s.index,\n        columns=other.index,\n    )\n", "path": "dtoolkit/geoaccessor/geoseries/geodistance_matrix.py"}]}