Added geospatial plus test cases

README.md

@@ -12,7 +12,7 @@ Either clone the repo directly or use pip install.
 
 Make sure to replace `<ssh-link>` with the correct URL.
 
-# Requirements 
+# Requirements
 
 To install all necessary requirements use:
 
@@ -22,6 +22,10 @@ To install all necessary requirements use:
 
 # Testing
 
-To run all available unit-tests: 
+To run all available unit-tests:
 
-`python -m pytest tests`
\ No newline at end of file
+`python -m pytest tests`
+
+to also output stdout use:
+
+`python -m pytest -s tests`

measprocess/__init__.py

 from .preprocess import link_dataframes
+from .geospatial import get_geoseries_streets, get_geoseries_blockages
+from .geospatial import project_onto_streets

measprocess/geospatial.py

-import pandas as pd
-import geopandas as gpd
 import overpy
+import geopandas as gpd
+from shapely.geometry import LineString, Polygon
+import networkx as nx
+from scipy.spatial.distance import cdist
+import numpy as np
+import warnings
+import time
+
+def make_overpy_request(request_body, retries):
+    for _ in range(retries):
+        try:
+            api = overpy.Overpass()
+            result = api.query(
+                request_body
+            )
+        except (overpy.exception.OverpassTooManyRequests, overpy.exception.OverpassGatewayTimeout):
+            warnings.warn("Overpass Too Many Requests --- Retry in 5 seconds")
+            time.sleep(5)
+        else:
+            return result
+
+def get_geoseries_streets(measurement_coords, retries=5):
+    '''
+    Obtain the street shapes in the area spanned by the measurement_coords.
+
+    :param measurement_coords: A geopandas geoseries of Points(lon, lat) representing the measurement positions in a EPSG:4326 projection
+
+    :return: A geopandas geoseries consiting of LineStrings representing the streets in EPSG:4326
+    '''
+
+    if measurement_coords.crs != "EPSG:4326":
+        raise ValueError("Make sure to pass data with EPSG:4326 projection")
+
+    long, lat = measurement_coords.x, measurement_coords.y
+    bounds = (lat.min(), long.min(), lat.max(), long.max())
+
+    result = make_overpy_request(
+        f"""
+        [out:json][timeout:25];
+        (
+        way['highway']['highway'!='footway']['highway'!='platform']['highway'!='steps']{bounds};
+        );
+        out body;
+        >;
+        out skel qt;
+        """,
+        retries
+    )
+
+    street_series = gpd.GeoSeries(
+        LineString(
+            ((node.lon, node.lat) for node in way.nodes)
+        ) for way in result.ways
+    ).set_crs("EPSG:4326")
+
+    return street_series
+
+def get_geoseries_blockages(measurement_coords, retries=5):
+    '''
+    Obtain the blockage shapes in the area spanned by the measurement_coords.
+
+    :param measurement_coords: A geopandas geoseries of Points(lon, lat) representing the measurement positions in a EPSG:4326 projection
+
+    :return: A geopandas geoseries consiting of Polygons representing the blockages in EPSG:4326
+    '''
+
+    if measurement_coords.crs != "EPSG:4326":
+        raise ValueError("Make sure to pass data with EPSG:4326 projection")
+
+    long, lat = measurement_coords.x, measurement_coords.y
+    bounds = (lat.min(), long.min(), lat.max(), long.max())
+
+    result = make_overpy_request(
+        f"""
+        [out:json][timeout:25];
+        (
+        way['building']{bounds};
+        relation["building"]{bounds};
+        );
+        out body;
+        >;
+        out skel qt;
+        """,
+        retries
+    )
+
+    blockages = gpd.GeoSeries(
+        Polygon(
+            ((node.lon, node.lat) for node in way.nodes)
+        ) for way in result.ways
+    ).set_crs("EPSG:4326")
+
+    return blockages
+
+def project_onto_streets(point_series, street_series, epsg="EPSG:31287", plot=False):
+    '''
+    Todo: Muss noch angepasst werden!
+    '''
+
+    if point_series.crs != epsg or street_series.crs != epsg:
+        raise ValueError("GeoSeries does not match value set by epsg argument (Default EPSG:31287)")
+
+    projected = []
+    for point in point_series:
+        street_ind = street_series.distance(point).argmin()
+        projected.append(
+            street_series.loc[street_ind].interpolate(street_series.loc[street_ind].project(point))
+        )
+
+    projected = gpd.GeoSeries(projected)
+    projected = projected.set_crs(epsg)
+    deviation_projection = projected.distance(point_series)
+
+    if plot:
+        plt.figure()
+        plt.title("Deviation Histogram")
+        deviation_projection.plot.hist(bins=50)
+        plt.show()
+
+        fig, (ax1, ax2) = plt.subplots(1,2, figsize=(20, 10))
+        ax1.set_title("Raw Measurements")
+        ax2.set_title("After Projection")
+        ax1.scatter(point_series.x, point_series.y, label="Point Series")
+        clb = ax2.scatter(projected.x, projected.y, c=deviation_projection, label="Projected Series")
+        for street in street_series:
+            ax1.plot(*street.xy, color="black")
+            ax2.plot(*street.xy, color="black")
+        bar = fig.colorbar(clb)
+        bar.ax.set_title('Deviation [Meters]')
+        ax1.legend()
+        ax2.legend()
+        plt.show()
 
-class OSM_API():
-    #Collect all the possible queries here 
-    pass
+    return projected, deviation_projection

tests/geospatial_test.py

 from context import measprocess as mpc
+import unittest
+import geopandas as gpd
+from shapely.geometry import Point, LineString
 
 class TestOSMAPI(unittest.TestCase):
-    pass
+    def setUp(self):
+        self._coords = [
+            (48.201914, 16.363859),
+            (48.194170, 16.385466)
+        ]
 
-class TestProjections(unittest.TestCase):
-    pass
+        self._measurement_coords_not_set = gpd.GeoSeries(
+            (Point(lon, lat) for lat, lon in self._coords)
+        )
+
+        self._measurement_coords_set = gpd.GeoSeries(
+            (Point(lon, lat) for lat, lon in self._coords)
+        ).set_crs("EPSG:4326")
+
+    def test_basic_API(self):
+        street_series = mpc.geospatial.get_geoseries_streets(
+            self._measurement_coords_set
+        )
+
+        blockage_series = mpc.geospatial.get_geoseries_blockages(
+            self._measurement_coords_set
+        )
+
+        self.assertTrue(
+            "EPSG:4326" == street_series.crs
+        )
+
+        self.assertTrue(
+            "EPSG:4326" == blockage_series.crs
+        )
+
+    def test_raise_error_if_crs_not_set(self):
+
+        with self.assertRaises(ValueError):
+            mpc.geospatial.get_geoseries_blockages(
+                self._measurement_coords_not_set
+            )
+
+            mpc.geospatial.get_geoseries_streets(
+                self._measurement_coords_not_set
+            )
+
+    def test_output_coordinates_streets(self):
+
+        street_series = mpc.geospatial.get_geoseries_streets(
+            self._measurement_coords_set
+        )
+
+        street_series = street_series.to_crs("EPSG:31287")
+        projected = self._measurement_coords_set.to_crs("EPSG:31287")
+
+        self.assertTrue(
+            street_series.distance(projected[0]).mean() < 1000
+        )
+
+    def test_output_coordinates_blockages(self):
+
+        blockage_series = mpc.geospatial.get_geoseries_blockages(
+            self._measurement_coords_set
+        )
+
+        blockage_series = blockage_series.to_crs("EPSG:31287")
+        projected = self._measurement_coords_set.to_crs("EPSG:31287")
+
+        self.assertTrue(
+            blockage_series.distance(projected[0]).mean() < 1000
+        )

tests/projection_test.py

+import unittest
+import geopandas as gpd
+from shapely.geometry import Point, LineString
+from context import measprocess as mpc
+import numpy as np
+
+class TestProjections(unittest.TestCase):
+    def setUp(self):
+        self._street_series = gpd.GeoSeries([
+                LineString([
+                    (0, 0),
+                    (10, 0),
+                    (10, 10)
+                ]),
+                LineString([
+                    (100, 100),
+                    (200, 200)
+                ])
+        ])
+
+        self._point_series = gpd.GeoSeries([
+            Point(0, -10), #Should be mapped to 0, 0
+            Point(10, -1), #Should be mapped to 10, 0
+            Point(20, 20) #Should be mapped to 10, 10
+        ])
+
+
+
+    def test_other_projection(self):
+        mpc.geospatial.project_onto_streets(
+            self._point_series.set_crs("EPSG:3416"),
+            self._street_series.set_crs("EPSG:3416"),
+            epsg="EPSG:3416"
+        )
+
+    def test_exception_epsg(self):
+        with self.assertRaises(ValueError):
+            mpc.geospatial.project_onto_streets(
+                self._point_series,
+                self._street_series,
+                epsg="EPSG:31287"
+            )
+
+            mpc.geospatial.project_onto_streets(
+                self._point_series.set_crs("EPSG:31287"),
+                self._street_series.set_crs("EPSG:31287"),
+                epsg="EPSG:3416"
+            )
+
+    def test_basic_projection(self):
+        projected, deviation = mpc.geospatial.project_onto_streets(
+            self._point_series.set_crs("EPSG:31287"),
+            self._street_series.set_crs("EPSG:31287")
+        )
+
+        distances = projected.distance(
+            gpd.GeoSeries((
+                Point(x, y) for x, y in [
+                    (0, 0),
+                    (10, 0),
+                    (10, 10)
+                ]
+            )).set_crs("EPSG:31287")
+        )
+
+        self.assertTrue(distances.max() < 1e-1)
+
+        self.assertTrue(
+            np.max(
+                deviation.values - np.array([10, 1, np.sqrt(10**2+10**2)])
+            ) < 1e-1
+        )