Abfrage und Visualisierung von OSM-Daten mit Python

/*
Table Tennis Tables
*/
[out:json];
node["sport"="table_tennis"]({{bbox}});
out body; >; out skel qt;

import plotly.express as px
import geopandas as gpd
import matplotlib.pyplot as plt
from collections import namedtuple
from osm2geojson import overpass_call, json2geojson

# BBox mit den Komponenten in der Reihenfolge, die von overpass.ql_query erwartet wird
Bbox = namedtuple("Bbox", ["south", "west", "north", "east"])
Tag = namedtuple("Tag", ["key", "values"])

def load_osm_from_overpass(bbox, tag, crs="epsg:4326") -> gpd.GeoDataFrame:
    geojson = load_osm_from_overpass_geojson(bbox, tag, crs)
    return gpd.GeoDataFrame.from_features(geojson, crs=crs)

def load_osm_from_overpass_geojson(bbox, tag, crs="epsg:4326"):
    values = None if "*" in tag.values else tag.values
    query = f"""
    [out:json];
    node["{tag.key}"="{values[0] if values else '*'}"]({bbox.south},{bbox.west},{bbox.north},{bbox.east});
    out body;
    >;
    out skel qt;
    """
    response = overpass_call(query)
    return json2geojson(response)

def plot_gdf_on_map(gdf, color="k", title=""):
    fig = px.scatter_mapbox(
        gdf,
        lat=gdf.geometry.y,
        lon=gdf.geometry.x,
        color_discrete_sequence=["fuchsia"],
        height=500,
        zoom=11,
    )
    fig.update_layout(title=title, mapbox_style="open-street-map")
    fig.update_layout(margin={"r": 0, "l": 0, "b": 0})
    return fig

def main():
    # Beispiel: Tischtennistische in Zürich
    tag = Tag(key="sport", values=["table_tennis"])
    bbox_zurich = Bbox(west=8.471668, south=47.348834, east=8.600454, north=47.434379)
    table_tennis_gdf = load_osm_from_overpass(bbox_zurich, tag)
    gdf = table_tennis_gdf.to_crs(epsg=3857)
    gdf.geometry = gdf.geometry.centroid
    gdf = gdf.to_crs(epsg=4236)
    fig = plot_gdf_on_map(gdf, title="Zurich Table Tennis")
    # Falls kein Fenster erscheint, die folgende Zeile auskommentieren und die Datei standalone_map.html im Browser öffnen
    # fig.write_html("standalone_map.html")
    fig.show()

if __name__ == "__main__":
    main()

CREATE LANGUAGE plpythonu;

CREATE OR REPLACE FUNCTION geocode(address text)
RETURNS text AS $$
"""Returns WKT geometry which is EMPTY if address was not found.

Please respect the Nominatim free usage policy:
https://operations.osmfoundation.org/policies/nominatim/
"""

import requests
from time import sleep

base_url = 'https://nominatim.openstreetmap.org/search?'
params = {
    'q': address,
    'format': 'json',
    'limit': 1,
    'User-Agent': 'geocode hack (pls. replace)'
}
response = requests.get(base_url, params=params)
loc = response.json()
if len(loc) == 0:
    wkt = 'SRID=4326;POINT EMPTY'
else:
    wkt = f'SRID=4326;POINT({loc[0]['lon']} {loc[0]['lat']})'
sleep(1)  # Throttle, to avoid overloading Nominatim
return wkt
$$ LANGUAGE plpythonu VOLATILE COST 1000;

-- Test:
SELECT ST_AsGeoJSON(geocode('Vulkanstrasse 106, Zürich'));

WITH my_customers_tmp(id, address, geom) AS (
  SELECT id, address, ST_Transform(ST_MakeValid(GeomFromEWKT(geocode(address))),3857) AS geom
  FROM (values
   (1,'Vulkanstrasse 106, Zürich'),
   (2,'Bellevue, 8001 Zürich'),
   (3,'Weihnachtsmann')
  ) address(id,address)
),
poi_tmp AS (
  SELECT
        osm_id,
        gtype,
        -- Die Adresse des POI müsste rückgeokodiert werden!
        tags,
        way AS geom FROM osm_poi
  WHERE tags->'sport' IN ('table_tennis')
)
SELECT
  my_customers_tmp.id as customer_id,
  ST_Distance(my_customers_tmp.geom, poi_tmp.geom)::integer AS distance,
  poi_tmp.osm_id,
  poi_tmp.gtype,
  ST_Transform(poi_tmp.geom,4326) AS geom
FROM my_customers_tmp, poi_tmp
ORDER BY my_customers_tmp.geom <-> poi_tmp.geom
LIMIT 5;