Implement web mercator projection

map.py

@@ -1,6 +1,7 @@
 from argparse import ArgumentParser
 from collections import namedtuple
 import json
+import math
 from pathlib import Path
 import xml.etree.ElementTree as ET
 
@@ -102,12 +103,12 @@ def main():
     parser = ArgumentParser()
     parser.add_argument('osm', type=Path, help='Path to OSM file to parse')
     parser.add_argument('--convert', action='store_true', default=False, help='Convert data into JSON')
-    parser.add_argument('--map-height', type=int, default=720, help='Height of output map in pixels')
+    parser.add_argument('--map-width', type=int, default=1200, help='Target width of output map in pixels')
     arguments = parser.parse_args()
 
     osm_path: Path = arguments.osm
     convert_data: bool = arguments.convert
-    map_height: int = arguments.map_height
+    map_target_width: int = arguments.map_width
     del arguments
 
     bounds, nodes, ways, relations = parse_xml(osm_path)
@@ -151,12 +152,40 @@ def main():
             json_file.write('}')
         return
 
-    lat_range = bounds.maxlat - bounds.minlat
-    lon_range = bounds.maxlon - bounds.minlon
-    map_ratio = lon_range / lat_range
-    map_width = int(map_height * map_ratio)
-    lat_scale = map_height / lat_range
-    lon_scale = map_width / lon_range
+    zoom_level: int = 0
+    zoom_scale: float = 0.
+
+    def get_map_raw_coords(lon: float, lat: float) -> tuple[int, int]:
+        return (
+            int(zoom_scale * math.radians(lon)),
+            int(zoom_scale * math.log(math.tan(math.pi / 4 + math.radians(lat) / 2)))
+        )
+
+    # def get_gps_raw_coords(x: float, y: float) -> tuple[float, float]:
+    #     return (
+    #         math.degrees(x / zoom_scale),
+    #         math.degrees(2 * math.atan(math.exp(y / zoom_scale)) - math.pi / 2.)
+    #     )
+
+    map_width = 0
+    map_height = 0
+    min_x = 0
+    max_y = 0
+    for zoom_level in range(50):
+        zoom_scale = 2**zoom_level / 2 * math.pi
+        min_x, min_y = get_map_raw_coords(bounds.minlon, bounds.minlat)
+        max_x, max_y = get_map_raw_coords(bounds.maxlon, bounds.maxlat)
+        map_width = max_x - min_x
+        if map_width >= map_target_width:
+            map_height = abs(max_y - min_y)
+            break
+
+    assert map_width > 0 and 0 < map_height < map_width * 3
+
+    def get_map_coords(lon: float, lat: float) -> tuple[int, int]:
+        nonlocal min_x, max_y
+        x, y = get_map_raw_coords(lon, lat)
+        return (x - min_x, max_y - y)
 
     background_color = (235, 235, 235)
     image = Image.new('RGB', (map_width, map_height), background_color)
@@ -164,22 +193,24 @@ def main():
 
     def draw_way_polygon(way: Way, color: tuple[int, ...] | str, outline: tuple[int, ...] | str | None = None,
                          width: int = 1):
-        nonlocal bounds, draw, lat_scale, lon_scale, map_height, map_width, nodes
+        nonlocal bounds, draw, map_height, map_width, nodes
         coords = []
         for node_id in way.nodes:
             node = nodes[node_id]
-            new_coord = ((node.lon - bounds.minlon) * lon_scale, map_height - ((node.lat - bounds.minlat) * lat_scale))
-            if not coords or new_coord != coords[-1]:
+            new_coord = get_map_coords(node.lon, node.lat)
+            # print(f'{new_coord=}')
+            if len(coords) < 2 or new_coord != coords[-1]:
                 coords.append(new_coord)
         draw.polygon(coords, fill=color, outline=outline, width=width)
 
     def draw_way_line(way: Way, width: int, color: tuple[int, ...] | str, joint: str | None = None):
-        nonlocal bounds, draw, lat_scale, lon_scale, map_height, map_width, nodes
+        nonlocal bounds, draw, map_height, map_width, nodes
         coords = []
         for node_id in way.nodes:
             node = nodes[node_id]
-            new_coord = ((node.lon - bounds.minlon) * lon_scale, map_height - ((node.lat - bounds.minlat) * lat_scale))
-            if not coords or new_coord != coords[-1]:
+            new_coord = get_map_coords(node.lon, node.lat)
+            # print(f'{new_coord=}')
+            if len(coords) < 2 or new_coord != coords[-1]:
                 coords.append(new_coord)
         draw.line(coords, width=width, fill=color, joint=joint)