kopia lustrzana https://github.com/marceloprates/prettymaps
106 wiersze
3.9 KiB
Python
106 wiersze
3.9 KiB
Python
# OpenStreetMap Networkx library to download data from OpenStretMap
|
|
import osmnx as ox
|
|
|
|
# CV2 & Scipy & Numpy & Pandas
|
|
import numpy as np
|
|
|
|
# Shapely
|
|
from shapely.geometry import *
|
|
from shapely.affinity import *
|
|
|
|
# Geopandas
|
|
from geopandas import GeoDataFrame
|
|
|
|
# Matplotlib
|
|
from matplotlib.path import Path
|
|
|
|
# etc
|
|
from collections.abc import Iterable
|
|
from functools import reduce
|
|
|
|
# Helper functions to fetch data from OSM
|
|
|
|
def ring_coding(ob):
|
|
codes = np.ones(len(ob.coords), dtype = Path.code_type) * Path.LINETO
|
|
codes[0] = Path.MOVETO
|
|
return codes
|
|
|
|
def pathify(polygon):
|
|
vertices = np.concatenate([np.asarray(polygon.exterior)] + [np.asarray(r) for r in polygon.interiors])
|
|
codes = np.concatenate([ring_coding(polygon.exterior)] + [ring_coding(r) for r in polygon.interiors])
|
|
return Path(vertices, codes)
|
|
|
|
def union(geometry):
|
|
geometry = np.concatenate([[x] if type(x) == Polygon else x for x in geometry if type(x) in [Polygon, MultiPolygon]])
|
|
geometry = reduce(lambda x, y: x.union(y), geometry[1:], geometry[0])
|
|
return geometry
|
|
|
|
def get_perimeter(query, by_osmid = False):
|
|
return ox.geocode_to_gdf(query, by_osmid = by_osmid)
|
|
|
|
def get_footprints(perimeter = None, point = None, radius = None, footprint = 'building'):
|
|
|
|
if perimeter is not None:
|
|
# Boundary defined by polygon (perimeter)
|
|
footprints = ox.geometries_from_polygon(union(perimeter.geometry), tags = {footprint: True} if type(footprint) == str else footprint)
|
|
perimeter = union(ox.project_gdf(perimeter).geometry)
|
|
|
|
elif (point is not None) and (radius is not None):
|
|
# Boundary defined by circle with radius 'radius' around point
|
|
footprints = ox.geometries_from_point(point, dist = radius, tags = {footprint: True} if type(footprint) == str else footprint)
|
|
perimeter = GeoDataFrame(geometry=[Point(point[::-1])], crs = footprints.crs)
|
|
perimeter = ox.project_gdf(perimeter).geometry[0].buffer(radius)
|
|
|
|
if len(footprints) > 0:
|
|
footprints = ox.project_gdf(footprints)
|
|
|
|
footprints = [
|
|
[x] if type(x) == Polygon else x
|
|
for x in footprints.geometry if type(x) in [Polygon, MultiPolygon]
|
|
]
|
|
footprints = list(np.concatenate(footprints)) if len(footprints) > 0 else []
|
|
footprints = [pathify(x) for x in footprints if x.within(perimeter)]
|
|
|
|
return footprints, perimeter
|
|
|
|
def get_streets(perimeter = None, point = None, radius = None, dilate = 6, custom_filter = None):
|
|
|
|
if perimeter is not None:
|
|
# Boundary defined by polygon (perimeter)
|
|
streets = ox.graph_from_polygon(union(perimeter.geometry), custom_filter = custom_filter)
|
|
streets = ox.project_graph(streets)
|
|
streets = ox.graph_to_gdfs(streets, nodes = False)
|
|
#streets = ox.project_gdf(streets)
|
|
streets = MultiLineString(list(streets.geometry)).buffer(dilate)
|
|
|
|
elif (point is not None) and (radius is not None):
|
|
# Boundary defined by polygon (perimeter)
|
|
|
|
streets = ox.graph_from_point(point, dist = radius, custom_filter = custom_filter)
|
|
crs = ox.graph_to_gdfs(streets, nodes = False).crs
|
|
streets = ox.project_graph(streets)
|
|
|
|
perimeter = GeoDataFrame(geometry=[Point(point[::-1])], crs = crs)
|
|
perimeter = ox.project_gdf(perimeter).geometry[0].buffer(radius)
|
|
|
|
streets = ox.graph_to_gdfs(streets, nodes = False)
|
|
|
|
streets = MultiLineString(list(
|
|
filter(
|
|
# Filter lines with at least 2 points
|
|
lambda line: len(line) >= 2,
|
|
# Iterate over lines in geometry
|
|
map(
|
|
# Filter points within perimeter
|
|
lambda line: list(filter(lambda xy: Point(xy).within(perimeter), zip(*line.xy))),
|
|
streets.geometry
|
|
)
|
|
)
|
|
)).buffer(dilate) # Dilate lines
|
|
|
|
if not isinstance(streets, Iterable):
|
|
streets = [streets]
|
|
|
|
streets = list(map(pathify, streets))
|
|
|
|
return streets, perimeter |