From 845a99eec4bc0cca8354abffd70d06157374ef82 Mon Sep 17 00:00:00 2001 From: Johannes Schauer Marin Rodrigues Date: Sun, 17 Oct 2021 20:10:18 +0200 Subject: [PATCH] using mapnik --- mapbender.py | 476 +++++++++++++++++++++++++++++++++------------------ 1 file changed, 310 insertions(+), 166 deletions(-) diff --git a/mapbender.py b/mapbender.py index dc70827..f09809e 100755 --- a/mapbender.py +++ b/mapbender.py @@ -26,6 +26,8 @@ import urllib.request import matplotlib.path import matplotlib.transforms import xml.etree.ElementTree as ET +import mapnik +import cairo TILESIZE = 256 EARTHRADIUS = 6378137 @@ -45,17 +47,12 @@ def haversine(lon1, lat1, lon2, lat2): return EARTHRADIUS * 2 * math.asin(math.sqrt(a)) -def lat2y(a, zoom): - return ( - (1.0 - math.asinh(math.tan(math.radians(a))) / math.pi) - / 2.0 - * 2 ** zoom - * TILESIZE - ) +def lat2y(lat, zoom): + return math.log(math.tan(math.pi / 4 + math.radians(lat) / 2)) * EARTHRADIUS -def lon2x(a, zoom): - return (a + 180.0) / 360.0 * (2 ** zoom * TILESIZE) +def lon2x(lon, zoom): + return math.radians(lon) * EARTHRADIUS def pairwise(iterable): @@ -106,10 +103,88 @@ def intersects(p0, p1, p2, p3): return True -def main(path, width, subdiv, zoom): +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] + halfwidth = width / 2.0 found_smoothing = False - for smoothing in [2 ** i for i in range(30)]: + #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) @@ -185,6 +260,7 @@ def main(path, width, subdiv, zoom): have_convex = True break if have_convex: + print("have convex") continue ## check for quads that look too much like a triangle # have_triangle = False @@ -215,6 +291,7 @@ def main(path, width, subdiv, zoom): ): found_smoothing = True break + print("doesn't contain path") if not found_smoothing: print("cannot find smoothing") exit(1) @@ -243,162 +320,229 @@ def main(path, width, subdiv, zoom): 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 = [] + 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)): - 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 + 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__": - 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)) - - 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 - widthpx = dpi / 0.0254 * paperwidthm - zoom = math.ceil( - math.log2( - 2 - * math.pi - * earth - * math.cos(math.radians((latmax + latmin) / 2)) - * widthpx - / (mapwidthm * TILESIZE) - ) - ) - subdiv = math.ceil(4*length/mapwidthm) - print("zoom:", zoom) - print("length:", length) - print("subdiv:", subdiv) - - path = [(lon2x(lon, zoom), lat2y(lat, zoom)) for lon, lat in path] - - main(path, widthpx, subdiv, zoom) + main()