33 lines
1.1 KiB
Python
33 lines
1.1 KiB
Python
import geopandas as gpd
|
|
from scipy.spatial import cKDTree
|
|
import numpy as np
|
|
import pyproj
|
|
from functools import lru_cache
|
|
|
|
# Constants
|
|
EPSG_WGS84 = 4326
|
|
EPSG_WebMercator = 3857
|
|
|
|
# Load the geodata.
|
|
gdf = gpd.read_file("berlin-latest-free.shp/gis_osm_pois_free_1.shp")
|
|
|
|
# Get all restaurant points and represent as meters on web mercator projection.
|
|
rest = gdf[gdf['fclass'].isin(['restaurant','cafe','fast_food','biergarten','pub'])].to_crs(epsg=EPSG_WGS84)
|
|
|
|
# Construct a cKDTree of all restaurant points for easy nearest-neighbor lookup.
|
|
coords_rest = np.vstack([rest.geometry.x, rest.geometry.y]).T
|
|
tree_rest = cKDTree(coords_rest)
|
|
|
|
# Construct a projection from lat/long to web mercator.
|
|
proj_rest = pyproj.Transformer.from_crs(EPSG_WGS84, EPSG_WebMercator, always_xy=True)
|
|
|
|
def distance_to_nearest_restaurant(lon: float, lat: float) -> float:
|
|
|
|
# Project coordinates to mercator.
|
|
x, y = proj.transform(lon, lat)
|
|
|
|
# Get the closest point from the tree.
|
|
dist, idx = tree.query([x, y])
|
|
|
|
# Returned distance is in meters (from mercator projection).
|
|
return float(dist) |