Repository: mie-lab/trackintel
Pull request #563 (2023-10-30 22:50:58)

ENH: add radius of gyration
I did not know where to put the function, so I put it so far into a new file in `analysis`.
Added branch to calculate planar distances as well. 

bifbof: PR updated @hongyeehh :)

Patch:
diff --git a/docs/modules/analysis.rst b/docs/modules/analysis.rst
index c3e9c5e..0375a31 100644
--- a/docs/modules/analysis.rst
+++ b/docs/modules/analysis.rst
@@ -30,4 +30,9 @@ Location Identification
 
 .. autofunction:: trackintel.analysis.location_identification.freq_method
 
-.. autofunction:: trackintel.analysis.location_identification.osna_method
\ No newline at end of file
+.. autofunction:: trackintel.analysis.location_identification.osna_method
+
+Metrics
+=======
+
+.. autofunction:: trackintel.analysis.radius_gyration
diff --git a/trackintel/analysis/__init__.py b/trackintel/analysis/__init__.py
index 987be42..6cabcc3 100644
--- a/trackintel/analysis/__init__.py
+++ b/trackintel/analysis/__init__.py
@@ -10,6 +10,8 @@ from .location_identification import location_identifier
 from .location_identification import pre_filter_locations
 from .location_identification import freq_method, osna_method
 
+from .metrics import radius_gyration
+
 __all__ = [
     "temporal_tracking_quality",
     "split_overlaps",
@@ -20,4 +22,5 @@ __all__ = [
     "pre_filter_locations",
     "freq_method",
     "osna_method",
+    "radius_gyration",
 ]
diff --git a/trackintel/analysis/metrics.py b/trackintel/analysis/metrics.py
new file mode 100644
index 0000000..0398ab1
--- /dev/null
+++ b/trackintel/analysis/metrics.py
@@ -0,0 +1,77 @@
+import numpy as np
+from tqdm import tqdm
+
+from trackintel.geogr import point_haversine_dist, check_gdf_planar
+
+
+def radius_gyration(sp, method="count", print_progress=False):
+    """
+    Radius of gyration for individual users.
+
+    Parameters
+    ----------
+    sp : Staypoints
+
+    method: string, {"count", "duration"}
+        Weighting for center of mass and average distance calculation.
+        - "count": assigns each Point the same weight of 1.
+        - "duration": assigns each Point a weight based on duration.
+           Additionally needs Timedelta column "duration" in sp.
+
+    print_progress: bool, default False
+        Show per-user progress if set to True.
+
+    Returns
+    -------
+    Series
+        Radius of gyration for individual users.
+
+    References
+    ----------
+    [1] Gonzalez, M. C., Hidalgo, C. A., & Barabasi, A. L. (2008).
+    Understanding individual human mobility patterns. Nature, 453(7196), 779-782.
+
+    """
+    if method not in ["count", "duration"]:
+        raise ValueError(f'Method unknown. Should be on of {{"count", "duration"}}. You passed "{method}"')
+
+    if print_progress:
+        tqdm.pandas(desc="User radius of gyration calculation")
+        s = sp.groupby("user_id").progress_apply(_radius_gyration_user, method=method)
+    else:
+        s = sp.groupby("user_id").apply(_radius_gyration_user, method=method)
+
+    s = s.rename("radius_gyration")
+    return s
+
+
+def _radius_gyration_user(sp, method):
+    """
+    User level radius of gyration calculation, see radius_gyration() for more details.
+
+    Parameters
+    ----------
+    sp : Staypoints
+    method: {"count", "duration"}
+
+    Returns
+    -------
+    float
+        The radius of gyration of the user
+    """
+    x = sp.geometry.x
+    y = sp.geometry.y
+
+    if method == "duration":
+        w = sp["duration"].dt.total_seconds()
+    else:  # method == count (check is done in upper level)
+        w = np.ones_like(x)
+
+    x_center = np.average(x, weights=w)
+    y_center = np.average(y, weights=w)
+    if check_gdf_planar(sp):
+        sq_dist = (x - x_center) ** 2 + (y - y_center) ** 2
+    else:
+        sq_dist = point_haversine_dist(x, y, x_center, y_center) ** 2
+    square_rg = np.average(sq_dist, weights=w)
+    return np.sqrt(square_rg)
diff --git a/trackintel/model/staypoints.py b/trackintel/model/staypoints.py
index 1fc7286..4dfac3c 100644
--- a/trackintel/model/staypoints.py
+++ b/trackintel/model/staypoints.py
@@ -164,3 +164,11 @@ class Staypoints(TrackintelBase, TrackintelGeoDataFrame):
         return ti.preprocessing.triplegs.generate_trips(
             self, triplegs, gap_threshold=gap_threshold, add_geometry=add_geometry
         )
+
+    def radius_gyration(self, method="count", print_progress=False):
+        """
+        Calculate radius for gyration for Staypoints
+
+        See :func:`trackintel.analysis.radius_gyration` for full documentation.
+        """
+        return ti.analysis.radius_gyration(self, method, print_progress)


Test patch:
diff --git a/tests/analysis/test_metrics.py b/tests/analysis/test_metrics.py
new file mode 100644
index 0000000..5b5f86f
--- /dev/null
+++ b/tests/analysis/test_metrics.py
@@ -0,0 +1,82 @@
+import pytest
+import numpy as np
+import pandas as pd
+from pandas.testing import assert_series_equal
+from shapely.geometry import Point
+
+import trackintel as ti
+from trackintel.analysis import radius_gyration
+from trackintel.geogr import point_haversine_dist
+
+
+@pytest.fixture
+def staypoints():
+    """Staypoints for two users.
+    First user has 5 unique points all with the same degree distance from each other, but different timestamps.
+    Second user has 2 unique points, but duration of one point is zero.
+    """
+    p1 = Point(0.0, 9.0)
+    p2 = Point(3.0, 12.0)
+    p3 = Point(6.0, 15.0)
+
+    t = pd.Timestamp("1971-01-01 00:00:00", tz="utc")
+    h = pd.Timedelta(hours=1)
+
+    list_dict = [
+        {"id": 1, "user_id": 0, "started_at": t, "finished_at": t + 1 * h, "geom": p1},
+        {"id": 2, "user_id": 0, "started_at": t, "finished_at": t + 2 * h, "geom": p2},
+        {"id": 3, "user_id": 0, "started_at": t, "finished_at": t + 1 * h, "geom": p3},
+        {"id": 4, "user_id": 1, "started_at": t, "finished_at": t + 1 * h, "geom": p1},
+        {"id": 7, "user_id": 1, "started_at": t, "finished_at": t + 1 * h, "geom": p1},
+        {"id": 8, "user_id": 1, "started_at": t, "finished_at": t + 0 * h, "geom": p3},
+    ]
+    sp = ti.Staypoints(data=list_dict, geometry="geom", crs="EPSG:2056")
+    return sp
+
+
+class TestRadius_gyration:
+    def test_unknown_method(self, staypoints):
+        """Test if unknown method raises a ValueError"""
+        with pytest.raises(
+            ValueError, match='Method unknown. Should be on of {"count", "duration"}. You passed "unknown"'
+        ):
+            radius_gyration(staypoints, method="unknown")
+
+    def test_count(self, staypoints):
+        """Test count-method with planar crs"""
+        s = radius_gyration(staypoints, method="count")
+        v1 = np.sqrt((4 * 3**2) / 3)  # center is (3, 12)
+        v2 = np.sqrt((2 * (2**3 + 4**2)) / 3)  # center is (2, 11)
+        assert_series_equal(s, pd.Series([v1, v2]), check_index=False, check_names=False)
+
+    def test_duration(self, staypoints):
+        """Test duration-method with a planar crs."""
+        staypoints["duration"] = staypoints["finished_at"] - staypoints["started_at"]
+        s = radius_gyration(staypoints, method="duration")
+        v1 = np.sqrt(4 * 3**2 / 4)  # center is (3, 12) weight is 4 due to duration
+        v2 = 0  # center lies on p1 and only p1 remains -> 0 variance
+        assert_series_equal(s, pd.Series([v1, v2]), check_index=False, check_names=False)
+
+    def test_haversine(self, staypoints):
+        """Test haversine distance calculations"""
+        staypoints = staypoints.set_crs(4326, allow_override=True)
+        s = radius_gyration(staypoints, method="count")
+        x = staypoints.geometry.x
+        y = staypoints.geometry.y
+        f1 = staypoints["user_id"] == 0
+        f2 = staypoints["user_id"] == 1
+        d1 = point_haversine_dist(x[f1], y[f1], 3, 12)
+        d2 = point_haversine_dist(x[f2], y[f2], 2, 11)
+        v1 = np.sqrt(np.mean(d1**2))
+        v2 = np.sqrt(np.mean(d2**2))
+        assert_series_equal(s, pd.Series([v1, v2]), check_index=False, check_names=False)
+
+    def test_tqdm(self, staypoints):
+        """Test if tqdm works fine"""
+        radius_gyration(staypoints, print_progress=True)
+
+    def test_staypoints_method(self, staypoints):
+        """Test if staypoint method returns same result"""
+        sfunc = radius_gyration(staypoints)
+        smeth = staypoints.radius_gyration()
+        assert_series_equal(sfunc, smeth)