You cannot select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

549 lines
19 KiB
Python

#!/usr/bin/env python3
#
# 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
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# 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
from scipy import interpolate
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
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 pairwise(iterable):
"s -> (s0,s1), (s1,s2), (s2,s3), ..."
a, b = tee(iterable, 2)
next(b, None)
return zip(a, b)
def triplewise(iterable):
"s -> (s0,s1,s2), (s1,s2,s3), (s2,s3,s4), ..."
a, b, c = tee(iterable, 3)
next(b, None)
next(c, None)
next(c, None)
return zip(a, b, c)
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():
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 ** 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()