import time import warnings import geopandas as gpd import overpy from shapely.geometry import LineString, Polygon import os import shutil import urllib.request import zipfile from typing import List import pandas as pd import tempfile def make_overpy_request(request_body: str, retries: int): 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 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 """ if measurement_coords.crs != "EPSG:4326": raise ValueError("Make sure to pass data with EPSG:4326 projection") # Find total bounds: (lat_min, long_min, lat_max, long_max) bounds = tuple(measurement_coords.total_bounds[[1, 0, 3, 2]]) 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 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) street_series = ( street_series.explode() ) # if any MultiLineStrings, they are transformed to LineStrings return street_series def geoseries_buildings( measurement_coords: gpd.GeoSeries, retries: int = 5 ) -> gpd.GeoSeries: """ 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") # Find total bounds: (lat_min, long_min, lat_max, long_max) bounds = tuple(measurement_coords.total_bounds[[1, 0, 3, 2]]) 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 if len(way.nodes) > 2 ) # Polygon needs at least 3 nodes ) for way in result.ways ).set_crs("EPSG:4326") return blockages def pull_from_geodatenviewer_list(squares: List[str]) -> gpd.GeoDataFrame: """ Downloads .zip files for specified squares from https://www.wien.gv.at/ma41datenviewer/public/, extracts and combines all .shp files to single geopandas geodataframe, then deletes the dir_path directory. :param dir_path: path to the directory where the files will be downloaded and extracted to :param squares: list of integers denoting the squares to be downloaded :return: geopandas geodataframe containing all data from .shp files """ dir_path = tempfile.mkdtemp() # create temporary directory # downloading .zip files # example of single square with data: https://www.wien.gv.at/ma41datenviewer/downloads/ma41/geodaten/fmzk_bkm/103081_bkm.zip for square_number in squares: try: urllib.request.urlretrieve( "https://www.wien.gv.at/ma41datenviewer/downloads/geodaten/fmzk_bkm/{}_bkm.zip".format( square_number ), os.path.join(dir_path, "{}.zip".format(square_number)), ) except: continue # extracting and deleting .zip files for item in os.listdir(dir_path): # loop through items in dir if item.endswith(".zip"): file_name = os.path.join(dir_path, item) zip_ref = zipfile.ZipFile(file_name) # create zipfile object zip_ref.extractall(dir_path) # extract file to dir zip_ref.close() os.remove(file_name) # delete zipped file # combine all .shp files geo_df_all = pd.DataFrame() for square in squares: try: geo_df = gpd.read_file(os.path.join(dir_path, square + "_bkm.shp")) geo_df_all = pd.concat([geo_df_all, geo_df], ignore_index=True) except: continue shutil.rmtree(dir_path) # deletes the directory containing all the files return geo_df_all def pull_from_geodatenviewer_meas( measurement_coords: gpd.GeoSeries, ) -> gpd.GeoDataFrame: """ Downloads raster .zip file for Vienna Austria, extracts all raster polygon IDs containing measurements, then calls pull_from_geodatenviewer_list() to extract building polygons. :param measurement_coords: geopandas geoseries containing measurements in EPSG:4326 projection :return: geopandas geodataframe containing all data from .shp files """ if measurement_coords.crs != "EPSG:4326": raise ValueError("Make sure to pass data with EPSG:4326 projection") raster_path = tempfile.mkdtemp() urllib.request.urlretrieve( "https://data.wien.gv.at/daten/geo?service=WFS&request=GetFeature&version=1.1.0&typeName=ogdwien:MZKBLATT1000OGD&srsName=EPSG:4326&outputFormat=shape-zip", os.path.join(raster_path, "raster.zip"), ) file_name = os.path.join(raster_path, "raster.zip") zip_ref = zipfile.ZipFile(file_name) # create zipfile object zip_ref.extractall(raster_path) # extract file to dir zip_ref.close() os.remove(file_name) # delete zipped file for file in os.listdir(raster_path): if file.endswith(".shp"): raster_gdf = gpd.read_file(os.path.join(raster_path, file)) raster_gdf.to_crs(4326,inplace=True) meas_gdf = gpd.GeoDataFrame(geometry=measurement_coords) poly_gdf = gpd.sjoin( raster_gdf, meas_gdf, op="contains" ) # find all polygons from raster that contain measurements polygonID_list = poly_gdf.MZK1000.unique() return_gdf = pull_from_geodatenviewer_list(polygonID_list) return return_gdf