fetch.py

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