diff --git a/README.md b/README.md index c1f4eb76d6cfe9d958722861d8fafaabc406548c..0afcc7212d3f6c180d82aab449052b5d3d63ad79 100644 --- a/README.md +++ b/README.md @@ -12,7 +12,7 @@ Either clone the repo directly or use pip install. Make sure to replace `` 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` diff --git a/measprocess/__init__.py b/measprocess/__init__.py index 4c390daf09603a9e85d6df34402a605cf6b90d3c..72b757876214c35fe7cbf18717dcc6ad7613be7a 100644 --- a/measprocess/__init__.py +++ b/measprocess/__init__.py @@ -1 +1,3 @@ from .preprocess import link_dataframes +from .geospatial import get_geoseries_streets, get_geoseries_blockages +from .geospatial import project_onto_streets diff --git a/measprocess/geospatial.py b/measprocess/geospatial.py index 2216ef267facddf1a0a40fb58f60ccc87c90299b..46260d1565644f85398001bdd0a13d28bc200d05 100644 --- a/measprocess/geospatial.py +++ b/measprocess/geospatial.py @@ -1,7 +1,135 @@ -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 diff --git a/tests/geospatial_test.py b/tests/geospatial_test.py index fb1b6f20cb6ee8198f801a50fcff349ecc96553b..9e815112be5b13c8116337bbff04f5cae7c3ee65 100644 --- a/tests/geospatial_test.py +++ b/tests/geospatial_test.py @@ -1,7 +1,73 @@ 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 + ) diff --git a/tests/projection_test.py b/tests/projection_test.py new file mode 100644 index 0000000000000000000000000000000000000000..3ad60771492f08008e73355feb1c059e5a335c2e --- /dev/null +++ b/tests/projection_test.py @@ -0,0 +1,72 @@ +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 + )