From d3692ff587584669bfd83acc61b3239c7670d6c1 Mon Sep 17 00:00:00 2001 From: Johannes Schauer Marin Rodrigues Date: Thu, 14 Oct 2021 20:43:32 +0200 Subject: [PATCH] improved algorithm --- README.md | 34 +++ mapbender.py | 667 +++++++++++++++++++++++++-------------------------- 2 files changed, 360 insertions(+), 341 deletions(-) diff --git a/README.md b/README.md index da3ec0c..83221ea 100644 --- a/README.md +++ b/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 diff --git a/mapbender.py b/mapbender.py index cc793f2..dc70827 100755 --- a/mapbender.py +++ b/mapbender.py @@ -1,4 +1,4 @@ -#!/usr/bin/env python +#!/usr/bin/env python3 # # Copyright (C) 2014 - 2021 Johannes Schauer Marin Rodrigues # @@ -15,34 +15,49 @@ # You should have received a copy of the GNU General Public License # along with this program. If not, see . +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 -from matplotlib.patches import Polygon -from matplotlib.collections import PatchCollection -import matplotlib -from PIL import Image +from PIL import Image, ImageDraw +import urllib.request +import matplotlib.path +import matplotlib.transforms +import xml.etree.ElementTree as ET +TILESIZE = 256 +EARTHRADIUS = 6378137 -def y2lat(a): - 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(a): +def lat2y(a, zoom): return ( - 180.0 - / math.pi - * math.log(math.tan(math.pi / 4.0 + a * (math.pi / 180.0) / 2.0)) + (1.0 - math.asinh(math.tan(math.radians(a))) / math.pi) + / 2.0 + * 2 ** zoom + * TILESIZE ) +def lon2x(a, zoom): + return (a + 180.0) / 360.0 * (2 ** zoom * TILESIZE) + + def pairwise(iterable): "s -> (s0,s1), (s1,s2), (s2,s3), ..." a, b = tee(iterable, 2) @@ -59,220 +74,150 @@ def triplewise(iterable): 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 - - -# 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): +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] + + denom = s10_x * s32_y - s32_x * s10_y + + 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(path, width, subdiv, zoom): 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(list(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(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 + found_smoothing = False + for smoothing in [2 ** i for i in range(30)]: + tck, u = interpolate.splprep(list(zip(*path)), s=smoothing) + unew = numpy.linspace(0, 1.0, subdiv + 1) + out = interpolate.splev(unew, tck) + # prepend and append a segment + out = ( + [2 * out[0][0] - out[0][1]] + list(out[0]) + [2 * out[0][-1] - out[0][-2]], + [2 * out[1][0] - out[1][1]] + list(out[1]) + [2 * out[1][-1] - out[1][-2]], + ) + heights = [] + offs = [] + height = 0.0 + for (ax, ay), (bx, by) in pairwise(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 - px.append(ux + cx / cl) - py.append(uy + cy / cl) - # and in the other direction - dx = ay + by - dy = -ax - bx + dx = out[1][1] - out[1][0] + dy = -out[0][1] + out[0][0] 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( - list(zip(px, py, qx, qy)) - ): - quads.append(((p0x, p0y), (p1x, p1y), (p2x, p2y), (p3x, p3y))) - polygon = Polygon(((p0x, p0y), (p1x, p1y), (p2x, p2y), (p3x, p3y)), True) - 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 zip(range(len(containingquad) - 1, -1, -1), reversed(containingquad)): - if q != None: + 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: + 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 - # 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) + if not found_smoothing: + print("cannot find smoothing") + exit(1) assert ( len(out[0]) == len(out[1]) @@ -284,136 +229,176 @@ def main(x, y, width, smoothing, subdiv): == len(heights) + 1 == len(offs) + 1 ) - for (rx, ry), i in zip(list(zip(x, y)), containingquad): - if i == None: - 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 + + 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)): - # 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])), + 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, + ), + ) ) - 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]))), + 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) - # np.random.seed(seed=0) - # colors = 100*np.random.rand(len(patches)//2)+100*np.random.rand(len(patches)//2) - # p = PatchCollection(patches, cmap=matplotlib.cm.jet, alpha=0.4) - # p.set_array(np.array(colors)) - # plt.figure() - # plt.axes().set_aspect('equal') - ##plt.axhspan(0, height, xmin=0, xmax=width) - # fig, ax = plt.subplots() - ##ax.add_collection(p) - # ax.set_aspect('equal') - # plt.axis((0,width,0,height)) - # plt.imshow(np.asarray(im_out),extent=[0,width,0,height]) - # plt.imshow(np.asarray(im),extent=[bbox[0],bbox[2],bbox[1],bbox[3]]) - # plt.plot(x,y,out[0],out[1],px,py,qx,qy,tx,ty) - # plt.show() return True if __name__ == "__main__": - x = [] - y = [] import sys - if len(sys.argv) != 5: - 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]) + 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: - for l in f: - a, b = l.split() + 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 - 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 + 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)