Created
May 28, 2019 07:25
-
-
Save mlankenau/cf6a7528020c77deb26899eb49544f95 to your computer and use it in GitHub Desktop.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
from PIL import Image, ImageDraw | |
# https://github.com/dezhin/osmread | |
from osmread import parse_file, Way, Node, Relation | |
points = {} | |
ways = {} | |
relations = [] | |
lat_from = 47.7192 | |
lon_from = 12.4236 | |
lat_to = 47.7525 | |
lon_to = 12.5034 | |
texture_size = 1024 | |
texture_border = 100 | |
def geo_to_img(map): | |
lat = map["lat"] | |
lon = map["lon"] | |
y = texture_size - (lat - lat_from) / (lat_to - lat_from) * texture_size | |
x = (lon - lon_from) / (lon_to - lon_from) * texture_size | |
return (int(texture_border+x), int(texture_border+y)) | |
# source: https://rosettacode.org/wiki/Sutherland-Hodgman_polygon_clipping#Python | |
def clip(subjectPolygon, clipPolygon): | |
def inside(p): | |
return(cp2[0]-cp1[0])*(p[1]-cp1[1]) > (cp2[1]-cp1[1])*(p[0]-cp1[0]) | |
def computeIntersection(): | |
dc = [ cp1[0] - cp2[0], cp1[1] - cp2[1] ] | |
dp = [ s[0] - e[0], s[1] - e[1] ] | |
n1 = cp1[0] * cp2[1] - cp1[1] * cp2[0] | |
n2 = s[0] * e[1] - s[1] * e[0] | |
n3 = 1.0 / (dc[0] * dp[1] - dc[1] * dp[0]) | |
return [(n1*dp[0] - n2*dc[0]) * n3, (n1*dp[1] - n2*dc[1]) * n3] | |
outputList = subjectPolygon | |
cp1 = clipPolygon[-1] | |
for clipVertex in clipPolygon: | |
cp2 = clipVertex | |
inputList = outputList | |
outputList = [] | |
if len(inputList) == 0: | |
return None | |
s = inputList[-1] | |
for subjectVertex in inputList: | |
e = subjectVertex | |
if inside(e): | |
if not inside(s): | |
outputList.append(computeIntersection()) | |
outputList.append(e) | |
elif inside(s): | |
outputList.append(computeIntersection()) | |
s = e | |
cp1 = cp2 | |
return (map(lambda x: tuple(x), outputList)) | |
current_way = None | |
current_relation = None | |
tile = 'data/alps.pbf' | |
tile = 'data/mini-uw.pbf' | |
for entity in parse_file(tile): | |
if isinstance(entity, Node): | |
points[entity.id] = {'lat': entity.lat, 'lon': entity.lon} | |
if isinstance(entity, Way): | |
ways[entity.id] = {"nodes": entity.nodes, "tags": entity.tags} | |
if isinstance(entity, Relation): | |
relations.append({"member": entity.members, "tags": entity.tags, "type": "multipolygon"}) | |
img = Image.new('RGB', (texture_size+200, texture_size+200), color=(255,255,255,255)) | |
draw = ImageDraw.Draw(img) | |
clip_rect = [ | |
(texture_border, texture_border ), | |
(texture_border + texture_size, texture_border ), | |
(texture_border + texture_size, texture_border + texture_size), | |
(texture_border , texture_border + texture_size)] | |
def draw_way(way, color): | |
draw_points = [] | |
for node in way["nodes"]: | |
point = geo_to_img(points[node]) | |
draw_points.append(point) | |
solution = clip(draw_points, clip_rect) | |
if solution != None and len(solution) > 2: | |
draw.polygon(tuple(solution), fill=color,outline=color) | |
for relation in relations: | |
if "landuse" in relation["tags"] and relation["tags"]["landuse"] == "forest": | |
for member in relation["member"]: | |
color = (0,200,0,255) | |
if member.role == u"inner": | |
color = (200,200,200,255) | |
way = member.member_id | |
if way in ways: | |
draw_way(ways[way], color) | |
for id in ways: | |
way = ways[id] | |
if 'landuse' in way['tags'] and way['tags'][u'landuse'] == u'forest': | |
draw_way(way, (0,200,0,255)) | |
img.save("test.png") |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment