Repository: mie-lab/trackintel
Pull request #66 (2020-11-20 21:46:20)

feat(filter): add spatial filter function
- add `spatial filter` function under `filter` module.
- register the function under the `staypoints`, `locations`, and `triplegs` accessors.
- add tests: the results are compared with ArcGIS manually obtained results

codecov[bot]: # [Codecov](https://codecov.io/gh/mie-lab/trackintel/pull/66?src=pr&el=h1) Report
> Merging [#66](https://codecov.io/gh/mie-lab/trackintel/pull/66?src=pr&el=desc) (80c5d4c) into [master](https://codecov.io/gh/mie-lab/trackintel/commit/734a46f98e6ff7b74a79e25a0ebb8929fb732521?el=desc) (734a46f) will **increase** coverage by `0.85%`.
> The diff coverage is `83.87%`.

[![Impacted file tree graph](https://codecov.io/gh/mie-lab/trackintel/pull/66/graphs/tree.svg?width=650&height=150&src=pr&token=ITDZEEE1MO)](https://codecov.io/gh/mie-lab/trackintel/pull/66?src=pr&el=tree)

```diff
@@            Coverage Diff             @@
##           master      #66      +/-   ##
==========================================
+ Coverage   51.56%   52.42%   +0.85%     
==========================================
  Files          29       30       +1     
  Lines        1247     1280      +33     
  Branches      201      207       +6     
==========================================
+ Hits          643      671      +28     
- Misses        542      544       +2     
- Partials       62       65       +3     
```


| [Impacted Files](https://codecov.io/gh/mie-lab/trackintel/pull/66?src=pr&el=tree) | Coverage Δ | |
|---|---|---|
| [trackintel/preprocessing/filter.py](https://codecov.io/gh/mie-lab/trackintel/pull/66/diff?src=pr&el=tree#diff-dHJhY2tpbnRlbC9wcmVwcm9jZXNzaW5nL2ZpbHRlci5weQ==) | `79.16% <79.16%> (ø)` | |
| [trackintel/model/locations.py](https://codecov.io/gh/mie-lab/trackintel/pull/66/diff?src=pr&el=tree#diff-dHJhY2tpbnRlbC9tb2RlbC9sb2NhdGlvbnMucHk=) | `75.00% <100.00%> (+3.57%)` | :arrow_up: |
| [trackintel/model/staypoints.py](https://codecov.io/gh/mie-lab/trackintel/pull/66/diff?src=pr&el=tree#diff-dHJhY2tpbnRlbC9tb2RlbC9zdGF5cG9pbnRzLnB5) | `84.84% <100.00%> (+1.51%)` | :arrow_up: |
| [trackintel/model/triplegs.py](https://codecov.io/gh/mie-lab/trackintel/pull/66/diff?src=pr&el=tree#diff-dHJhY2tpbnRlbC9tb2RlbC90cmlwbGVncy5weQ==) | `81.48% <100.00%> (+2.31%)` | :arrow_up: |

------

[Continue to review full report at Codecov](https://codecov.io/gh/mie-lab/trackintel/pull/66?src=pr&el=continue).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/mie-lab/trackintel/pull/66?src=pr&el=footer). Last update [734a46f...80c5d4c](https://codecov.io/gh/mie-lab/trackintel/pull/66?src=pr&el=lastupdated). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments).


Patch:
diff --git a/trackintel/model/locations.py b/trackintel/model/locations.py
index 96579b2..bc0f5e5 100644
--- a/trackintel/model/locations.py
+++ b/trackintel/model/locations.py
@@ -3,6 +3,7 @@ import trackintel as ti
 
 import trackintel.visualization.staypoints
 import trackintel.visualization.locations
+import trackintel.preprocessing.filter
 
 
 @pd.api.extensions.register_dataframe_accessor("as_locations")
@@ -56,3 +57,8 @@ class LocationsAccessor(object):
         See :func:`trackintel.io.postgis.write_locations_postgis`."""
         ti.io.postgis.write_locations_postgis(self._obj, conn_string, table_name, 
             schema, sql_chunksize, if_exists)
+        
+    def spatial_filter(self, *args, **kwargs):
+        """Filter locations with a geo extent.
+        See :func:`trackintel.preprocessing.filter.spatial_filter`."""
+        return ti.preprocessing.filter.spatial_filter(self._obj, *args, **kwargs)
\ No newline at end of file
diff --git a/trackintel/model/staypoints.py b/trackintel/model/staypoints.py
index 937b7bd..1e93cc2 100644
--- a/trackintel/model/staypoints.py
+++ b/trackintel/model/staypoints.py
@@ -2,6 +2,8 @@ import pandas as pd
 
 import trackintel as ti
 import trackintel.preprocessing.staypoints
+import trackintel.preprocessing.filter
+
 import trackintel.visualization.staypoints
 
 
@@ -67,7 +69,10 @@ class StaypointsAccessor(object):
         See :func:`trackintel.preprocessing.staypoints.create_activity_flag`."""
         return ti.preprocessing.staypoints.create_activity_flag(self._obj, *args, **kwargs)
 
-
+    def spatial_filter(self, *args, **kwargs):
+        """Filter staypoints with a geo extent.
+        See :func:`trackintel.preprocessing.filter.spatial_filter`."""
+        return ti.preprocessing.filter.spatial_filter(self._obj, *args, **kwargs)
 
     def plot(self, *args, **kwargs):
         """Plots this collection of staypoints. 
diff --git a/trackintel/model/triplegs.py b/trackintel/model/triplegs.py
index 95f45df..f53746d 100644
--- a/trackintel/model/triplegs.py
+++ b/trackintel/model/triplegs.py
@@ -2,6 +2,7 @@ import pandas as pd
 
 import trackintel as ti
 import trackintel.visualization.triplegs
+import trackintel.preprocessing.filter
 
 
 @pd.api.extensions.register_dataframe_accessor("as_triplegs")
@@ -67,4 +68,9 @@ class TriplegsAccessor(object):
         """Calculate pair-wise distance among triplegs (x) or to other triplegs (y)
         See :func:`trackintel.geogr.distances.calculate_distance_matrix`.
         """
-        return ti.geogr.distances.calculate_distance_matrix(self._obj, *args, **kwargs)
\ No newline at end of file
+        return ti.geogr.distances.calculate_distance_matrix(self._obj, *args, **kwargs)
+    
+    def spatial_filter(self, *args, **kwargs):
+        """Filter triplegs with a geo extent.
+        See :func:`trackintel.preprocessing.filter.spatial_filter`."""
+        return ti.preprocessing.filter.spatial_filter(self._obj, *args, **kwargs)
\ No newline at end of file
diff --git a/trackintel/preprocessing/filter.py b/trackintel/preprocessing/filter.py
new file mode 100644
index 0000000..df92dfe
--- /dev/null
+++ b/trackintel/preprocessing/filter.py
@@ -0,0 +1,74 @@
+import geopandas as gpd
+
+
+
+def spatial_filter(source, areas, method="within", re_project=False):
+    """Filter staypoints, locations or triplegs with a geo extent. 
+
+    Parameters
+    ----------
+    source : GeoDataFrame
+        The source feature to perform the spatial filtering 
+
+    areas : GeoDataFrame
+        The areas used to perform the spatial filtering. Note, you can have multiple Polygons 
+        and it will return all the features intersect with ANY of those geometries.
+
+    method : str, one of {'within', 'intersects', 'crosses'}, default 'within'
+        'within'    : return instances in 'source' where no points of these instances lies in the 
+            exterior of the 'areas' and at least one point of the interior of these instances lies 
+            in the interior of 'areas'.
+            
+        'intersects': return instances in 'source' where the boundary or interior of these instances 
+            intersect in any way with those of the 'areas'
+            
+        'crosses'   : return instances in 'source' where the interior of these instances intersects 
+            the interior of the 'areas' but does not contain it, and the dimension of the intersection 
+            is less than the dimension of the one of the 'areas'.
+            
+    re_project : bool, default False
+        If this is set to True, the 'source' will be projected to the coordinate reference system of 'areas' 
+    
+    Returns
+    -------
+    GeoDataFrame
+        A new GeoDataFrame containing the features after the spatial filtering.
+        
+    Examples
+    --------
+    >>> spts.as_staypoints.spatial_filter(areas, method="within", re_project=False)
+    """
+    
+    gdf = source.copy()
+    
+    if re_project:
+        init_crs = gdf.crs
+        gdf = gdf.to_crs(areas.crs)
+    
+    # build spatial index for pre filtering
+    source_sindex = gdf.sindex
+    possible_matches_index = []
+    for other in areas.itertuples():
+        bounds = other.geometry.bounds
+        c = list(source_sindex.intersection(bounds))
+        possible_matches_index += c
+
+    # Get unique candidates
+    unique_candidate_matches = list(set(possible_matches_index))
+    possible_matches = gdf.iloc[unique_candidate_matches]
+    
+    # get final result
+    if method == 'within':
+        ret_gdf = possible_matches.loc[possible_matches.within(areas.unary_union)]
+    elif method == 'intersects':
+        ret_gdf = possible_matches.loc[possible_matches.intersects(areas.unary_union)]
+    elif method == 'crosses':
+        ret_gdf = possible_matches.loc[possible_matches.crosses(areas.unary_union)]
+    else:
+        raise AttributeError("method unknown. We only support ['within', 'intersects', 'crosses']. "
+                                f"You passed {method}")
+        
+    if re_project:
+        return ret_gdf.to_crs(init_crs)
+    else:
+        return ret_gdf
\ No newline at end of file


Test patch:
diff --git a/tests/data/area/tsinghua.geojson b/tests/data/area/tsinghua.geojson
new file mode 100644
index 0000000..a80335c
--- /dev/null
+++ b/tests/data/area/tsinghua.geojson
@@ -0,0 +1,7 @@
+{
+"type": "FeatureCollection",
+"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:EPSG::32649" } },
+"features": [
+{ "type": "Feature", "properties": { }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 954219.888899999670684, 4442967.642000000923872 ], [ 954386.277999999932945, 4442542.3092 ], [ 954945.964599999599159, 4441979.919399999082088 ], [ 954996.083200000226498, 4441528.671399999409914 ], [ 955040.927600000053644, 4441008.1053 ], [ 954724.604399999603629, 4440586.730499999597669 ], [ 954831.445299999788404, 4440316.122700000181794 ], [ 954476.167700000107288, 4440326.947799999266863 ], [ 953910.890100000426173, 4440237.836100000888109 ], [ 953781.917200000025332, 4440398.779799999669194 ], [ 953775.824599999934435, 4440783.211200000718236 ], [ 953461.738400000147521, 4440766.668800000101328 ], [ 953289.028699999675155, 4440846.307800000533462 ], [ 953228.022499999962747, 4441168.712200000882149 ], [ 953349.642699999734759, 4441263.926200000569224 ], [ 953398.595800000242889, 4441414.5483 ], [ 953254.71999999973923, 4441357.575400000438094 ], [ 953210.768600000068545, 4441572.42710000090301 ], [ 953473.077800000086427, 4441625.664300000295043 ], [ 953591.298100000247359, 4441918.094399999827147 ], [ 953442.578200000338256, 4442374.0351 ], [ 953774.034900000318885, 4442569.08699999935925 ], [ 954219.888899999670684, 4442967.642000000923872 ] ] ] } }
+]
+}
diff --git a/tests/preprocessing/test_filter.py b/tests/preprocessing/test_filter.py
new file mode 100644
index 0000000..c1853fe
--- /dev/null
+++ b/tests/preprocessing/test_filter.py
@@ -0,0 +1,77 @@
+import pytest
+import os
+import geopandas as gpd
+
+import trackintel as ti
+
+
+class TestFilter():
+    
+    def test_filter_staypoints(self):
+        # read staypoints and area file
+        spts = ti.read_staypoints_csv(os.path.join('tests', 'data', 'geolife', 'geolife_staypoints.csv'))
+        extent = gpd.read_file(os.path.join('tests', 'data', 'area', 'tsinghua.geojson'))
+        
+        # the projection needs to be defined: WGS84
+        spts.crs = 'epsg:4326'
+        within_spts = spts.as_staypoints.spatial_filter(areas=extent, method="within", re_project=True)
+        intersects_spts = spts.as_staypoints.spatial_filter(areas=extent, method="intersects", re_project=True)
+        crosses_spts = spts.as_staypoints.spatial_filter(areas=extent, method="crosses", re_project=True)
+        
+        # the result obtained from ArcGIS
+        gis_within_num = 13
+        
+        assert len(within_spts) == gis_within_num, "The spatial filtered sp number should be the same as" + \
+            "the one from the result with ArcGIS"
+        assert all(within_spts.geometry == intersects_spts.geometry), "For sp the result of within and" + \
+            "intersects should be the same"
+        assert len(crosses_spts) == 0, "There will be no point crossing area"
+        
+    def test_filter_triplegs(self):
+        # read triplegs and area file
+        tl = ti.read_triplegs_csv(os.path.join('tests', 'data', 'geolife', 'geolife_triplegs.csv'))
+        extent = gpd.read_file(os.path.join('tests', 'data', 'area', 'tsinghua.geojson'))
+        
+        # the projection needs to be defined: WGS84
+        tl.crs = 'epsg:4326'
+        within_tl = tl.as_triplegs.spatial_filter(areas=extent, method="within", re_project=True)
+        intersects_tl = tl.as_triplegs.spatial_filter(areas=extent, method="intersects", re_project=True)
+        crosses_tl = tl.as_triplegs.spatial_filter(areas=extent, method="crosses", re_project=True)
+        
+        # the result obtained from ArcGIS
+        gis_within_num = 9
+        gis_intersects_num = 20
+        
+        assert len(within_tl) == gis_within_num, "The within tripleg number should be the same as" + \
+            "the one from the result with ArcGIS"
+        assert len(intersects_tl) == gis_intersects_num, "The intersects tripleg number should be " + \
+            "the same as the one from the result with ArcGIS"
+        assert len(crosses_tl) == len(intersects_tl) - len(within_tl), "The crosses tripleg number" + \
+            "should equal the number of intersect triplegs minus the number of within triplegs"
+    
+    def test_filter_locations(self):
+        # read staypoints and area file
+        spts = ti.read_staypoints_csv(os.path.join('tests', 'data', 'geolife', 'geolife_staypoints.csv'))
+        extent = gpd.read_file(os.path.join('tests', 'data', 'area', 'tsinghua.geojson'))
+        
+        # cluster staypoints to locations
+        _, locs = spts.as_staypoints.extract_locations(method='dbscan', epsilon=10, 
+                                                       num_samples=0, distance_matrix_metric='haversine',
+                                                       agg_level='dataset')
+        
+        # the projection needs to be defined: WGS84
+        locs.crs = 'epsg:4326'
+        
+        # filter locations with the area
+        within_loc = locs.as_locations.spatial_filter(areas=extent, method="within", re_project=True)
+        intersects_loc = locs.as_locations.spatial_filter(areas=extent, method="intersects", re_project=True)
+        crosses_loc = locs.as_locations.spatial_filter(areas=extent, method="crosses", re_project=True)
+        
+        # the result obtained from ArcGIS
+        gis_within_num = 12
+        
+        assert len(within_loc) == gis_within_num, "The spatial filtered location number should be the same as" + \
+            "the one from the result with ArcGIS"
+        assert all(within_loc.geometry == intersects_loc.geometry), "For location the result of within and" + \
+            "intersects should be the same"
+        assert len(crosses_loc) == 0, "There will be no point crossing area"
\ No newline at end of file