266 lines
8.0 KiB
Python
266 lines
8.0 KiB
Python
from shapely import geometry, ops, affinity
|
|
|
|
import axi
|
|
import math
|
|
import osm2shapely
|
|
import sys
|
|
|
|
import logging
|
|
logging.basicConfig()
|
|
|
|
HOME = 35.768616, -78.844005
|
|
RALEIGH = 35.777486, -78.635794
|
|
DURHAM = 35.999392, -78.919217
|
|
PARIS = 48.856744, 2.351248
|
|
HOLBROOK = 32.422131, -80.713289
|
|
FORMLABS = 42.374414, -71.087908
|
|
PINE_RIDGE_ROAD = 42.427164, -71.143553
|
|
TIMES_SQUARE = 40.758582, -73.985066
|
|
|
|
LAT, LNG = TIMES_SQUARE
|
|
LANDSCAPE = False
|
|
PAGE_WIDTH_IN = 12
|
|
PAGE_HEIGHT_IN = 8.5
|
|
MAP_WIDTH_KM = 1.61
|
|
LANE_WIDTH_M = 3.7
|
|
EARTH_RADIUS_KM = 6371
|
|
|
|
WEIGHTS = {
|
|
'motorway': 2,
|
|
'motorway_link': 2,
|
|
'trunk_link': 2,
|
|
'trunk': 2,
|
|
'primary_link': 1.75,
|
|
'primary': 1.75,
|
|
'secondary': 1.5,
|
|
'secondary_link': 1.5,
|
|
'tertiary_link': 1.25,
|
|
'tertiary': 1.25,
|
|
'living_street': 1,
|
|
'unclassified': 1,
|
|
'residential': 1,
|
|
# 'service': 0,
|
|
# 'railway': 0,
|
|
# 'pedestrian': 0,
|
|
# 'footway': 0,
|
|
# 'natural': 0,
|
|
}
|
|
|
|
def paths_to_shapely(paths):
|
|
return geometry.MultiLineString(paths)
|
|
|
|
def shapely_to_paths(g):
|
|
if isinstance(g, geometry.Point):
|
|
return []
|
|
elif isinstance(g, geometry.LineString):
|
|
return [list(g.coords)]
|
|
elif isinstance(g, (geometry.MultiPoint, geometry.MultiLineString, geometry.MultiPolygon, geometry.collection.GeometryCollection)):
|
|
paths = []
|
|
for x in g:
|
|
paths.extend(shapely_to_paths(x))
|
|
return paths
|
|
elif isinstance(g, geometry.Polygon):
|
|
paths = []
|
|
paths.append(list(g.exterior.coords))
|
|
for interior in g.interiors:
|
|
paths.extend(shapely_to_paths(interior))
|
|
return paths
|
|
else:
|
|
raise Exception('unhandled shapely geometry: %s' % type(g))
|
|
|
|
def follow(g, step):
|
|
result = []
|
|
d = step
|
|
e = 1e-3
|
|
l = g.length
|
|
while d < l:
|
|
p = g.interpolate(d)
|
|
a = g.interpolate(d - e)
|
|
b = g.interpolate(d + e)
|
|
angle = math.atan2(b.y - a.y, b.x - a.x)
|
|
result.append((p.x, p.y, angle))
|
|
d += step
|
|
return result
|
|
|
|
def hatch(g, angle, step):
|
|
print g.area
|
|
x0, y0, x1, y1 = g.bounds
|
|
d = max(x1 - x0, y1 - y0) * 2
|
|
lines = []
|
|
x = 0
|
|
while x < d:
|
|
lines.append(geometry.LineString([(x - d / 2, -d / 2), (x - d / 2, d / 2)]))
|
|
x += step
|
|
if not lines:
|
|
return None
|
|
lines = geometry.collection.GeometryCollection(lines)
|
|
lines = affinity.rotate(lines, angle)
|
|
lines = affinity.translate(lines, (x0 + x1) / 2, (y0 + y1) / 2)
|
|
return g.intersection(lines)
|
|
|
|
def crop(g, w, h):
|
|
return g.intersection(geometry.Polygon(box(w, h)))
|
|
|
|
def box(w, h):
|
|
w *= 0.5
|
|
h *= 0.5
|
|
return [(-w, -h), (w, -h), (w, h), (-w, h), (-w, -h)]
|
|
|
|
def haversine(lat1, lng1, lat2, lng2):
|
|
lng1, lat1, lng2, lat2 = map(math.radians, [lng1, lat1, lng2, lat2])
|
|
dlng = lng2 - lng1
|
|
dlat = lat2 - lat1
|
|
a = math.sin(dlat/2)**2 + math.cos(lat1) * math.cos(lat2) * math.sin(dlng/2)**2
|
|
return math.asin(math.sqrt(a)) * 2 * EARTH_RADIUS_KM
|
|
|
|
class LambertAzimuthalEqualAreaProjection(object):
|
|
def __init__(self, lat, lng):
|
|
self.lat = lat
|
|
self.lng = lng
|
|
self.angle = math.radians(29)
|
|
self.scale = 1
|
|
self.scale = self.kilometer_scale()
|
|
def project(self, lng, lat):
|
|
lng, lat = math.radians(lng), math.radians(lat)
|
|
clng, clat = math.radians(self.lng), math.radians(self.lat)
|
|
k = math.sqrt(2 / (1 + math.sin(clat)*math.sin(lat) + math.cos(clat)*math.cos(lat)*math.cos(lng-clng)))
|
|
x = k * math.cos(lat) * math.sin(lng-clng)
|
|
y = k * (math.cos(clat)*math.sin(lat) - math.sin(clat)*math.cos(lat)*math.cos(lng-clng))
|
|
rx = x * math.cos(self.angle) - y * math.sin(self.angle)
|
|
ry = y * math.cos(self.angle) + x * math.sin(self.angle)
|
|
x = rx
|
|
y = ry
|
|
s = self.scale
|
|
return (x * s, -y * s)
|
|
def kilometer_scale(self):
|
|
e = 1e-3
|
|
lat, lng = self.lat, self.lng
|
|
km_per_lat = haversine(lat - e, lng, lat + e, lng) / (2 * e)
|
|
km_per_lng = haversine(lat, lng - e, lat, lng + e) / (2 * e)
|
|
x1, y1 = self.project(lng - 1 / km_per_lng, lat - 1 / km_per_lat)
|
|
x2, y2 = self.project(lng + 1 / km_per_lng, lat + 1 / km_per_lat)
|
|
sx = 2 / (x2 - x1)
|
|
sy = 2 / (y1 - y2)
|
|
return (sx + sy) / 2
|
|
def transform(self, g):
|
|
result = ops.transform(self.project, g)
|
|
result.tags = g.tags
|
|
return result
|
|
|
|
def circle(cx, cy, r, n, revs=5):
|
|
points = []
|
|
n *= revs
|
|
for i in range(n + 1):
|
|
a = revs * math.pi * i / n
|
|
x = cx + math.cos(a) * r
|
|
y = cy + math.sin(a) * r
|
|
points.append((x, y))
|
|
return points
|
|
|
|
def main():
|
|
# parse osm file into shapely geometries
|
|
geometries = osm2shapely.parse(sys.argv[1])
|
|
|
|
# setup map projection
|
|
projection = LambertAzimuthalEqualAreaProjection(LAT, LNG)
|
|
geometries = [projection.transform(g) for g in geometries]
|
|
|
|
# determine width and height of map
|
|
w = MAP_WIDTH_KM
|
|
if LANDSCAPE:
|
|
h = float(w) * PAGE_HEIGHT_IN / PAGE_WIDTH_IN
|
|
else:
|
|
h = float(w) * PAGE_WIDTH_IN / PAGE_HEIGHT_IN
|
|
|
|
# process geometries
|
|
gs = []
|
|
roads = []
|
|
for g in geometries:
|
|
highway = g.tags.get('highway')
|
|
if highway and highway in WEIGHTS:
|
|
weight = WEIGHTS[highway]
|
|
if weight:
|
|
g = g.buffer(LANE_WIDTH_M / 1000.0 * weight)
|
|
g = crop(g, w * 1.1, h * 1.1)
|
|
roads.append(g)
|
|
# elif 'natural' in g.tags:
|
|
# gs.append(g)
|
|
elif 'building' in g.tags:
|
|
gs.append(g)
|
|
|
|
print 'crop'
|
|
gs = [crop(g, w * 1.1, h * 1.1) for g in gs]
|
|
roads = [crop(g, w * 1.1, h * 1.1) for g in roads]
|
|
|
|
# gs = []
|
|
# for key, ways in handler.ways.items():
|
|
# if key not in WEIGHTS:
|
|
# print 'skip', key, len(ways)
|
|
# continue
|
|
# print 'layer', key, len(ways)
|
|
# ggs = []
|
|
# for way in ways:
|
|
# coords = [projection.project(*handler.coords[x]) for x in way]
|
|
# if key == 'natural' and coords[0] == coords[-1]:
|
|
# ggs.append(geometry.Polygon(coords))
|
|
# else:
|
|
# ggs.append(geometry.LineString(coords))
|
|
# # g = paths_to_shapely(paths)
|
|
# g = geometry.collection.GeometryCollection(ggs)
|
|
# if key == 'railway':
|
|
# paths = shapely_to_paths(g)
|
|
# g = paths_to_shapely(paths)
|
|
# points = follow(g, 20 / 1000.0)
|
|
# s = 4 / 1000.0
|
|
# for x, y, a in points:
|
|
# x1 = x + math.cos(a + math.pi / 2) * s
|
|
# y1 = y + math.sin(a + math.pi / 2) * s
|
|
# x2 = x + math.cos(a - math.pi / 2) * s
|
|
# y2 = y + math.sin(a - math.pi / 2) * s
|
|
# paths.append([(x1, y1), (x2, y2)])
|
|
# g = paths_to_shapely(paths)
|
|
# if key == 'natural':
|
|
# gs.append(crop(paths_to_shapely(shapely_to_paths(g)), w * 1.1, h * 1.1))
|
|
# paths = hatch(g, 45, 10 / 1000.0)
|
|
# g = paths_to_shapely(paths)
|
|
# weight = WEIGHTS[key]
|
|
# if weight:
|
|
# g = g.buffer(LANE_WIDTH_M / 1000.0 * weight)
|
|
# g = crop(g, w * 1.1, h * 1.1)
|
|
# gs.append(g)
|
|
|
|
print 'union'
|
|
roads = ops.cascaded_union(roads)
|
|
all_roads = []
|
|
while not roads.is_empty:
|
|
all_roads.append(roads)
|
|
roads = roads.buffer(-LANE_WIDTH_M / 1000.0 / 3)
|
|
g = geometry.collection.GeometryCollection(gs + all_roads)
|
|
g = paths_to_shapely(shapely_to_paths(g))
|
|
|
|
print 'crop'
|
|
g = crop(g, w, h)
|
|
|
|
paths = shapely_to_paths(g)
|
|
|
|
# dot at origin
|
|
# s = 3 / 1000.0 / 3
|
|
# for i in range(4):
|
|
# paths.append(circle(0, 0, i * s, 360))
|
|
|
|
# border around map
|
|
s = 6 / 1000.0 / 2
|
|
for m in range(1):
|
|
paths.append(box(w - s * m, h - s * m))
|
|
|
|
print 'axi'
|
|
d = axi.Drawing(paths)
|
|
d = d.rotate_and_scale_to_fit(PAGE_WIDTH_IN, PAGE_HEIGHT_IN, step=90)
|
|
d = d.sort_paths().join_paths(0.002).simplify_paths(0.002)
|
|
im = d.render()
|
|
im.write_to_png('out.png')
|
|
axi.draw(d)
|
|
|
|
if __name__ == '__main__':
|
|
main()
|