josch/mapbender
1
0
Fork 0
Code Issues Pull requests Projects Releases Wiki Activity

Compare commits

base: josch:845a99eec4bc0cca8354abffd70d06157374ef82
Branches Tags
josch:main
..
compare: josch:a4ced15215d32ef96f5a35e5691d9870bca2d920
Branches Tags
josch:main

No commits in common. "845a99eec4bc0cca8354abffd70d06157374ef82" and "a4ced15215d32ef96f5a35e5691d9870bca2d920" have entirely different histories.
845a99eec4 ... a4ced15215

2 changed files with 353 additions and 516 deletions
Show stats Expand all files Collapse all files

34
README.md
View file

@ -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

835
mapbender.py
View file

@ -1,4 +1,4 @@
#!/usr/bin/env python3
#!/usr/bin/env python
#
# Copyright (C) 2014 - 2021 Johannes Schauer Marin Rodrigues <josch@mister-muffin.de>
#
@ -15,44 +15,32 @@
# 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 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
from matplotlib.patches import Polygon
from matplotlib.collections import PatchCollection
import matplotlib
from PIL import Image
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
def y2lat(a):
return (
180.0
/ math.pi
* (2.0 * math.atan(math.exp(a * math.pi / 180.0)) - math.pi / 2.0)
)
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):
@ -71,230 +59,220 @@ def triplewise(iterable):
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)
)
# 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]
)
subdiv = math.ceil(4*length/mapwidthm)
subdiv = 40
width = 15000
print("zoom:", zoom)
print("length:", length)
print("subdiv:", subdiv)
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
path = [(lon2x(lon, zoom), lat2y(lat, zoom)) for lon, lat in path]
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):
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]
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
cl = sqrt(cx * cx + cy * cy) / halfwidth
dx = out[1][1] - out[1][0]
dy = -out[0][1] + out[0][0]
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
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
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
print("doesn't contain path")
if not found_smoothing:
print("cannot find smoothing")
exit(1)
# find the last missing ones
for j, q in zip(range(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])
@ -306,243 +284,136 @@ def main():
== 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 (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
data = []
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))
# 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")
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
# 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__":
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])
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