improved algorithm

README.md

@@ -14,3 +14,37 @@ which will individually be transformed into rectangles which are also plotted.
 On Debian systems you need the following packages:
 
 	apt-get install python python-pil python-scipy python-tk python-matplotlib python-numpy
+
+
+how to setup postgresql+postgis locally without root
+----------------------------------------------------
+
+/usr/lib/postgresql/14/bin/initdb -D /tmp/postgres
+/usr/lib/postgresql/14/bin/postgres --port=5433 --unix_socket_directories=/tmp/postgres -D /tmp/postgres &
+/usr/lib/postgresql/14/bin/createdb --port=5433 --host=/tmp/postgres gis
+/usr/lib/postgresql/14/bin/psql --port=5433 --host=/tmp/postgres gis -c 'CREATE EXTENSION postgis;' -c 'CREATE EXTENSION hstore;'
+
+osm2pgsql --port=5433 --host=/tmp/postgres -d gis --create --slim -G --hstore --tag-transform-script /tmp/openstreetmap-carto/openstreetmap-carto.lua -S /tmp/openstreetmap-carto/openstreetmap-carto.style ~/Downloads/map.xml
+/usr/lib/postgresql/14/bin/psql --port=5433 --host=/tmp/postgres gis -f /tmp/openstreetmap-carto/indexes.sql
+
+openstreetmap-carto:
+    scripts/get-external-data.py --port=5433 --host=/tmp/postgres
+
+diff --git a/project.mml b/project.mml
+index 7fb3d47..d8014f8 100644
+--- a/project.mml
++++ b/project.mml
+@@ -27,6 +27,8 @@ _parts:
+   osm2pgsql: &osm2pgsql
+     type: "postgis"
+     dbname: "gis"
++    port: "5433"
++    host: "/tmp/postgres"
+     key_field: ""
+     geometry_field: "way"
+     extent: "-20037508,-20037508,20037508,20037508"
+
+
+# the database connection settings are part of the style xml!
+carto project.mml > mapnik.xml
+nik4 --url https://www.openstreetmap.org/\#map\=12/49.7731/9.6726 /tmp/openstreetmap-carto/mapnik.xml screenshot.svg

mapbender.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#!/usr/bin/env python3
 #
 # Copyright (C) 2014 - 2021 Johannes Schauer Marin Rodrigues <josch@mister-muffin.de>
 #
@@ -15,32 +15,47 @@
 # You should have received a copy of the GNU General Public License
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
+import os
 import math
 from math import sqrt
-import numpy as np
+import numpy
-import matplotlib.pyplot as plt
 from scipy import interpolate
 from itertools import tee
-from matplotlib.patches import Polygon
+from PIL import Image, ImageDraw
-from matplotlib.collections import PatchCollection
+import urllib.request
-import matplotlib
+import matplotlib.path
-from PIL import Image
+import matplotlib.transforms
+import xml.etree.ElementTree as ET
+
+TILESIZE = 256
+EARTHRADIUS = 6378137
 
 
-def y2lat(a):
+def haversine(lon1, lat1, lon2, lat2):
+    lon1 = math.radians(lon1)
+    lat1 = math.radians(lat1)
+    lon2 = math.radians(lon2)
+    lat2 = math.radians(lat2)
+    dlon = lon2 - lon1
+    dlat = lat2 - lat1
+    a = (
+        math.sin(dlat / 2) ** 2
+        + math.cos(lat1) * math.cos(lat2) * math.sin(dlon / 2) ** 2
+    )
+    return EARTHRADIUS * 2 * math.asin(math.sqrt(a))
+
+
+def lat2y(a, zoom):
     return (
-        180.0
+        (1.0 - math.asinh(math.tan(math.radians(a))) / math.pi)
-        / math.pi
+        / 2.0
-        * (2.0 * math.atan(math.exp(a * math.pi / 180.0)) - math.pi / 2.0)
+        * 2 ** zoom
+        * TILESIZE
     )
 
 
-def lat2y(a):
+def lon2x(a, zoom):
-    return (
+    return (a + 180.0) / 360.0 * (2 ** zoom * TILESIZE)
-        180.0
-        / math.pi
-        * math.log(math.tan(math.pi / 4.0 + a * (math.pi / 180.0) / 2.0))
-    )
 
 
 def pairwise(iterable):
@@ -59,120 +74,54 @@ def triplewise(iterable):
     return zip(a, b, c)
 
 
-# using barycentric coordinates
+def intersects(p0, p1, p2, p3):
-def ptInTriangle(p, p0, p1, p2):
+    s10_x = p1[0] - p0[0]
-    A = 0.5 * (
+    s10_y = p1[1] - p0[1]
-        -p1[1] * p2[0]
+    s32_x = p3[0] - p2[0]
-        + p0[1] * (-p1[0] + p2[0])
+    s32_y = p3[1] - p2[1]
-        + p0[0] * (p1[1] - p2[1])
+
-        + p1[0] * p2[1]
+    denom = s10_x * s32_y - s32_x * s10_y
-    )
+
-    sign = -1 if A < 0 else 1
+    if denom == 0:
-    s = (
+        return False  # collinear
-        p0[1] * p2[0] - p0[0] * p2[1] + (p2[1] - p0[1]) * p[0] + (p0[0] - p2[0]) * p[1]
+
-    ) * sign
+    denom_is_positive = denom > 0
-    t = (
+
-        p0[0] * p1[1] - p0[1] * p1[0] + (p0[1] - p1[1]) * p[0] + (p1[0] - p0[0]) * p[1]
+    s02_x = p0[0] - p2[0]
-    ) * sign
+    s02_y = p0[1] - p2[1]
-    return s >= 0 and t >= 0 and (s + t) <= 2 * A * sign
+
+    s_numer = s10_x * s02_y - s10_y * s02_x
+
+    if (s_numer < 0) == denom_is_positive:
+        return False  # no collision
+
+    t_numer = s32_x * s02_y - s32_y * s02_x
+
+    if (t_numer < 0) == denom_is_positive:
+        return False  # no collision
+
+    if (s_numer > denom) == denom_is_positive or (t_numer > denom) == denom_is_positive:
+        return False  # no collision
+
+    return True
 
 
-def getxing(p0, p1, p2, p3):
+def main(path, width, subdiv, zoom):
-    ux = p1[0] - p0[0]
-    uy = p1[1] - p0[1]
-    vx = p2[0] - p3[0]
-    vy = p2[1] - p3[1]
-    # get multiplicity of u at which u meets v
-    a = vy * ux - vx * uy
-    if a == 0:
-        # lines are parallel and never meet
-        return None
-    s = (vy * (p3[0] - p0[0]) + vx * (p0[1] - p3[1])) / a
-    if 0.0 < s < 1.0:
-        return (p0[0] + s * ux, p0[1] + s * uy)
-    else:
-        return None
-
-
-# the line p0-p1 is the upper normal to the path
-# the line p2-p3 is the lower normal to the path
-#
-#  |        |        |
-# p0--------|--------p1
-#  |        |        |
-#  |        |        |
-# p3--------|--------p2
-#  |        |        |
-def ptInQuadrilateral(p, p0, p1, p2, p3):
-    # it might be that the two normals cross at some point
-    # in that case the two triangles are created differently
-    cross = getxing(p0, p1, p2, p3)
-    if cross:
-        return ptInTriangle(p, p0, cross, p3) or ptInTriangle(p, p2, cross, p1)
-    else:
-        return ptInTriangle(p, p0, p1, p2) or ptInTriangle(p, p2, p3, p0)
-
-
-def get_st(Ax, Ay, Bx, By, Cx, Cy, Dx, Dy, Xx, Xy):
-    d = Bx - Ax - Cx + Dx
-    e = By - Ay - Cy + Dy
-    l = Dx - Ax
-    g = Dy - Ay
-    h = Cx - Dx
-    m = Cy - Dy
-    i = Xx - Dx
-    j = Xy - Dy
-    n = g * h - m * l
-    # calculation for s
-    a1 = m * d - h * e
-    b1 = n - j * d + i * e
-    c1 = l * j - g * i
-    # calculation for t
-    a2 = g * d - l * e
-    b2 = n + j * d - i * e
-    c2 = h * j - m * i
-    s = []
-    if a1 == 0:
-        s.append(-c1 / b1)
-    else:
-        r1 = b1 * b1 - 4 * a1 * c1
-        if r1 >= 0:
-            r11 = (-b1 + sqrt(r1)) / (2 * a1)
-            if -0.0000000001 <= r11 <= 1.0000000001:
-                s.append(r11)
-            r12 = (-b1 - sqrt(r1)) / (2 * a1)
-            if -0.0000000001 <= r12 <= 1.0000000001:
-                s.append(r12)
-    t = []
-    if a2 == 0:
-        t.append(-c2 / b2)
-    else:
-        r2 = b2 * b2 - 4 * a2 * c2
-        if r2 >= 0:
-            r21 = (-b2 + sqrt(r2)) / (2 * a2)
-            if -0.0000000001 <= r21 <= 1.0000000001:
-                t.append(r21)
-            r22 = (-b2 - sqrt(r2)) / (2 * a2)
-            if -0.0000000001 <= r22 <= 1.0000000001:
-                t.append(r22)
-    if not s or not t:
-        return [], []
-    if len(s) == 1 and len(t) == 2:
-        s = [s[0], s[0]]
-    if len(s) == 2 and len(t) == 1:
-        t = [t[0], t[0]]
-    return s, t
-
-
-def main(x, y, width, smoothing, subdiv):
     halfwidth = width / 2.0
-    tck, u = interpolate.splprep([x, y], s=smoothing)
+    found_smoothing = False
-    unew = np.linspace(0, 1.0, subdiv + 1)
+    for smoothing in [2 ** i for i in range(30)]:
+        tck, u = interpolate.splprep(list(zip(*path)), s=smoothing)
+        unew = numpy.linspace(0, 1.0, subdiv + 1)
         out = interpolate.splev(unew, tck)
+        # prepend and append a segment
+        out = (
+            [2 * out[0][0] - out[0][1]] + list(out[0]) + [2 * out[0][-1] - out[0][-2]],
+            [2 * out[1][0] - out[1][1]] + list(out[1]) + [2 * out[1][-1] - out[1][-2]],
+        )
         heights = []
         offs = []
         height = 0.0
-    for (ax, ay), (bx, by) in pairwise(list(zip(*out))):
+        for (ax, ay), (bx, by) in pairwise(zip(*out)):
             s = ax - bx
             t = ay - by
             l = sqrt(s * s + t * t)
@@ -190,7 +139,7 @@ def main(x, y, width, smoothing, subdiv):
         py = [out[1][0] + cy / cl]
         qx = [out[0][0] + dx / dl]
         qy = [out[1][0] + dy / dl]
-    for (ubx, uby), (ux, uy), (uax, uay) in triplewise(list(zip(*out))):
+        for (ubx, uby), (ux, uy), (uax, uay) in triplewise(zip(*out)):
             # get adjacent line segment vectors
             ax = ux - ubx
             ay = uy - uby
@@ -227,52 +176,48 @@ def main(x, y, width, smoothing, subdiv):
         qx.append(out[0][-1] + dx / dl)
         qy.append(out[1][-1] + dy / dl)
         quads = []
-    patches = []
+        for (p2x, p2y, p1x, p1y), (p3x, p3y, p0x, p0y) in pairwise(zip(px, py, qx, qy)):
-    for (p3x, p3y, p2x, p2y), (p0x, p0y, p1x, p1y) in pairwise(
-        list(zip(px, py, qx, qy))
-    ):
             quads.append(((p0x, p0y), (p1x, p1y), (p2x, p2y), (p3x, p3y)))
-        polygon = Polygon(((p0x, p0y), (p1x, p1y), (p2x, p2y), (p3x, p3y)), True)
+        # check for convex quads (sides intersect)
-        patches.append(polygon)
+        have_convex = False
-    containingquad = []
+        for (p0, p1, p2, p3) in quads:
-    for pt in zip(x, y):
+            if intersects(p0, p3, p1, p2):
-        # for each point, find the quadrilateral that contains it
+                have_convex = True
-        found = []
-        for i, (p0, p1, p2, p3) in enumerate(quads):
-            if ptInQuadrilateral(pt, p0, p1, p2, p3):
-                found.append(i)
-        if found:
-            if len(found) > 1:
-                print("point found in two quads")
-                return None
-            containingquad.append(found[0])
-        else:
-            containingquad.append(None)
-    # check if the only points for which no quad could be found are in the
-    # beginning or in the end
-    # find the first missing ones:
-    for i, q in enumerate(containingquad):
-        if q != None:
                 break
-    # find the last missing ones
+        if have_convex:
-    for j, q in zip(range(len(containingquad) - 1, -1, -1), reversed(containingquad)):
+            continue
-        if q != None:
+        ## check for quads that look too much like a triangle
+        # have_triangle = False
+        # for ((p00, p01), (p10, p11), (p20, p21), (p30, p31)) in quads:
+        #    len1 = sqrt((p10-p00)**2+(p11-p01)**2)
+        #    len2 = sqrt((p30-p20)**2+(p31-p21)**2)
+        #    if len1/len2 > 100:
+        #        have_triangle = True
+        #        break
+        #    if len2/len1 < 1/100:
+        #        have_triangle = True
+        #        break
+        # if have_triangle:
+        #    continue
+        # draw
+        polygon = []
+        for ((p00, p01), (p10, p11), (p20, p21), (p30, p31)) in quads:
+            polygon.append((p10, p11))
+            polygon.append((p00, p01))
+        for ((p00, p01), (p10, p11), (p20, p21), (p30, p31)) in reversed(quads):
+            polygon.append((p30, p31))
+            polygon.append((p20, p21))
+        polygon.append(polygon[0])
+
+        # check if path is inside polygon
+        if matplotlib.path.Path(polygon).contains_path(
+            matplotlib.path.Path(path)
+        ):
+            found_smoothing = True
             break
-    # remove the first and last missing ones
+    if not found_smoothing:
-    if i != 0 or j != len(containingquad) - 1:
+        print("cannot find smoothing")
-        containingquad = containingquad[i : j + 1]
+        exit(1)
-        x = x[i : j + 1]
-        y = y[i : j + 1]
-    # check if there are any remaining missing ones:
-    if None in containingquad:
-        print("cannot find quad for point")
-        return None
-    for off, h in zip(offs, heights):
-        targetquad = ((0, off + h), (width, off + h), (width, off), (0, off))
-        patches.append(Polygon(targetquad, True))
-    tx = []
-    ty = []
-    assert len(containingquad) == len(x) == len(y)
     assert (
         len(out[0])
         == len(out[1])
@@ -284,136 +229,176 @@ def main(x, y, width, smoothing, subdiv):
         == len(heights) + 1
         == len(offs) + 1
     )
-    for (rx, ry), i in zip(list(zip(x, y)), containingquad):
+
-        if i == None:
+    minx = math.inf
+    maxx = -1
+    miny = math.inf
+    maxy = -1
+    for (xi, yi) in polygon:
+        if xi < minx:
+            minx = xi
+        if xi > maxx:
+            maxx = xi
+        if yi < miny:
+            miny = yi
+        if yi > maxy:
+            maxy = yi
+    im1 = Image.new("RGB", (int(maxx - minx), int(maxy - miny)))
+    im2 = Image.new("RGB", (int(maxx - minx), int(maxy - miny)))
+    opener = urllib.request.build_opener()
+    opener.addheaders = [("User-agent", "mapbender")]
+    urllib.request.install_opener(opener)
+    todl = []
+    for i in range(int(minx / TILESIZE) - 1, int(maxx / TILESIZE) + 2):
+        for j in range(int(miny / TILESIZE) - 1, int(maxy / TILESIZE) + 2):
+            os.makedirs("%d/%d" % (zoom, i), exist_ok=True)
+            fname = "%d/%d/%d.png" % (zoom, i, j)
+            if not matplotlib.path.Path(numpy.array(polygon)).intersects_bbox(
+                matplotlib.transforms.Bbox(
+                    [
+                        (i * TILESIZE, j * TILESIZE),
+                        (
+                            (i + 1) * TILESIZE,
+                            (j + 1) * TILESIZE,
+                        ),
+                    ]
+                )
+            ):
                 continue
-        (ax, ay), (bx, by), (cx, cy), (dx, dy) = quads[i]
+            if not os.path.exists(fname):
-        s, t = get_st(ax, ay, bx, by, cx, cy, dx, dy, rx, ry)
+                todl.append((i, j))
-        # if more than one solution, take second
+    for n, (i, j) in enumerate(todl):
-        # TODO: investigate if this is always the right solution
+        print("%d/%d" % (n, len(todl)))
-        if len(s) != 1 or len(t) != 1:
+        fname = "%d/%d/%d.png" % (zoom, i, j)
-            s = s[1]
+        urllib.request.urlretrieve(
-            t = t[1]
+            #"https://tile.openstreetmap.org/%d/%d/%d.png" % (zoom, i, j),
-        else:
+            #https://a.tile.thunderforest.com/cycle/17/68690/44518.png?apikey=6170aad10dfd42a38d4d8c709a53
-            s = s[0]
+            "https://tile.thunderforest.com/cycle/%d/%d/%d.png?apikey=d8f470ce7a8e4dd0acf39cc8fd3cf979" % (zoom, i, j),
-            t = t[0]
+            #"https://tile.thunderforest.com/outdoors/%d/%d/%d.png?apikey=d8f470ce7a8e4dd0acf39cc8fd3cf979" % (zoom, i, j),
-        u = s * width
+            #"https://tile.thunderforest.com/landscape/%d/%d/%d.png?apikey=d8f470ce7a8e4dd0acf39cc8fd3cf979" % (zoom, i, j),
-        v = offs[i] + t * heights[i]
+            #"https://tile.thunderforest.com/atlas/%d/%d/%d.png?apikey=d8f470ce7a8e4dd0acf39cc8fd3cf979" % (zoom, i, j),
-        tx.append(u)
+            filename=fname,
-        ty.append(v)
+        )
-    # sx = []
+    for i in range(int(minx / TILESIZE) - 1, int(maxx / TILESIZE) + 2):
-    # sy = []
+        for j in range(int(miny / TILESIZE) - 1, int(maxy / TILESIZE) + 2):
-    # for ((x1,y1),(x2,y2)),((ax,ay),(bx,by),(cx,cy),(dx,dy)),off,h in zip(pairwise(zip(*out)),quads,offs,heights):
+            if not matplotlib.path.Path(numpy.array(polygon)).intersects_bbox(
-    #    s,t = get_st(ax,ay,bx,by,cx,cy,dx,dy,x1,y1)
+                matplotlib.transforms.Bbox(
-    #    if len(s) != 1 or len(t) != 1:
+                    [
-    #        return None
+                        (i * TILESIZE, j * TILESIZE),
-    #    u = s[0]*width
+                        (
-    #    v = off+t[0]*h
+                            (i + 1) * TILESIZE,
-    #    sx.append(u)
+                            (j + 1) * TILESIZE,
-    #    sy.append(v)
+                        ),
-    #    s,t = get_st(ax,ay,bx,by,cx,cy,dx,dy,x2,y2)
+                    ]
-    #    if len(s) != 1 or len(t) != 1:
+                )
-    #        return None
+            ):
-    #    u = s[0]*width
+                continue
-    #    v = off+t[0]*h
+            fname = "%d/%d/%d.png" % (zoom, i, j)
-    #    sx.append(u)
+            with Image.open(fname) as tile:
-    #    sy.append(v)
+                im1.paste(tile, (int(i * TILESIZE - minx), int(j * TILESIZE - miny)))
-    # create map with
+                im2.paste(tile, (int(i * TILESIZE - minx), int(j * TILESIZE - miny)))
-    # python -c 'import logging; logging.basicConfig(level=logging.DEBUG); from landez import ImageExporter; ie = ImageExporter(tiles_url="http://{s}.tile.opencyclemap.org/cycle/{z}/{x}/{y}.png"); ie.export_image(bbox=(8.0419921875,51.25160146817652,10.074462890625,54.03681240523652), zoomlevel=14, imagepath="image.png")'
+    draw2 = ImageDraw.Draw(im2)
-    im = Image.open("map.png")
+    draw2.line([(xi - minx, yi - miny) for xi, yi in path], fill=(255, 0, 0), width=4)
-    bbox = [8.0419921875, 51.25160146817652, 10.074462890625, 54.03681240523652]
+    draw1 = ImageDraw.Draw(im1)
-    # apply mercator projection
+    draw1.line([(xi - minx, yi - miny) for xi, yi in path], fill=(255, 0, 0), width=4)
-    bbox[1] = lat2y(bbox[1])
+    draw1.line([(xi - minx, yi - miny) for xi, yi in zip(*out)], fill=(0, 255, 0))
-    bbox[3] = lat2y(bbox[3])
+    for ((p00, p01), (p10, p11), (p20, p21), (p30, p31)) in quads:
-    iw, ih = im.size
+        draw1.polygon(
+            [
+                (p00 - minx, p01 - miny),
+                (p10 - minx, p11 - miny),
+                (p20 - minx, p21 - miny),
+                (p30 - minx, p31 - miny),
+            ]
+        )
+    draw1.polygon([(xi - minx, yi - miny) for xi, yi in polygon], outline=(0, 0, 255))
+    im1.save("out2.png")
+
     data = []
     for i, (off, h, (p0, p1, p2, p3)) in enumerate(zip(offs, heights, quads)):
-        # first, account for the offset of the input image
+        data.append(
-        p0 = p0[0] - bbox[0], p0[1] - bbox[1]
+            (
-        p1 = p1[0] - bbox[0], p1[1] - bbox[1]
+                (
-        p2 = p2[0] - bbox[0], p2[1] - bbox[1]
-        p3 = p3[0] - bbox[0], p3[1] - bbox[1]
-        # PIL expects coordinates in counter clockwise order
-        p1, p3 = p3, p1
-        #   x          lon
-        # ----- =     -----
-        #   w     bbox[2]-bbox[0]
-        # translate to pixel coordinates
-        p0 = (iw * p0[0]) / (bbox[2] - bbox[0]), (ih * p0[1]) / (bbox[3] - bbox[1])
-        p1 = (iw * p1[0]) / (bbox[2] - bbox[0]), (ih * p1[1]) / (bbox[3] - bbox[1])
-        p2 = (iw * p2[0]) / (bbox[2] - bbox[0]), (ih * p2[1]) / (bbox[3] - bbox[1])
-        p3 = (iw * p3[0]) / (bbox[2] - bbox[0]), (ih * p3[1]) / (bbox[3] - bbox[1])
-        # PIL starts coordinate system at the upper left corner, swap y coord
-        p0 = int(p0[0]), int(ih - p0[1])
-        p1 = int(p1[0]), int(ih - p1[1])
-        p2 = int(p2[0]), int(ih - p2[1])
-        p3 = int(p3[0]), int(ih - p3[1])
-        box = (
                     0,
-            int(ih * (height - off - h) / (bbox[3] - bbox[1])),
+                    int(height - offs[i] - heights[i]),
-            int(iw * width / (bbox[2] - bbox[0])),
+                    int(width),
-            int(ih * (height - off) / (bbox[3] - bbox[1])),
+                    int(height - offs[i]),
+                ),
+                (
+                    p0[0] - minx,
+                    p0[1] - miny,
+                    p1[0] - minx,
+                    p1[1] - miny,
+                    p2[0] - minx,
+                    p2[1] - miny,
+                    p3[0] - minx,
+                    p3[1] - miny,
+                ),
             )
-        quad = (p0[0], p0[1], p1[0], p1[1], p2[0], p2[1], p3[0], p3[1])
+        )
-        data.append((box, quad))
+    im_out = im2.transform(
-    im_out = im.transform(
+        (int(width), int(height)),
-        (int(iw * width / (bbox[2] - bbox[0])), int(ih * height / (bbox[3] - bbox[1]))),
         Image.MESH,
         data,
         Image.BICUBIC,
     )
     im_out.save("out.png")
+    im_out.save("out.jpg", quality=95)
 
-    # np.random.seed(seed=0)
-    # colors = 100*np.random.rand(len(patches)//2)+100*np.random.rand(len(patches)//2)
-    # p = PatchCollection(patches, cmap=matplotlib.cm.jet, alpha=0.4)
-    # p.set_array(np.array(colors))
-    # plt.figure()
-    # plt.axes().set_aspect('equal')
-    ##plt.axhspan(0, height, xmin=0, xmax=width)
-    # fig, ax = plt.subplots()
-    ##ax.add_collection(p)
-    # ax.set_aspect('equal')
-    # plt.axis((0,width,0,height))
-    # plt.imshow(np.asarray(im_out),extent=[0,width,0,height])
-    # plt.imshow(np.asarray(im),extent=[bbox[0],bbox[2],bbox[1],bbox[3]])
-    # plt.plot(x,y,out[0],out[1],px,py,qx,qy,tx,ty)
-    # plt.show()
     return True
 
 
 if __name__ == "__main__":
-    x = []
-    y = []
     import sys
 
-    if len(sys.argv) != 5:
+    if len(sys.argv) != 4:
-        print("usage: %s data.csv width smoothing N" % sys.argv[0])
+        print("usage: %s data.gpx mapwidth paperwidth" % sys.argv[0])
-        print("")
-        print(
-            " data.csv     whitespace delimited lon/lat pairs of points along the path"
-        )
-        print(" width        width of the resulting map in degrees")
-        print(
-            " smoothing    curve smoothing from 0 (exact fit) to higher values (looser fit)"
-        )
-        print(" N            amount of quads to split the path into")
-        print("")
-        print(" example usage:")
-        print("              %s Weser-Radweg-Hauptroute.csv 0.286 6 20" % sys.argv[0])
         exit(1)
+
+    zoom = 10
+    latmin = math.inf
+    latmax = -1
+
+    path = []
     with open(sys.argv[1]) as f:
-        for l in f:
+        root = ET.parse(f)
-            a, b = l.split()
+        for trkpt in root.findall(
+            "./gpx:trk/gpx:trkseg/gpx:trkpt",
+            {"gpx": "http://www.topografix.com/GPX/1/1"},
+        ):
+            lat = float(trkpt.attrib["lat"])
+            lon = float(trkpt.attrib["lon"])
+            if lat < latmin:
+                latmin = lat
+            if lat > latmax:
+                latmax = lat
             # apply mercator projection
-            b = lat2y(float(b))
+            path.append((lon, lat))
-            x.append(float(a))
+
-            y.append(b)
+    length = 0
-    width = float(sys.argv[2])
+    for (lon1, lat1), (lon2, lat2) in pairwise(path):
-    smoothing = float(sys.argv[3])
+        length += haversine(lon1, lat1, lon2, lat2)
-    N = int(sys.argv[4])
+
-    main(x, y, width, smoothing, N)
+    dpi = 96  # because we use bitmap tiles instead of vectors
-    # for smoothing in [1,2,4,8,12]:
+    mapwidthm = float(sys.argv[2])  # map width in m
-    #    for subdiv in range(10,30):
+    paperwidthm = float(sys.argv[3])  # paper width in m
-    #        if main(x,y,width,smoothing,subdiv):
+    earth = 6378137  # earth equator radius in m
-    #            print width,smoothing,subdiv
+    widthpx = dpi / 0.0254 * paperwidthm
+    zoom = math.ceil(
+        math.log2(
+            2
+            * math.pi
+            * earth
+            * math.cos(math.radians((latmax + latmin) / 2))
+            * widthpx
+            / (mapwidthm * TILESIZE)
+        )
+    )
+    subdiv = math.ceil(4*length/mapwidthm)
+    print("zoom:", zoom)
+    print("length:", length)
+    print("subdiv:", subdiv)
+
+    path = [(lon2x(lon, zoom), lat2y(lat, zoom)) for lon, lat in path]
+
+    main(path, widthpx, subdiv, zoom)