|
| 1 | +#!/usr/bin/env python |
| 2 | + |
| 3 | +# Author: Florian Winkler (Fju) 2018 |
| 4 | +# This file is part of the overpass-api-python-wrapper project |
| 5 | +# which is licensed under Apache 2.0. |
| 6 | +# See LICENSE.txt for the full license text. |
| 7 | + |
| 8 | +import os |
| 9 | +import xml.etree.ElementTree |
| 10 | +import matplotlib.pyplot as plt |
| 11 | +import numpy as np |
| 12 | +import overpass |
| 13 | + |
| 14 | +MAX_SIZE = 400 |
| 15 | +XML_FILE = 'state_border_saxony.xml' |
| 16 | +QUERY = """area[name="Sachsen"][type="boundary"]->.saxony; |
| 17 | + rel(area.saxony)[admin_level=4][type="boundary"][boundary="administrative"]; |
| 18 | + out geom;""" |
| 19 | + |
| 20 | +# the query is can be quite big (multiple MB) so we save the servers response in a file |
| 21 | +# that can be re-opened the next time |
| 22 | +if not(os.path.isfile(XML_FILE)): |
| 23 | + # increase timeout because the response can be pretty heavy |
| 24 | + api = overpass.API(timeout=60) |
| 25 | + |
| 26 | + response = api.get(QUERY, responseformat="xml") |
| 27 | + |
| 28 | + # write xml file |
| 29 | + f = open(XML_FILE, 'w') |
| 30 | + # encode response in UTF-8 because name translation contain non-ascii characters |
| 31 | + f.write(response.encode('utf-8')) |
| 32 | + f.close() |
| 33 | + |
| 34 | + # free up memory |
| 35 | + #del response |
| 36 | + |
| 37 | +# open xml file |
| 38 | +root = xml.etree.ElementTree.parse(XML_FILE).getroot() |
| 39 | + |
| 40 | +# bounds element contains information of the width and height |
| 41 | +bounds = root.find('relation').find('bounds') |
| 42 | +min_lat = float(bounds.get('minlat')) |
| 43 | +min_lon = float(bounds.get('minlon')) |
| 44 | +max_lat = float(bounds.get('maxlat')) |
| 45 | +max_lon = float(bounds.get('maxlon')) |
| 46 | + |
| 47 | +# longitude: east - west (x) |
| 48 | +# latitude: north - south (y) |
| 49 | +box_width = max_lon - min_lon |
| 50 | +box_height = max_lat - min_lat |
| 51 | + |
| 52 | +# compute scale factors so that the biggest distance is equal to `MAX_SIZE` |
| 53 | +scale_x = int(MAX_SIZE * box_width / max(box_width, box_height)) |
| 54 | +scale_y = int(MAX_SIZE * box_height / max(box_width, box_height)) |
| 55 | + |
| 56 | +def outside(x1, y1, x2, y2, tolerance): |
| 57 | + # check if distance between two points is bigger than the given tolerance |
| 58 | + # since tolerance is a constant it doesn't have to be squared |
| 59 | + dx = x1 - x2 |
| 60 | + dy = y1 - y2 |
| 61 | + return (dx*dx + dy*dy) > tolerance |
| 62 | + |
| 63 | + |
| 64 | +point_count = 0 |
| 65 | +# look through all `member` nodes |
| 66 | +for member in root.iter('member'): |
| 67 | + m_type, m_role = member.get('type'), member.get('role') |
| 68 | + |
| 69 | + # check if the element belongs to the outer border |
| 70 | + if m_type == 'way' and m_role == 'outer': |
| 71 | + # `nd` elements contain lon and lat coordinates |
| 72 | + m_points = member.findall('nd') |
| 73 | + |
| 74 | + prev_x, prev_y = -1, -1 |
| 75 | + |
| 76 | + index = 0 |
| 77 | + for mp in m_points: |
| 78 | + x, y = float(mp.get('lon')), float(mp.get('lat')) |
| 79 | + |
| 80 | + # convert lon and lat coordinates to pixel coordinates |
| 81 | + x = (x - min_lon) / box_width * scale_x |
| 82 | + y = (y - min_lat) / box_height * scale_y |
| 83 | + |
| 84 | + if index == 0: |
| 85 | + # first point |
| 86 | + prev_x = x |
| 87 | + prev_y = y |
| 88 | + elif outside(x, y, prev_x, prev_y, 1) or index == len(m_points) - 1: |
| 89 | + # check if points are not too close to each other (too much detail) or if it's the last point of the section |
| 90 | + # if the last point was ignored, there would be holes |
| 91 | + |
| 92 | + # draw line from current point to previous point |
| 93 | + plt.plot([x, prev_x], [y, prev_y], 'k-') |
| 94 | + point_count += 1 |
| 95 | + # save coordinates |
| 96 | + prev_x = x |
| 97 | + prev_y = y |
| 98 | + |
| 99 | + index += 1 |
| 100 | + |
| 101 | +print(point_count) |
| 102 | +# show plot |
| 103 | +plt.show() |
0 commit comments