diff --git a/measprocess/geospatial.py b/measprocess/geospatial.py index d426cc683e392af6809242dac002799785f0624c..4a3b30b3a1fc0aee42883f576e1ca53d790b5aae 100644 --- a/measprocess/geospatial.py +++ b/measprocess/geospatial.py @@ -1,10 +1,11 @@ import overpy import geopandas as gpd -from shapely.geometry import LineString, Polygon +from shapely.geometry import LineString, Polygon, box from scipy.spatial.distance import cdist import numpy as np import warnings import time +import matplotlib.pyplot as plt def make_overpy_request(request_body : str, retries : int): for _ in range(retries): @@ -19,11 +20,13 @@ def make_overpy_request(request_body : str, retries : int): else: return result -def get_geoseries_streets(measurement_coords : gpd.GeoSeries, retries : int = 5) -> gpd.GeoSeries: +def get_geoseries_streets(measurement_coords : gpd.GeoSeries, retries : int = 5, crop : bool = False) -> gpd.GeoSeries: ''' 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 + :param retries: Number of retries for overpy requests, default:5 + :param crop: Set to True if the returned geoseries should be fixed to bounding box determined by measurement_coords, default: False :return: A geopandas geoseries consiting of LineStrings representing the streets in EPSG:4326 ''' @@ -49,10 +52,13 @@ def get_geoseries_streets(measurement_coords : gpd.GeoSeries, retries : int = 5) street_series = gpd.GeoSeries( LineString( - ((node.lon, node.lat) for node in way.nodes) + ((node.lon, node.lat) for node in way.nodes if len(way.nodes)>1) ) for way in result.ways ).set_crs("EPSG:4326") + if crop == True: + street_series = crop_geoseries_by_bounds(street_series, measurement_coords) + return street_series def get_geoseries_blockages(measurement_coords : gpd.GeoSeries, retries : int = 5) -> gpd.GeoSeries: @@ -132,3 +138,22 @@ def project_onto_streets(point_series : gpd.GeoSeries, street_series : gpd.GeoSe plt.show() return projected, deviation_projection + +def crop_geoseries_by_bounds(original_geoseries : gpd.GeoSeries, filter_geoseries : gpd.GeoSeries) -> gpd.GeoSeries: + """ + Returns cut off geoseries from 'original_geoseries' which lies within the bounds on the filter_geoseries + + :param original_geoseries: A geopandas geoseries of points/streets/polygons to be filtered, represented in a EPSG:4326 projection + :param filter_geoseries: A geopandas geoseries based on which a bounding box is for filter is made + + :return: geopandas geoseries of same format as original_geoseries + """ + if original_geoseries.crs != "EPSG:4326" or filter_geoseries.crs != "EPSG:4326" : + raise ValueError("Make sure to pass data with EPSG:4326 projection") + + filter_bounds = filter_geoseries.total_bounds + bound_box = box(filter_bounds[0], filter_bounds[1],filter_bounds[2],filter_bounds[3]) + bound_box = gpd.GeoDataFrame(geometry = gpd.GeoSeries(bound_box)).set_crs("EPSG:4326") # create geodataframe for overlay + filtered_series = (gpd.overlay(gpd.GeoDataFrame(geometry = original_geoseries),bound_box, how='intersection')).geometry + + return filtered_series \ No newline at end of file