Repository: mie-lab/trackintel

Issue:
Staypoint generation can have invalid id assignment for positionfixes (causes error in tripleg generation)
When generating staypoints it can happen that a single positionfix in a group of positionfixes that belong to a single staypoint does not get assigned the staypoint id but stays nan. This causes errors in the trip generation as this positionfix leads to the generation of an invalid tripleg with length 1 (The detection mechanism for invalid triplegs depends on the staypoint id). This problem occurs for Geolife users 165, 59, 65, 10, 163, 141. Therefore I am quite sure that this problem has to be related to a recent PR because I have tested the tripleg generation with all geolife users in the past. 

Here is the code to reproduce this problem for geolife user 65: 
```Python
import trackintel as ti
import datetime
import pytz
# the path to a folder that contains geolife user 65
geolife_path = r"xxx"

pfs, mode_labels =ti.io.dataset_reader.read_geolife(geolife_path)
pfs, spts = pfs.as_positionfixes.generate_staypoints(gap_threshold=15, include_last=True, print_progress=True)
start = datetime.datetime(year=2011, month=9, day=3, hour=23, minute=21, second=55, tzinfo=pytz.utc)
end = datetime.datetime(year=2011, month=9, day=3, hour=23, minute=22, second=15, tzinfo=pytz.utc)
pfs_ = pfs[(pfs['tracked_at'] > start) & (pfs['tracked_at'] < end) & (pfs['user_id'] == 65)]

print(pfs_['staypoint_id'])
```
Output: 

|tracked_at | staypoint id |
|---|---|
|234703 2011-09-03 23:21:56+00:00  |        <NA>|
|234704 2011-09-03 23:21:56+00:00  |        <NA>|
|234705 2011-09-03 23:21:56+00:00  |        <NA>|
|234706 2011-09-03 23:21:56+00:00  |        <NA>|
|234707 2011-09-03 23:21:56+00:00  |        <NA>|
|234708 2011-09-03 23:21:57+00:00  |         533|
|234709 2011-09-03 23:21:57+00:00  |         533|
|234710 2011-09-03 23:21:57+00:00  |         533|
|234711 2011-09-03 23:21:57+00:00  |         533|
|234712 2011-09-03 23:21:57+00:00  |        <NA>|
|234713 2011-09-03 23:21:58+00:00  |         533|
|234714 2011-09-03 23:22:14+00:00  |         533|
|234715 2011-09-03 23:22:14+00:00  |         533|
|234716 2011-09-03 23:22:14+00:00  |         533|
|234717 2011-09-03 23:22:14+00:00  |         533|
|234718 2011-09-03 23:22:14+00:00  |         533|



The expected output should be that there is no gap in the staypoint id between positionfixes that belong to the same staypoint. 

Discussion:
hong2223: I encountered these problems when I tried to prepare the GIScience paper's code on the Geolife dataset. These fixes work very well :)
henrymartin1: @hong2223 : Ready for review

Patch:
diff --git a/trackintel/preprocessing/positionfixes.py b/trackintel/preprocessing/positionfixes.py
index e2381b4..ef872d8 100644
--- a/trackintel/preprocessing/positionfixes.py
+++ b/trackintel/preprocessing/positionfixes.py
@@ -1,4 +1,5 @@
 import datetime
+import warnings
 
 import geopandas as gpd
 import numpy as np
@@ -15,9 +16,10 @@ def generate_staypoints(
     distance_metric="haversine",
     dist_threshold=100,
     time_threshold=5.0,
-    gap_threshold=1e6,
+    gap_threshold=15.0,
     include_last=False,
     print_progress=False,
+    exclude_duplicate_pfs=True,
 ):
     """
     Generate staypoints from positionfixes.
@@ -41,18 +43,23 @@ def generate_staypoints(
     time_threshold : float, default 5.0 (minutes)
         The time threshold for the 'sliding' method in minutes.
 
-    gap_threshold : float, default 1e6 (minutes)
-        The time threshold of determine whether a gap exists between consecutive pfs. Staypoints
-        will not be generated between gaps. Only valid in 'sliding' method.
+    gap_threshold : float, default 15.0 (minutes)
+        The time threshold of determine whether a gap exists between consecutive pfs. Consecutive pfs with
+        temporal gaps larger than 'gap_threshold' will be excluded from staypoints generation.
+        Only valid in 'sliding' method.
 
     include_last: boolen, default False
-        The original algorithm (see Li et al. (2008)) only detects staypoint if the user steps out
+        The algorithm in Li et al. (2008) only detects staypoint if the user steps out
         of that staypoint. This will omit the last staypoint (if any). Set 'include_last'
         to True to include this last staypoint.
 
     print_progress: boolen, default False
         Show per-user progress if set to True.
 
+    exclude_duplicate_pfs: boolean, default True
+        Filters duplicate positionfixes before generating staypoints. Duplicates can lead to problems in later
+        processing steps (e.g., when generating triplegs). It is not recommended to set this to False.
+
     Returns
     -------
     pfs: GeoDataFrame (as trackintel positionfixes)
@@ -65,9 +72,9 @@ def generate_staypoints(
     -----
     The 'sliding' method is adapted from Li et al. (2008). In the original algorithm, the 'finished_at'
     time for the current staypoint lasts until the 'tracked_at' time of the first positionfix outside
-    this staypoint. This implies potential tracking gaps may be included in staypoints, and users
-    are assumed to be stationary during this missing period. To avoid including too large gaps, set
-    'gap_threshold' parameter to a small value, e.g., 15 min.
+    this staypoint. Users are assumed to be stationary during this missing period and potential tracking
+    gaps may be included in staypoints. To avoid including too large missing signal gaps, set 'gap_threshold'
+    to a small value, e.g., 15 min.
 
     Examples
     --------
@@ -85,6 +92,17 @@ def generate_staypoints(
     # copy the original pfs for adding 'staypoint_id' column
     pfs = pfs_input.copy()
 
+    if exclude_duplicate_pfs:
+        len_org = pfs.shape[0]
+        pfs = pfs.drop_duplicates()
+        nb_dropped = len_org - pfs.shape[0]
+        if nb_dropped > 0:
+            warn_str = (
+                f"{nb_dropped} duplicates were dropped from your positionfixes. Dropping duplicates is"
+                + " recommended but can be prevented using the 'exclude_duplicate_pfs' flag."
+            )
+            warnings.warn(warn_str)
+
     # if the positionfixes already have a column "staypoint_id", we drop it
     if "staypoint_id" in pfs:
         pfs.drop(columns="staypoint_id", inplace=True)
@@ -169,7 +187,12 @@ def generate_staypoints(
     return pfs, stps
 
 
-def generate_triplegs(pfs_input, stps_input, method="between_staypoints", gap_threshold=15):
+def generate_triplegs(
+    pfs_input,
+    stps_input,
+    method="between_staypoints",
+    gap_threshold=15,
+):
     """Generate triplegs from positionfixes.
 
     Parameters
@@ -365,6 +388,7 @@ def generate_triplegs(pfs_input, stps_input, method="between_staypoints", gap_th
         tpls.crs = pfs.crs
 
         # assert validity of triplegs
+        tpls, pfs = _drop_invalid_triplegs(tpls, pfs)
         tpls.as_triplegs
 
         if case == 2:
@@ -393,7 +417,7 @@ def _generate_staypoints_sliding_user(
     else:
         raise AttributeError("distance_metric unknown. We only support ['haversine']. " f"You passed {distance_metric}")
 
-    df.sort_values("tracked_at", inplace=True)
+    df = df.sort_index(kind="mergesort").sort_values(by=["tracked_at"], kind="mergesort")
     # pfs id should be in index, create separate idx for storing the matching
     pfs = df.to_dict("records")
     idx = df.index.to_list()
@@ -413,7 +437,8 @@ def _generate_staypoints_sliding_user(
             # the duration of gap in the last two pfs
             gap_t = (pfs[curr]["tracked_at"] - pfs[curr - 1]["tracked_at"]).total_seconds()
 
-            # we want the spt to have long duration, but the gap of missing signal should not be too long
+            # we want the staypoint to have long duration,
+            # but the gap of two consecutive positionfixes should not be too long
             if (delta_t >= (time_threshold * 60)) and (gap_t < gap_threshold * 60):
                 new_stps = __create_new_staypoints(start, curr, pfs, idx, elevation_flag, geo_col)
                 # add staypoint
@@ -461,3 +486,43 @@ def __create_new_staypoints(start, end, pfs, idx, elevation_flag, geo_col, last_
     new_stps["pfs_id"] = [idx[k] for k in range(start, end)]
 
     return new_stps
+
+
+def _drop_invalid_triplegs(tpls, pfs):
+    """Remove triplegs with invalid geometries. Also remove the corresponding invalid tripleg ids from positionfixes.
+
+    Parameters
+    ----------
+    tpls : GeoDataFrame (as trackintel triplegs)
+    pfs : GeoDataFrame (as trackintel positionfixes)
+
+    Returns
+    -------
+    tpls: GeoDataFrame (as trackintel triplegs)
+        original tpls with invalid geometries removed.
+
+    pfs: GeoDataFrame (as trackintel positionfixes)
+        original pfs with invalid tripleg id set to pd.NA.
+
+    Notes
+    -----
+    Valid is defined using shapely (https://shapely.readthedocs.io/en/stable/manual.html#object.is_valid) via
+    the geopandas accessor.
+    """
+    invalid_tpls = tpls[~tpls.geometry.is_valid]
+    if invalid_tpls.shape[0] > 0:
+        # identify invalid tripleg ids
+        invalid_tpls_ids = invalid_tpls.index.to_list()
+
+        # reset tpls id in pfs
+        invalid_pfs_ixs = pfs[pfs.tripleg_id.isin(invalid_tpls_ids)].index
+        pfs.loc[invalid_pfs_ixs, "tripleg_id"] = pd.NA
+        warn_string = (
+            f"The positionfixes with ids {invalid_pfs_ixs.values} lead to invalid tripleg geometries. The "
+            f"resulting triplegs were omitted and the tripleg id of the positionfixes was set to nan"
+        )
+        warnings.warn(warn_string)
+
+        # return valid triplegs
+        tpls = tpls[tpls.geometry.is_valid]
+    return tpls, pfs


Test patch:
diff --git a/tests/preprocessing/test_positionfixes.py b/tests/preprocessing/test_positionfixes.py
index e16b9a9..f5a9d14 100644
--- a/tests/preprocessing/test_positionfixes.py
+++ b/tests/preprocessing/test_positionfixes.py
@@ -3,11 +3,11 @@ import os
 import sys
 
 import geopandas as gpd
-from geopandas.testing import assert_geodataframe_equal
-
 import numpy as np
 import pandas as pd
 import pytest
+from geopandas.testing import assert_geodataframe_equal
+from shapely.geometry import Point
 
 import trackintel as ti
 
@@ -20,9 +20,102 @@ def geolife_pfs_stps_long():
     return pfs, stps
 
 
+@pytest.fixture
+def example_positionfixes():
+    """Positionfixes for tests."""
+    p1 = Point(8.5067847, 47.4)
+    p2 = Point(8.5067847, 47.5)
+    p3 = Point(8.5067847, 47.6)
+
+    t1 = pd.Timestamp("1971-01-01 00:00:00", tz="utc")
+    t2 = pd.Timestamp("1971-01-01 05:00:00", tz="utc")
+    t3 = pd.Timestamp("1971-01-02 07:00:00", tz="utc")
+
+    list_dict = [
+        {"user_id": 0, "tracked_at": t1, "geometry": p1},
+        {"user_id": 0, "tracked_at": t2, "geometry": p2},
+        {"user_id": 1, "tracked_at": t3, "geometry": p3},
+    ]
+    pfs = gpd.GeoDataFrame(data=list_dict, geometry="geometry", crs="EPSG:4326")
+    pfs.index.name = "id"
+
+    # assert validity of positionfixes.
+    pfs.as_positionfixes
+    return pfs
+
+
+@pytest.fixture
+def isolated_positionfixes():
+    """Positionfixes with two isolated positionfixes that have the same geometry but different timestamps."""
+    p1 = Point(8.5067847, 47.4)
+    p2 = Point(8.5067847, 47.5)
+    p3 = Point(8.5067847, 47.6)
+    p4 = Point(8.5067847, 47.6)
+
+    t1 = pd.Timestamp("1971-01-01 00:00:00", tz="utc")
+    t2 = pd.Timestamp("1971-01-02 01:01:00", tz="utc")
+    t3 = pd.Timestamp("1971-01-02 01:02:00", tz="utc")
+    t4 = pd.Timestamp("1971-01-03 03:00:00", tz="utc")
+
+    list_dict = [
+        {"user_id": 0, "tracked_at": t1, "geometry": p1, "staypoint_id": 0},
+        {"user_id": 0, "tracked_at": t2, "geometry": p2, "staypoint_id": pd.NA},
+        {"user_id": 0, "tracked_at": t3, "geometry": p2, "staypoint_id": pd.NA},
+        {"user_id": 0, "tracked_at": t4, "geometry": p3, "staypoint_id": 1},
+        {"user_id": 1, "tracked_at": t1, "geometry": p1, "staypoint_id": 2},
+        {"user_id": 1, "tracked_at": t2, "geometry": p2, "staypoint_id": pd.NA},
+        {"user_id": 1, "tracked_at": t3, "geometry": p3, "staypoint_id": pd.NA},
+        {"user_id": 1, "tracked_at": t4, "geometry": p4, "staypoint_id": 3},
+    ]
+    pfs = gpd.GeoDataFrame(data=list_dict, geometry="geometry", crs="EPSG:4326")
+    pfs.index.name = "id"
+
+    # assert validity of positionfixes.
+    pfs.as_positionfixes
+    return pfs
+
+
 class TestGenerate_staypoints:
     """Tests for generate_staypoints() method."""
 
+    def test_duplicate_pfs_warning(self, example_positionfixes):
+        """Calling generate_staypoints with duplicate positionfixes should raise a warning."""
+        pfs_duplicate_loc = example_positionfixes.copy()
+        pfs_duplicate_loc.loc[0, "geometry"] = pfs_duplicate_loc.loc[1, "geometry"]
+
+        pfs_duplicate_t = example_positionfixes.copy()
+        pfs_duplicate_t.loc[0, "tracked_at"] = pfs_duplicate_t.loc[1, "tracked_at"]
+
+        # 0 and 1 pfs are now identical
+        pfs_duplicate_all = pfs_duplicate_loc.copy()
+        pfs_duplicate_all.loc[0, "tracked_at"] = pfs_duplicate_all.loc[1, "tracked_at"]
+
+        with pytest.warns(None) as record:
+            example_positionfixes.as_positionfixes.generate_staypoints()
+            pfs_duplicate_loc.as_positionfixes.generate_staypoints()
+            pfs_duplicate_t.as_positionfixes.generate_staypoints()
+
+            # assert that no warning is raised
+            assert len(record) == 0
+
+        warn_string = ".* duplicates were dropped from your positionfixes."
+        with pytest.warns(UserWarning, match=warn_string):
+            pfs_duplicate_all.as_positionfixes.generate_staypoints()
+
+    def test_duplicate_pfs_filtering(self, example_positionfixes):
+        """Test that duplicate positionfixes are filtered in generate_staypoints."""
+        pfs = example_positionfixes
+
+        # set 0 and 1 pfs identical
+        pfs.loc[0, "geometry"] = pfs.loc[1, "geometry"]
+        pfs.loc[0, "tracked_at"] = pfs.loc[1, "tracked_at"]
+
+        with pytest.warns(UserWarning):
+            pfs_out, _ = pfs.as_positionfixes.generate_staypoints()
+
+        # drop staypoint_id column of pfs_out and ensure that the second duplicate (id=1) is filtered
+        assert_geodataframe_equal(pfs_out.drop(columns="staypoint_id"), pfs.iloc[[0, 2]])
+
     def test_duplicate_columns(self):
         """Test if running the function twice, the generated column does not yield exception in join statement"""
 
@@ -110,33 +203,44 @@ class TestGenerate_staypoints:
 
         # the duration should be not longer than time_threshold
         time_threshold_ls = [3, 5, 10]
+
+        # the missing time should not exceed gap_threshold
+        gap_threshold_ls = [5, 10, 15, 20]
         for time_threshold in time_threshold_ls:
-            _, stps = pfs_input.as_positionfixes.generate_staypoints(
-                method="sliding", dist_threshold=50, time_threshold=time_threshold
-            )
+            for gap_threshold in gap_threshold_ls:
+                pfs, stps = pfs_input.as_positionfixes.generate_staypoints(
+                    method="sliding", dist_threshold=50, time_threshold=time_threshold, gap_threshold=gap_threshold
+                )
 
-            duration_stps_min = (stps["finished_at"] - stps["started_at"]).dt.total_seconds() / 60
-            # all durations should be longer than the time_threshold
-            assert (duration_stps_min > time_threshold).all()
+                # all durations should be longer than the time_threshold
+                duration_stps_min = (stps["finished_at"] - stps["started_at"]).dt.total_seconds() / 60
+                assert (duration_stps_min > time_threshold).all()
 
-        # the missing time should not exceed gap_threshold
-        gap_threshold_ls = [10, 15, 20]
-        for gap_threshold in gap_threshold_ls:
-            pfs, _ = pfs_input.as_positionfixes.generate_staypoints(
-                method="sliding", dist_threshold=50, time_threshold=time_threshold, gap_threshold=gap_threshold
-            )
-            # get the difference between pfs tracking time, and assign back to the previous pfs
-            pfs["diff"] = ((pfs["tracked_at"] - pfs["tracked_at"].shift(1)).dt.total_seconds() / 60).shift(-1)
-            # get the last pf of stps and check the gap size
-            pfs.dropna(subset=["staypoint_id"], inplace=True)
-            pfs.drop_duplicates(subset=["staypoint_id"], keep="last", inplace=True)
-            # all last pfs should be shorter than the gap_threshold
-            assert (pfs["diff"] < gap_threshold).all()
+                # all last pfs should be shorter than the gap_threshold
+                # get the difference between pfs tracking time, and assign back to the previous pfs
+                pfs["diff"] = ((pfs["tracked_at"] - pfs["tracked_at"].shift(1)).dt.total_seconds() / 60).shift(-1)
+                # get the last pf of stps and check the gap size
+                pfs.dropna(subset=["staypoint_id"], inplace=True)
+                pfs.drop_duplicates(subset=["staypoint_id"], keep="last", inplace=True)
+                assert (pfs["diff"] < gap_threshold).all()
 
 
 class TestGenerate_triplegs:
     """Tests for generate_triplegs() method."""
 
+    def test_invalid_isolates(self, isolated_positionfixes):
+        """Triplegs generated from isolated duplicates are dropped."""
+        pfs = isolated_positionfixes
+
+        # a warning shall raise due to the duplicated positionfixes
+        warn_string = "The positionfixes with ids .* lead to invalid tripleg geometries."
+        with pytest.warns(UserWarning, match=warn_string):
+            pfs, tpls = pfs.as_positionfixes.generate_triplegs()
+
+        # tripleg of user 0 is dropped. Tripleg of user 1 is available.
+        assert len(tpls) == 1
+        assert tpls.user_id.iloc[0] == 1
+
     def test_duplicate_columns(self, geolife_pfs_stps_long):
         """Test if running the function twice, the generated column does not yield exception in join statement"""
 
diff --git a/tests/preprocessing/test_staypoints.py b/tests/preprocessing/test_staypoints.py
index f1c5339..691a3e4 100644
--- a/tests/preprocessing/test_staypoints.py
+++ b/tests/preprocessing/test_staypoints.py
@@ -108,7 +108,9 @@ class TestGenerate_locations:
         """Test with small epsilon parameter."""
         pfs_file = os.path.join("tests", "data", "positionfixes.csv")
         pfs = ti.read_positionfixes_csv(pfs_file, sep=";", tz="utc", index_col="id")
-        _, stps = pfs.as_positionfixes.generate_staypoints(method="sliding", dist_threshold=0, time_threshold=0)
+        _, stps = pfs.as_positionfixes.generate_staypoints(
+            method="sliding", gap_threshold=1e6, dist_threshold=0, time_threshold=0
+        )
         _, locs_user = stps.as_staypoints.generate_locations(
             method="dbscan", epsilon=1e-18, num_samples=0, agg_level="user"
         )
@@ -123,7 +125,9 @@ class TestGenerate_locations:
         """Test with large epsilon parameter."""
         pfs_file = os.path.join("tests", "data", "positionfixes.csv")
         pfs = ti.read_positionfixes_csv(pfs_file, sep=";", tz="utc", index_col="id")
-        _, stps = pfs.as_positionfixes.generate_staypoints(method="sliding", dist_threshold=0, time_threshold=0)
+        _, stps = pfs.as_positionfixes.generate_staypoints(
+            method="sliding", gap_threshold=1e6, dist_threshold=0, time_threshold=0
+        )
         _, locs_user = stps.as_staypoints.generate_locations(
             method="dbscan", epsilon=1e18, num_samples=1000, agg_level="user"
         )
@@ -138,7 +142,9 @@ class TestGenerate_locations:
         """Test nan is assigned for missing link between stps and locs."""
         pfs_file = os.path.join("tests", "data", "positionfixes.csv")
         pfs = ti.read_positionfixes_csv(pfs_file, sep=";", tz="utc", index_col="id")
-        _, stps = pfs.as_positionfixes.generate_staypoints(method="sliding", dist_threshold=0, time_threshold=0)
+        _, stps = pfs.as_positionfixes.generate_staypoints(
+            method="sliding", gap_threshold=1e6, dist_threshold=0, time_threshold=0
+        )
         stps, _ = stps.as_staypoints.generate_locations(
             method="dbscan", epsilon=1e18, num_samples=1000, agg_level="user"
         )
diff --git a/tests/preprocessing/test_triplegs.py b/tests/preprocessing/test_triplegs.py
index 80e8279..5d59730 100644
--- a/tests/preprocessing/test_triplegs.py
+++ b/tests/preprocessing/test_triplegs.py
@@ -56,7 +56,9 @@ class TestGenerate_trips:
 
         # create trips from geolife (based on positionfixes)
         pfs, _ = ti.io.dataset_reader.read_geolife(os.path.join("tests", "data", "geolife_long"))
-        pfs, stps = pfs.as_positionfixes.generate_staypoints(method="sliding", dist_threshold=25, time_threshold=5)
+        pfs, stps = pfs.as_positionfixes.generate_staypoints(
+            method="sliding", dist_threshold=25, time_threshold=5, gap_threshold=1e6
+        )
         stps = stps.as_staypoints.create_activity_flag(time_threshold=15)
         pfs, tpls = pfs.as_positionfixes.generate_triplegs(stps)
 
@@ -231,7 +233,9 @@ class TestGenerate_trips:
 
         # create trips from geolife (based on positionfixes)
         pfs, _ = ti.io.dataset_reader.read_geolife(os.path.join("tests", "data", "geolife_long"))
-        pfs, stps = pfs.as_positionfixes.generate_staypoints(method="sliding", dist_threshold=25, time_threshold=5)
+        pfs, stps = pfs.as_positionfixes.generate_staypoints(
+            method="sliding", dist_threshold=25, time_threshold=5, gap_threshold=1e6
+        )
         stps = stps.as_staypoints.create_activity_flag(time_threshold=15)
         pfs, tpls = pfs.as_positionfixes.generate_triplegs(stps)