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,6 +1,6 @@
-#!/usr/bin/env python3
+#!/usr/bin/env python
 #
-# Copyright (C) 2014 - 2021 Johannes Schauer Marin Rodrigues <josch@mister-muffin.de>
+# Copyright (C) 2014  Johannes Schauer <j.schauer@email.de>
 #
 # This program is free software: you can redistribute it and/or modify
 # it under the terms of the GNU General Public License as published by
@@ -15,534 +15,351 @@
 # 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 itertools import tee, izip
-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
+def y2lat(a):
-EARTHRADIUS = 6378137
+    return 180.0/math.pi*(2.0*math.atan(math.exp(a*math.pi/180.0))-math.pi/2.0)
-
-
-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(lat, zoom):
-    return math.log(math.tan(math.pi / 4 + math.radians(lat) / 2)) * EARTHRADIUS
-
-
-def lon2x(lon, zoom):
-    return math.radians(lon) * EARTHRADIUS
 
+def lat2y(a):
+    return 180.0/math.pi*math.log(math.tan(math.pi/4.0+a*(math.pi/180.0)/2.0))
 
 def pairwise(iterable):
     "s -> (s0,s1), (s1,s2), (s2,s3), ..."
     a, b = tee(iterable, 2)
     next(b, None)
-    return zip(a, b)
+    return izip(a, b)
-
 
 def triplewise(iterable):
     "s -> (s0,s1,s2), (s1,s2,s3), (s2,s3,s4), ..."
-    a, b, c = tee(iterable, 3)
+    a,b,c = tee(iterable, 3)
     next(b, None)
     next(c, None)
     next(c, None)
-    return zip(a, b, c)
+    return izip(a,b,c)
 
+# using barycentric coordinates
+def ptInTriangle(p, p0, p1, p2):
+    A = 0.5 * (-p1[1] * p2[0] + p0[1] * (-p1[0] + p2[0]) + p0[0] * (p1[1] - p2[1]) + p1[0] * p2[1]);
+    sign = -1 if A < 0 else 1;
+    s = (p0[1] * p2[0] - p0[0] * p2[1] + (p2[1] - p0[1]) * p[0] + (p0[0] - p2[0]) * p[1]) * sign;
+    t = (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;
 
-def intersects(p0, p1, p2, p3):
+def getxing(p0, p1, p2, p3):
-    s10_x = p1[0] - p0[0]
+    ux = p1[0]-p0[0]
-    s10_y = p1[1] - p0[1]
+    uy = p1[1]-p0[1]
-    s32_x = p3[0] - p2[0]
+    vx = p2[0]-p3[0]
-    s32_y = p3[1] - p2[1]
+    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
 
-    denom = s10_x * s32_y - s32_x * s10_y
+# 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)
 
-    if denom == 0:
+def get_st(Ax,Ay,Bx,By,Cx,Cy,Dx,Dy,Xx,Xy):
-        return False  # collinear
+    d = Bx-Ax-Cx+Dx
-
+    e = By-Ay-Cy+Dy
-    denom_is_positive = denom > 0
+    l = Dx-Ax
-
+    g = Dy-Ay
-    s02_x = p0[0] - p2[0]
+    h = Cx-Dx
-    s02_y = p0[1] - p2[1]
+    m = Cy-Dy
-
+    i = Xx-Dx
-    s_numer = s10_x * s02_y - s10_y * s02_x
+    j = Xy-Dy
-
+    n = g*h-m*l
-    if (s_numer < 0) == denom_is_positive:
+    # calculation for s
-        return False  # no collision
+    a1 = m*d-h*e
-
+    b1 = n-j*d+i*e
-    t_numer = s32_x * s02_y - s32_y * s02_x
+    c1 = l*j-g*i
-
+    # calculation for t
-    if (t_numer < 0) == denom_is_positive:
+    a2 = g*d-l*e
-        return False  # no collision
+    b2 = n+j*d-i*e
-
+    c2 = h*j-m*i
-    if (s_numer > denom) == denom_is_positive or (t_numer > denom) == denom_is_positive:
+    s = []
-        return False  # no collision
+    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)
+    unew = np.linspace(0,1.0,subdiv+1)
+    out = interpolate.splev(unew,tck)
+    heights = []
+    offs = []
+    height = 0.0
+    for (ax,ay),(bx,by) in pairwise(zip(*out)):
+        s = ax-bx
+        t = ay-by
+        l = sqrt(s*s+t*t)
+        offs.append(height)
+        height += l
+        heights.append(l)
+    # the border of the first segment is just perpendicular to the path
+    cx = -out[1][1]+out[1][0]
+    cy = out[0][1]-out[0][0]
+    cl = sqrt(cx*cx+cy*cy)/halfwidth
+    dx = out[1][1]-out[1][0]
+    dy = -out[0][1]+out[0][0]
+    dl = sqrt(dx*dx+dy*dy)/halfwidth
+    px = [out[0][0]+cx/cl]
+    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(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
+        px.append(ux+cx/cl)
+        py.append(uy+cy/cl)
+        # and in the other direction
+        dx = ay+by
+        dy = -ax-bx
+        dl = sqrt(dx*dx+dy*dy)/halfwidth
+        qx.append(ux+dx/dl)
+        qy.append(uy+dy/dl)
+    # the border of the last segment is just perpendicular to the path
+    cx = -out[1][-1]+out[1][-2]
+    cy = out[0][-1]-out[0][-2]
+    cl = sqrt(cx*cx+cy*cy)/halfwidth
+    dx = out[1][-1]-out[1][-2]
+    dy = -out[0][-1]+out[0][-2]
+    dl = sqrt(dx*dx+dy*dy)/halfwidth
+    px.append(out[0][-1]+cx/cl)
+    py.append(out[1][-1]+cy/cl)
+    qx.append(out[0][-1]+dx/dl)
+    qy.append(out[1][-1]+dy/dl)
+    quads = []
+    patches = []
+    for (p3x,p3y,p2x,p2y),(p0x,p0y,p1x,p1y) in pairwise(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)
+        patches.append(polygon)
+    containingquad = []
+    for pt in zip(x,y):
+        # for each point, find the quadrilateral that contains it
+        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
+    for j,q in izip(xrange(len(containingquad)-1, -1, -1), reversed(containingquad)):
+        if q != None:
+            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]) == len(px) == len(py) == len(qx) == len(qy) == len(quads)+1 == len(heights)+1 == len(offs)+1
+    for (rx,ry),i in zip(zip(x,y),containingquad):
+        if i == None:
+            continue
+        (ax,ay),(bx,by),(cx,cy),(dx,dy) = quads[i]
+        s,t = get_st(ax,ay,bx,by,cx,cy,dx,dy,rx,ry)
+        # if more than one solution, take second 
+        # TODO: investigate if this is always the right solution
+        if len(s) != 1 or len(t) != 1:
+            s = s[1]
+            t = t[1]
+        else:
+            s = s[0]
+            t = t[0]
+        u = s*width
+        v = offs[i]+t*heights[i]
+        tx.append(u)
+        ty.append(v)
+    #sx = []
+    #sy = []
+    #for ((x1,y1),(x2,y2)),((ax,ay),(bx,by),(cx,cy),(dx,dy)),off,h in zip(pairwise(zip(*out)),quads,offs,heights):
+    #    s,t = get_st(ax,ay,bx,by,cx,cy,dx,dy,x1,y1)
+    #    if len(s) != 1 or len(t) != 1:
+    #        return None
+    #    u = s[0]*width
+    #    v = off+t[0]*h
+    #    sx.append(u)
+    #    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
+    #    v = off+t[0]*h
+    #    sx.append(u)
+    #    sy.append(v)
+    # create map with
+    # 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")'
+    im = Image.open("map.png")
+    bbox = [8.0419921875,51.25160146817652,10.074462890625,54.03681240523652]
+    # apply mercator projection
+    bbox[1] = lat2y(bbox[1])
+    bbox[3] = lat2y(bbox[3])
+    iw,ih = im.size
+    data = []
+    for i,(off,h,(p0,p1,p2,p3)) in enumerate(zip(offs,heights,quads)):
+        # first, account for the offset of the input image
+        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(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")
+    #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__':
-def 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:
-        root = ET.parse(f)
+        for l in f:
-        for trkpt in root.findall(
+            a,b = l.split()
-            "./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))
+            b = lat2y(float(b))
-
+            x.append(float(a))
-    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')
+            y.append(b)
-    longlat = mapnik.Projection('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')
+    width = float(sys.argv[2])
-    #m = mapnik.Map(800, 600)
+    smoothing = float(sys.argv[3])
-    #mapnik.load_map_from_string(m, open("/tmp/openstreetmap-carto/mapnik.xml").read(), False, "/tmp/openstreetmap-carto/")
+    N = int(sys.argv[4])
-    #gpxstyle = mapnik.Style()
+    main(x,y,width,smoothing,N)
-    #gpxrule = mapnik.Rule()
+    #for smoothing in [1,2,4,8,12]:
-    #gpxsym = mapnik.LineSymbolizer()
+    #    for subdiv in range(10,30):
-    #gpxsym.stroke = mapnik.Color('red')
+    #        if main(x,y,width,smoothing,subdiv):
-    #gpxsym.stroke_width = 5
+    #            print width,smoothing,subdiv
-    #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)
-    subdiv = 40
-    width = 15000
-    print("zoom:", zoom)
-    print("length:", length)
-    print("subdiv:", subdiv)
-
-    path = [(lon2x(lon, zoom), lat2y(lat, zoom)) for lon, lat in path]
-
-    halfwidth = width / 2.0
-    found_smoothing = False
-    #for smoothing in [2 ** i for i in range(30)]:
-    for smoothing in [2 ** 35]:
-        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(zip(*out)):
-            s = ax - bx
-            t = ay - by
-            l = sqrt(s * s + t * t)
-            offs.append(height)
-            height += l
-            heights.append(l)
-        # the border of the first segment is just perpendicular to the path
-        cx = -out[1][1] + out[1][0]
-        cy = out[0][1] - out[0][0]
-        cl = sqrt(cx * cx + cy * cy) / halfwidth
-        dx = out[1][1] - out[1][0]
-        dy = -out[0][1] + out[0][0]
-        dl = sqrt(dx * dx + dy * dy) / halfwidth
-        px = [out[0][0] + cx / cl]
-        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(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
-            px.append(ux + cx / cl)
-            py.append(uy + cy / cl)
-            # and in the other direction
-            dx = ay + by
-            dy = -ax - bx
-            dl = sqrt(dx * dx + dy * dy) / halfwidth
-            qx.append(ux + dx / dl)
-            qy.append(uy + dy / dl)
-        # the border of the last segment is just perpendicular to the path
-        cx = -out[1][-1] + out[1][-2]
-        cy = out[0][-1] - out[0][-2]
-        cl = sqrt(cx * cx + cy * cy) / halfwidth
-        dx = out[1][-1] - out[1][-2]
-        dy = -out[0][-1] + out[0][-2]
-        dl = sqrt(dx * dx + dy * dy) / halfwidth
-        px.append(out[0][-1] + cx / cl)
-        py.append(out[1][-1] + cy / cl)
-        qx.append(out[0][-1] + dx / dl)
-        qy.append(out[1][-1] + dy / dl)
-        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")
-    if not found_smoothing:
-        print("cannot find smoothing")
-        exit(1)
-    assert (
-        len(out[0])
-        == len(out[1])
-        == len(px)
-        == len(py)
-        == len(qx)
-        == len(qy)
-        == len(quads) + 1
-        == len(heights) + 1
-        == len(offs) + 1
-    )
-
-    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
-    print(minx, maxx, miny, maxy)
-
-    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('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.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.zoom_to_box(mapnik.Box2d(mapnik.Coord(minx, miny), mapnik.Coord(maxx, maxy)))
-    m.srs = merc.params()
-    surface = cairo.SVGSurface("out.svg", 800, 600)
-    ctx = cairo.Context(surface)
-    mapnik.render(m, ctx)
-    for pos in zip(*out):
-        pos = m.view_transform().forward(mapnik.Coord(*pos))
-        ctx.line_to(pos.x, pos.y)
-    ctx.set_source_rgb(0,0,0)
-    ctx.stroke()
-    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
-        v2 = (p1[0], p1[1], width, off)   # bottom right
-        v3 = (p2[0], p2[1], 0, off)       # bottom left
-        v4 = (p3[0], p3[1], 0, off+h)     # top left
-        vertices.extend([v3, v2, v1])
-        triangles.append((len(vertices)-3, len(vertices)-2, len(vertices)-1))
-        vertices.extend([v3, v1, v4])
-        triangles.append((len(vertices)-3, len(vertices)-2, len(vertices)-1))
-
-    with open("/tmp/tinshift.json", "w") as f:
-        print("""
-{
-  "file_type": "triangulation_file",
-  "format_version": "1.0",
-  "transformed_components": [ "horizontal" ],
-  "fallback_strategy": "nearest",
-  "vertices_columns": [ "source_x", "source_y", "target_x", "target_y" ],
-  "triangles_columns": [ "idx_vertex1", "idx_vertex2", "idx_vertex3" ],
-  "vertices": [ %s ],
-  "triangles": [ %s ]
-}
-""" % (','.join(["[ %f, %f, %f, %f ]"%v for v in vertices]), ','.join(["[%d, %d, %d]"%t for t in triangles])), file=f)
-
-    tinshift = mapnik.Projection("+proj=pipeline +step +proj=webmerc +step +proj=tinshift +file=/tmp/tinshift.json")
-    m = mapnik.Map(1600, int(1600*height/width))
-    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()