using mapnik

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,6 +1,6 @@
-#!/usr/bin/env python
+#!/usr/bin/env python3
 #
-# Copyright (C) 2014  Johannes Schauer <j.schauer@email.de>
+# Copyright (C) 2014 - 2021 Johannes Schauer Marin Rodrigues <josch@mister-muffin.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,28 +15,52 @@
 # 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 matplotlib.pyplot as plt
+import numpy
 from scipy import interpolate
-from itertools import tee, izip
-from matplotlib.patches import Polygon
-from matplotlib.collections import PatchCollection
-import matplotlib
-from PIL import Image
+from itertools import tee
+from PIL import Image, ImageDraw
+import urllib.request
+import matplotlib.path
+import matplotlib.transforms
+import xml.etree.ElementTree as ET
+import mapnik
+import cairo
 
-def y2lat(a):
-    return 180.0/math.pi*(2.0*math.atan(math.exp(a*math.pi/180.0))-math.pi/2.0)
+TILESIZE = 256
+EARTHRADIUS = 6378137
+
+
+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 izip(a, b)
+    return zip(a, b)
+
 
 def triplewise(iterable):
     "s -> (s0,s1,s2), (s1,s2,s3), (s2,s3,s4), ..."
@@ -44,105 +68,131 @@ def triplewise(iterable):
     next(b, None)
     next(c, None)
     next(c, None)
-    return izip(a,b,c)
+    return zip(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 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
+def intersects(p0, p1, p2, p3):
+    s10_x = p1[0] - p0[0]
+    s10_y = p1[1] - p0[1]
+    s32_x = p3[0] - p2[0]
+    s32_y = p3[1] - p2[1]
 
-# 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)
+    denom = s10_x * s32_y - s32_x * s10_y
 
-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
+    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)
+    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]
 
-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)
+    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
@@ -201,165 +251,298 @@ def main(x,y,width,smoothing,subdiv):
         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)):
+        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)))
-        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:
+        # 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
-    # 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:
+        if have_convex:
+            print("have convex")
             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
+        ## 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])
 
-if __name__ == '__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]
+        # 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)
-    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
+    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()