README.md

@@ -14,37 +14,3 @@ 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 python3
+#!/usr/bin/env python
 #
 # Copyright (C) 2014 - 2021 Johannes Schauer Marin Rodrigues <josch@mister-muffin.de>
 #
@@ -15,44 +15,32 @@
 # 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
+import numpy as np
+import matplotlib.pyplot as plt
 from scipy import interpolate
 from itertools import tee
-from PIL import Image, ImageDraw
+from matplotlib.patches import Polygon
-import urllib.request
+from matplotlib.collections import PatchCollection
-import matplotlib.path
+import matplotlib
-import matplotlib.transforms
+from PIL import Image
-import xml.etree.ElementTree as ET
-import mapnik
-import cairo
-
-TILESIZE = 256
-EARTHRADIUS = 6378137
 
 
-def haversine(lon1, lat1, lon2, lat2):
+def y2lat(a):
-    lon1 = math.radians(lon1)
+    return (
-    lat1 = math.radians(lat1)
+        180.0
-    lon2 = math.radians(lon2)
+        / math.pi
-    lat2 = math.radians(lat2)
+        * (2.0 * math.atan(math.exp(a * math.pi / 180.0)) - math.pi / 2.0)
-    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(lat, zoom):
+def lat2y(a):
-    return math.log(math.tan(math.pi / 4 + math.radians(lat) / 2)) * EARTHRADIUS
+    return (
-
+        180.0
-
+        / math.pi
-def lon2x(lon, zoom):
+        * math.log(math.tan(math.pi / 4.0 + a * (math.pi / 180.0) / 2.0))
-    return math.radians(lon) * EARTHRADIUS
+    )
 
 
 def pairwise(iterable):
@@ -71,230 +59,220 @@ def triplewise(iterable):
     return zip(a, b, c)
 
 
-def intersects(p0, p1, p2, p3):
+# using barycentric coordinates
-    s10_x = p1[0] - p0[0]
+def ptInTriangle(p, p0, p1, p2):
-    s10_y = p1[1] - p0[1]
+    A = 0.5 * (
-    s32_x = p3[0] - p2[0]
+        -p1[1] * p2[0]
-    s32_y = p3[1] - p2[1]
+        + p0[1] * (-p1[0] + p2[0])
-
+        + p0[0] * (p1[1] - p2[1])
-    denom = s10_x * s32_y - s32_x * s10_y
+        + p1[0] * p2[1]
-
-    if denom == 0:
-        return False  # collinear
-
-    denom_is_positive = denom > 0
-
-    s02_x = p0[0] - p2[0]
-    s02_y = p0[1] - p2[1]
-
-    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 main():
-    import sys
-
-    if len(sys.argv) != 4:
-        print("usage: %s data.gpx mapwidth paperwidth" % sys.argv[0])
-        exit(1)
-
-    zoom = 10
-    latmin = math.inf
-    latmax = -1
-
-    path = []
-    with open(sys.argv[1]) as f:
-        root = ET.parse(f)
-        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
-            path.append((lon, lat))
-
-    merc = mapnik.Projection('+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs +over')
-    longlat = mapnik.Projection('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')
-    #m = mapnik.Map(800, 600)
-    #mapnik.load_map_from_string(m, open("/tmp/openstreetmap-carto/mapnik.xml").read(), False, "/tmp/openstreetmap-carto/")
-    #gpxstyle = mapnik.Style()
-    #gpxrule = mapnik.Rule()
-    #gpxsym = mapnik.LineSymbolizer()
-    #gpxsym.stroke = mapnik.Color('red')
-    #gpxsym.stroke_width = 5
-    #gpxsym.stroke_opacity = 0.5
-    #gpxrule.symbols.append(gpxsym)
-    #gpxstyle.rules.append(gpxrule)
-    #m.append_style('GPXStyle', gpxstyle)
-    #gpxlayer = mapnik.Layer('GPXLayer')
-    ##gpxlayer.srs =
-    #gpxlayer.datasource = mapnik.Ogr(file = sys.argv[1], layer = 'tracks')
-    #gpxlayer.styles.append('GPXStyle')
-    #m.layers.append(gpxlayer)
-    #m.aspect_fix_mode = mapnik.aspect_fix_mode.GROW_BBOX
-    #m.srs = merc.params()
-    #m.zoom_to_box(mapnik.ProjTransform(longlat, merc).forward(gpxlayer.envelope()))
-    #mapnik.render_to_file(m, "out.png", "png", 1.0)
-
-    length = 0
-    for (lon1, lat1), (lon2, lat2) in pairwise(path):
-        length += haversine(lon1, lat1, lon2, lat2)
-
-    dpi = 96  # because we use bitmap tiles instead of vectors
-    mapwidthm = float(sys.argv[2])  # map width in m
-    paperwidthm = float(sys.argv[3])  # paper width in m
-    earth = 6378137  # earth equator radius in m
-    width = dpi / 0.0254 * paperwidthm # width in px
-    zoom = math.ceil(
-        math.log2(
-            2
-            * math.pi
-            * earth
-            * math.cos(math.radians((latmax + latmin) / 2))
-            * width
-            / (mapwidthm * TILESIZE)
-        )
     )
-    subdiv = math.ceil(4*length/mapwidthm)
+    sign = -1 if A < 0 else 1
-    subdiv = 40
+    s = (
-    width = 15000
+        p0[1] * p2[0] - p0[0] * p2[1] + (p2[1] - p0[1]) * p[0] + (p0[0] - p2[0]) * p[1]
-    print("zoom:", zoom)
+    ) * sign
-    print("length:", length)
+    t = (
-    print("subdiv:", subdiv)
+        p0[0] * p1[1] - p0[1] * p1[0] + (p0[1] - p1[1]) * p[0] + (p1[0] - p0[0]) * p[1]
+    ) * sign
+    return s >= 0 and t >= 0 and (s + t) <= 2 * A * sign
 
-    path = [(lon2x(lon, zoom), lat2y(lat, zoom)) for lon, lat in path]
 
+def getxing(p0, p1, p2, p3):
+    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
-    found_smoothing = False
+    tck, u = interpolate.splprep([x, y], s=smoothing)
-    #for smoothing in [2 ** i for i in range(30)]:
+    unew = np.linspace(0, 1.0, subdiv + 1)
-    for smoothing in [2 ** 35]:
+    out = interpolate.splev(unew, tck)
-        tck, u = interpolate.splprep(list(zip(*path)), s=smoothing)
+    heights = []
-        unew = numpy.linspace(0, 1.0, subdiv + 1)
+    offs = []
-        out = interpolate.splev(unew, tck)
+    height = 0.0
-        # prepend and append a segment
+    for (ax, ay), (bx, by) in pairwise(list(zip(*out))):
-        out = (
+        s = ax - bx
-            [2 * out[0][0] - out[0][1]] + list(out[0]) + [2 * out[0][-1] - out[0][-2]],
+        t = ay - by
-            [2 * out[1][0] - out[1][1]] + list(out[1]) + [2 * out[1][-1] - out[1][-2]],
+        l = sqrt(s * s + t * t)
-        )
+        offs.append(height)
-        heights = []
+        height += l
-        offs = []
+        heights.append(l)
-        height = 0.0
+    # the border of the first segment is just perpendicular to the path
-        for (ax, ay), (bx, by) in pairwise(zip(*out)):
+    cx = -out[1][1] + out[1][0]
-            s = ax - bx
+    cy = out[0][1] - out[0][0]
-            t = ay - by
+    cl = sqrt(cx * cx + cy * cy) / halfwidth
-            l = sqrt(s * s + t * t)
+    dx = out[1][1] - out[1][0]
-            offs.append(height)
+    dy = -out[0][1] + out[0][0]
-            height += l
+    dl = sqrt(dx * dx + dy * dy) / halfwidth
-            heights.append(l)
+    px = [out[0][0] + cx / cl]
-        # the border of the first segment is just perpendicular to the path
+    py = [out[1][0] + cy / cl]
-        cx = -out[1][1] + out[1][0]
+    qx = [out[0][0] + dx / dl]
-        cy = out[0][1] - out[0][0]
+    qy = [out[1][0] + dy / dl]
+    for (ubx, uby), (ux, uy), (uax, uay) in triplewise(list(zip(*out))):
+        # get adjacent line segment vectors
+        ax = ux - ubx
+        ay = uy - uby
+        bx = uax - ux
+        by = uay - uy
+        # normalize length
+        al = sqrt(ax * ax + ay * ay)
+        bl = sqrt(bx * bx + by * by)
+        ax = ax / al
+        ay = ay / al
+        bx = bx / bl
+        by = by / bl
+        # get vector perpendicular to sum
+        cx = -ay - by
+        cy = ax + bx
         cl = sqrt(cx * cx + cy * cy) / halfwidth
-        dx = out[1][1] - out[1][0]
+        px.append(ux + cx / cl)
-        dy = -out[0][1] + out[0][0]
+        py.append(uy + cy / cl)
+        # and in the other direction
+        dx = ay + by
+        dy = -ax - bx
         dl = sqrt(dx * dx + dy * dy) / halfwidth
-        px = [out[0][0] + cx / cl]
+        qx.append(ux + dx / dl)
-        py = [out[1][0] + cy / cl]
+        qy.append(uy + dy / dl)
-        qx = [out[0][0] + dx / dl]
+    # the border of the last segment is just perpendicular to the path
-        qy = [out[1][0] + dy / dl]
+    cx = -out[1][-1] + out[1][-2]
-        for (ubx, uby), (ux, uy), (uax, uay) in triplewise(zip(*out)):
+    cy = out[0][-1] - out[0][-2]
-            # get adjacent line segment vectors
+    cl = sqrt(cx * cx + cy * cy) / halfwidth
-            ax = ux - ubx
+    dx = out[1][-1] - out[1][-2]
-            ay = uy - uby
+    dy = -out[0][-1] + out[0][-2]
-            bx = uax - ux
+    dl = sqrt(dx * dx + dy * dy) / halfwidth
-            by = uay - uy
+    px.append(out[0][-1] + cx / cl)
-            # normalize length
+    py.append(out[1][-1] + cy / cl)
-            al = sqrt(ax * ax + ay * ay)
+    qx.append(out[0][-1] + dx / dl)
-            bl = sqrt(bx * bx + by * by)
+    qy.append(out[1][-1] + dy / dl)
-            ax = ax / al
+    quads = []
-            ay = ay / al
+    patches = []
-            bx = bx / bl
+    for (p3x, p3y, p2x, p2y), (p0x, p0y, p1x, p1y) in pairwise(
-            by = by / bl
+        list(zip(px, py, qx, qy))
-            # get vector perpendicular to sum
+    ):
-            cx = -ay - by
+        quads.append(((p0x, p0y), (p1x, p1y), (p2x, p2y), (p3x, p3y)))
-            cy = ax + bx
+        polygon = Polygon(((p0x, p0y), (p1x, p1y), (p2x, p2y), (p3x, p3y)), True)
-            cl = sqrt(cx * cx + cy * cy) / halfwidth
+        patches.append(polygon)
-            px.append(ux + cx / cl)
+    containingquad = []
-            py.append(uy + cy / cl)
+    for pt in zip(x, y):
-            # and in the other direction
+        # for each point, find the quadrilateral that contains it
-            dx = ay + by
+        found = []
-            dy = -ax - bx
+        for i, (p0, p1, p2, p3) in enumerate(quads):
-            dl = sqrt(dx * dx + dy * dy) / halfwidth
+            if ptInQuadrilateral(pt, p0, p1, p2, p3):
-            qx.append(ux + dx / dl)
+                found.append(i)
-            qy.append(uy + dy / dl)
+        if found:
-        # the border of the last segment is just perpendicular to the path
+            if len(found) > 1:
-        cx = -out[1][-1] + out[1][-2]
+                print("point found in two quads")
-        cy = out[0][-1] - out[0][-2]
+                return None
-        cl = sqrt(cx * cx + cy * cy) / halfwidth
+            containingquad.append(found[0])
-        dx = out[1][-1] - out[1][-2]
+        else:
-        dy = -out[0][-1] + out[0][-2]
+            containingquad.append(None)
-        dl = sqrt(dx * dx + dy * dy) / halfwidth
+    # check if the only points for which no quad could be found are in the
-        px.append(out[0][-1] + cx / cl)
+    # beginning or in the end
-        py.append(out[1][-1] + cy / cl)
+    # find the first missing ones:
-        qx.append(out[0][-1] + dx / dl)
+    for i, q in enumerate(containingquad):
-        qy.append(out[1][-1] + dy / dl)
+        if q != None:
-        quads = []
-        for (p2x, p2y, p1x, p1y), (p3x, p3y, p0x, p0y) in pairwise(zip(px, py, qx, qy)):
-            quads.append(((p0x, p0y), (p1x, p1y), (p2x, p2y), (p3x, p3y)))
-        # check for convex quads (sides intersect)
-        have_convex = False
-        for (p0, p1, p2, p3) in quads:
-            if intersects(p0, p3, p1, p2):
-                have_convex = True
-                break
-        if have_convex:
-            print("have convex")
-            continue
-        ## 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
-        print("doesn't contain path")
+    # find the last missing ones
-    if not found_smoothing:
+    for j, q in zip(range(len(containingquad) - 1, -1, -1), reversed(containingquad)):
-        print("cannot find smoothing")
+        if q != None:
-        exit(1)
+            break
+    # remove the first and last missing ones
+    if i != 0 or j != len(containingquad) - 1:
+        containingquad = containingquad[i : j + 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])
@@ -306,243 +284,136 @@ def main():
         == len(heights) + 1
         == len(offs) + 1
     )
-
+    for (rx, ry), i in zip(list(zip(x, y)), containingquad):
-    minx = math.inf
+        if i == None:
-    maxx = -1
+            continue
-    miny = math.inf
+        (ax, ay), (bx, by), (cx, cy), (dx, dy) = quads[i]
-    maxy = -1
+        s, t = get_st(ax, ay, bx, by, cx, cy, dx, dy, rx, ry)
-    for (xi, yi) in polygon:
+        # if more than one solution, take second
-        if xi < minx:
+        # TODO: investigate if this is always the right solution
-            minx = xi
+        if len(s) != 1 or len(t) != 1:
-        if xi > maxx:
+            s = s[1]
-            maxx = xi
+            t = t[1]
-        if yi < miny:
+        else:
-            miny = yi
+            s = s[0]
-        if yi > maxy:
+            t = t[0]
-            maxy = yi
+        u = s * width
-    print(minx, maxx, miny, maxy)
+        v = offs[i] + t * heights[i]
-
+        tx.append(u)
-    m = mapnik.Map(800, 600)
+        ty.append(v)
-    mapnik.load_map_from_string(m, open("/tmp/openstreetmap-carto/mapnik.xml").read(), False, "/tmp/openstreetmap-carto/")
+    # sx = []
-    gpxstyle = mapnik.Style()
+    # sy = []
-    gpxrule = mapnik.Rule()
+    # for ((x1,y1),(x2,y2)),((ax,ay),(bx,by),(cx,cy),(dx,dy)),off,h in zip(pairwise(zip(*out)),quads,offs,heights):
-    gpxsym = mapnik.LineSymbolizer()
+    #    s,t = get_st(ax,ay,bx,by,cx,cy,dx,dy,x1,y1)
-    gpxsym.stroke = mapnik.Color('blue')
+    #    if len(s) != 1 or len(t) != 1:
-    gpxsym.stroke_width = 5
+    #        return None
-    gpxsym.stroke_opacity = 0.5
+    #    u = s[0]*width
-    gpxrule.symbols.append(gpxsym)
+    #    v = off+t[0]*h
-    gpxstyle.rules.append(gpxrule)
+    #    sx.append(u)
-    m.append_style('GPXStyle', gpxstyle)
+    #    sy.append(v)
-    gpxlayer = mapnik.Layer('GPXLayer')
+    #    s,t = get_st(ax,ay,bx,by,cx,cy,dx,dy,x2,y2)
-    #gpxlayer.srs =
+    #    if len(s) != 1 or len(t) != 1:
-    gpxlayer.datasource = mapnik.Ogr(file = sys.argv[1], layer = 'tracks')
+    #        return None
-    gpxlayer.styles.append('GPXStyle')
+    #    u = s[0]*width
-    m.layers.append(gpxlayer)
+    #    v = off+t[0]*h
-    m.aspect_fix_mode = mapnik.aspect_fix_mode.GROW_BBOX
+    #    sx.append(u)
-    m.zoom_to_box(mapnik.Box2d(mapnik.Coord(minx, miny), mapnik.Coord(maxx, maxy)))
+    #    sy.append(v)
-    m.srs = merc.params()
+    # create map with
-    surface = cairo.SVGSurface("out.svg", 800, 600)
+    # 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")'
-    ctx = cairo.Context(surface)
+    im = Image.open("map.png")
-    mapnik.render(m, ctx)
+    bbox = [8.0419921875, 51.25160146817652, 10.074462890625, 54.03681240523652]
-    for pos in zip(*out):
+    # apply mercator projection
-        pos = m.view_transform().forward(mapnik.Coord(*pos))
+    bbox[1] = lat2y(bbox[1])
-        ctx.line_to(pos.x, pos.y)
+    bbox[3] = lat2y(bbox[3])
-    ctx.set_source_rgb(0,0,0)
+    iw, ih = im.size
-    ctx.stroke()
+    data = []
-    for i, (p0, p1, p2, p3) in enumerate(quads):
-        p0 = m.view_transform().forward(mapnik.Coord(*p0))
-        p1 = m.view_transform().forward(mapnik.Coord(*p1))
-        p2 = m.view_transform().forward(mapnik.Coord(*p2))
-        p3 = m.view_transform().forward(mapnik.Coord(*p3))
-        if i == 1:
-            ctx.arc(p0.x, p0.y, 10, 0, 2*math.pi)
-            ctx.set_source_rgb(1,0,0)
-            ctx.fill()
-            ctx.arc(p1.x, p1.y, 10, 0, 2*math.pi)
-            ctx.set_source_rgb(0,1,0)
-            ctx.fill()
-            ctx.arc(p2.x, p2.y, 10, 0, 2*math.pi)
-            ctx.set_source_rgb(0,0,1)
-            ctx.fill()
-            ctx.arc(p3.x, p3.y, 10, 0, 2*math.pi)
-            ctx.set_source_rgb(1,1,0)
-            ctx.fill()
-        ctx.move_to(p0.x, p0.y)
-        ctx.line_to(p1.x, p1.y)
-        ctx.line_to(p2.x, p2.y)
-        ctx.line_to(p3.x, p3.y)
-        ctx.close_path()
-        ctx.set_source_rgb(1,0,0)
-        ctx.stroke()
-    surface.finish()
-
-    vertices = []
-    triangles = []
     for i, (off, h, (p0, p1, p2, p3)) in enumerate(zip(offs, heights, quads)):
-        v1 = (p0[0], p0[1], width, off+h) # top right
+        # first, account for the offset of the input image
-        v2 = (p1[0], p1[1], width, off)   # bottom right
+        p0 = p0[0] - bbox[0], p0[1] - bbox[1]
-        v3 = (p2[0], p2[1], 0, off)       # bottom left
+        p1 = p1[0] - bbox[0], p1[1] - bbox[1]
-        v4 = (p3[0], p3[1], 0, off+h)     # top left
+        p2 = p2[0] - bbox[0], p2[1] - bbox[1]
-        vertices.extend([v3, v2, v1])
+        p3 = p3[0] - bbox[0], p3[1] - bbox[1]
-        triangles.append((len(vertices)-3, len(vertices)-2, len(vertices)-1))
+        # PIL expects coordinates in counter clockwise order
-        vertices.extend([v3, v1, v4])
+        p1, p3 = p3, p1
-        triangles.append((len(vertices)-3, len(vertices)-2, len(vertices)-1))
+        #   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(iw * width / (bbox[2] - bbox[0])),
+            int(ih * (height - off) / (bbox[3] - bbox[1])),
+        )
+        quad = (p0[0], p0[1], p1[0], p1[1], p2[0], p2[1], p3[0], p3[1])
+        data.append((box, quad))
+    im_out = im.transform(
+        (int(iw * width / (bbox[2] - bbox[0])), int(ih * height / (bbox[3] - bbox[1]))),
+        Image.MESH,
+        data,
+        Image.BICUBIC,
+    )
+    im_out.save("out.png")
 
-    with open("/tmp/tinshift.json", "w") as f:
+    # np.random.seed(seed=0)
-        print("""
+    # 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)
-  "file_type": "triangulation_file",
+    # p.set_array(np.array(colors))
-  "format_version": "1.0",
+    # plt.figure()
-  "transformed_components": [ "horizontal" ],
+    # plt.axes().set_aspect('equal')
-  "fallback_strategy": "nearest",
+    ##plt.axhspan(0, height, xmin=0, xmax=width)
-  "vertices_columns": [ "source_x", "source_y", "target_x", "target_y" ],
+    # fig, ax = plt.subplots()
-  "triangles_columns": [ "idx_vertex1", "idx_vertex2", "idx_vertex3" ],
+    ##ax.add_collection(p)
-  "vertices": [ %s ],
+    # ax.set_aspect('equal')
-  "triangles": [ %s ]
+    # plt.axis((0,width,0,height))
-}
+    # plt.imshow(np.asarray(im_out),extent=[0,width,0,height])
-""" % (','.join(["[ %f, %f, %f, %f ]"%v for v in vertices]), ','.join(["[%d, %d, %d]"%t for t in triangles])), file=f)
+    # 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)
-    tinshift = mapnik.Projection("+proj=pipeline +step +proj=webmerc +step +proj=tinshift +file=/tmp/tinshift.json")
+    # plt.show()
-    m = mapnik.Map(1600, int(1600*height/width))
+    return True
-    mapnik.load_map_from_string(m, open("/tmp/openstreetmap-carto/mapnik.xml").read(), False, "/tmp/openstreetmap-carto/")
-    gpxstyle = mapnik.Style()
-    gpxrule = mapnik.Rule()
-    gpxsym = mapnik.LineSymbolizer()
-    gpxsym.stroke = mapnik.Color('blue')
-    gpxsym.stroke_width = 5
-    gpxsym.stroke_opacity = 0.5
-    gpxrule.symbols.append(gpxsym)
-    gpxstyle.rules.append(gpxrule)
-    m.append_style('GPXStyle', gpxstyle)
-    gpxlayer = mapnik.Layer('GPXLayer')
-    gpxlayer.datasource = mapnik.Ogr(file = sys.argv[1], layer = 'tracks')
-    gpxlayer.styles.append('GPXStyle')
-    m.layers.append(gpxlayer)
-    m.aspect_fix_mode = mapnik.aspect_fix_mode.GROW_BBOX
-    m.zoom_to_box(mapnik.Box2d(mapnik.Coord(0, 0), mapnik.Coord(width, height)))
-    m.srs = tinshift.params()
-    surface = cairo.PDFSurface("out.pdf", 1600, int(1600*height/width))
-    mapnik.render(m, surface)
-    surface.finish()
-
-
-    #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
-    #        if not os.path.exists(fname):
-    #            todl.append((i, j))
-    #for n, (i, j) in enumerate(todl):
-    #    print("%d/%d" % (n, len(todl)))
-    #    fname = "%d/%d/%d.png" % (zoom, i, j)
-    #    urllib.request.urlretrieve(
-    #        #"https://tile.openstreetmap.org/%d/%d/%d.png" % (zoom, i, j),
-    #        #https://a.tile.thunderforest.com/cycle/17/68690/44518.png?apikey=6170aad10dfd42a38d4d8c709a53
-    #        "https://tile.thunderforest.com/cycle/%d/%d/%d.png?apikey=d8f470ce7a8e4dd0acf39cc8fd3cf979" % (zoom, i, j),
-    #        #"https://tile.thunderforest.com/outdoors/%d/%d/%d.png?apikey=d8f470ce7a8e4dd0acf39cc8fd3cf979" % (zoom, i, j),
-    #        #"https://tile.thunderforest.com/landscape/%d/%d/%d.png?apikey=d8f470ce7a8e4dd0acf39cc8fd3cf979" % (zoom, i, j),
-    #        #"https://tile.thunderforest.com/atlas/%d/%d/%d.png?apikey=d8f470ce7a8e4dd0acf39cc8fd3cf979" % (zoom, i, j),
-    #        filename=fname,
-    #    )
-    #for i in range(int(minx / TILESIZE) - 1, int(maxx / TILESIZE) + 2):
-    #    for j in range(int(miny / TILESIZE) - 1, int(maxy / TILESIZE) + 2):
-    #        if not matplotlib.path.Path(numpy.array(polygon)).intersects_bbox(
-    #            matplotlib.transforms.Bbox(
-    #                [
-    #                    (i * TILESIZE, j * TILESIZE),
-    #                    (
-    #                        (i + 1) * TILESIZE,
-    #                        (j + 1) * TILESIZE,
-    #                    ),
-    #                ]
-    #            )
-    #        ):
-    #            continue
-    #        fname = "%d/%d/%d.png" % (zoom, i, j)
-    #        with Image.open(fname) as tile:
-    #            im1.paste(tile, (int(i * TILESIZE - minx), int(j * TILESIZE - miny)))
-    #            im2.paste(tile, (int(i * TILESIZE - minx), int(j * TILESIZE - miny)))
-    #draw2 = ImageDraw.Draw(im2)
-    #draw2.line([(xi - minx, yi - miny) for xi, yi in path], fill=(255, 0, 0), width=4)
-    #draw1 = ImageDraw.Draw(im1)
-    #draw1.line([(xi - minx, yi - miny) for xi, yi in path], fill=(255, 0, 0), width=4)
-    #draw1.line([(xi - minx, yi - miny) for xi, yi in zip(*out)], fill=(0, 255, 0))
-    #for ((p00, p01), (p10, p11), (p20, p21), (p30, p31)) in quads:
-    #    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)):
-    #    data.append(
-    #        (
-    #            (
-    #                0,
-    #                int(height - offs[i] - heights[i]),
-    #                int(width),
-    #                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,
-    #            ),
-    #        )
-    #    )
-    #im_out = im2.transform(
-    #    (int(width), int(height)),
-    #    Image.MESH,
-    #    data,
-    #    Image.BICUBIC,
-    #)
-    #im_out.save("out.png")
-    #im_out.save("out.jpg", quality=95)
-
-    #return True
 
 
 if __name__ == "__main__":
-    main()
+    x = []
+    y = []
+    import sys
+
+    if len(sys.argv) != 5:
+        print("usage: %s data.csv width smoothing N" % 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)
+    with open(sys.argv[1]) as f:
+        for l in f:
+            a, b = l.split()
+            # apply mercator projection
+            b = lat2y(float(b))
+            x.append(float(a))
+            y.append(b)
+    width = float(sys.argv[2])
+    smoothing = float(sys.argv[3])
+    N = int(sys.argv[4])
+    main(x, y, width, smoothing, N)
+    # for smoothing in [1,2,4,8,12]:
+    #    for subdiv in range(10,30):
+    #        if main(x,y,width,smoothing,subdiv):
+    #            print width,smoothing,subdiv