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

Compare commits

base: josch:main
Branches Tags
josch:main
..
compare: josch:ddf598c8d9c207e9463f6ae2a731cbb64e1c5003
Branches Tags
josch:main

No commits in common. "main" and "ddf598c8d9c207e9463f6ae2a731cbb64e1c5003" have entirely different histories.
main ... ddf598c8d9

2 changed files with 323 additions and 540 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: On Debian systems you need the following packages:
apt-get install python python-pil python-scipy python-tk python-matplotlib python-numpy 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

711
mapbender.py
View file

@ -1,6 +1,6 @@
#!/usr/bin/env python3 #!/usr/bin/env python
# #
# Copyright (C) 2014 - 2021 Johannes Schauer Marin Rodrigues <josch@mister-muffin.de> # Copyright (C) 2014 Johannes Schauer <j.schauer@email.de>
# #
# This program is free software: you can redistribute it and/or modify # 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 # it under the terms of the GNU General Public License as published by
@ -15,52 +15,28 @@
# You should have received a copy of the GNU General Public License # You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>. # along with this program. If not, see <http://www.gnu.org/licenses/>.
import os
import math import math
from math import sqrt from math import sqrt
import numpy import numpy as np
import matplotlib.pyplot as plt
from scipy import interpolate from scipy import interpolate
from itertools import tee from itertools import tee, izip
from PIL import Image, ImageDraw from matplotlib.patches import Polygon
import urllib.request from matplotlib.collections import PatchCollection
import matplotlib.path import matplotlib
import matplotlib.transforms from PIL import Image
import xml.etree.ElementTree as ET
import mapnik
import cairo
TILESIZE = 256 def y2lat(a):
EARTHRADIUS = 6378137 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(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): def pairwise(iterable):
"s -> (s0,s1), (s1,s2), (s2,s3), ..." "s -> (s0,s1), (s1,s2), (s2,s3), ..."
a, b = tee(iterable, 2) a, b = tee(iterable, 2)
next(b, None) next(b, None)
return zip(a, b) return izip(a, b)
def triplewise(iterable): def triplewise(iterable):
"s -> (s0,s1,s2), (s1,s2,s3), (s2,s3,s4), ..." "s -> (s0,s1,s2), (s1,s2,s3), (s2,s3,s4), ..."
@ -68,131 +44,105 @@ def triplewise(iterable):
next(b, None) next(b, None)
next(c, None) next(c, None)
next(c, None) next(c, None)
return zip(a, b, c) return izip(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 intersects(p0, p1, p2, p3): def getxing(p0, p1, p2, p3):
s10_x = p1[0] - p0[0] ux = p1[0]-p0[0]
s10_y = p1[1] - p0[1] uy = p1[1]-p0[1]
s32_x = p3[0] - p2[0] vx = p2[0]-p3[0]
s32_y = p3[1] - p2[1] 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
denom = s10_x * s32_y - s32_x * s10_y # 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)
if denom == 0: def get_st(Ax,Ay,Bx,By,Cx,Cy,Dx,Dy,Xx,Xy):
return False # collinear d = Bx-Ax-Cx+Dx
e = By-Ay-Cy+Dy
denom_is_positive = denom > 0 l = Dx-Ax
g = Dy-Ay
s02_x = p0[0] - p2[0] h = Cx-Dx
s02_y = p0[1] - p2[1] m = Cy-Dy
i = Xx-Dx
s_numer = s10_x * s02_y - s10_y * s02_x j = Xy-Dy
n = g*h-m*l
if (s_numer < 0) == denom_is_positive: # calculation for s
return False # no collision a1 = m*d-h*e
b1 = n-j*d+i*e
t_numer = s32_x * s02_y - s32_y * s02_x c1 = l*j-g*i
# calculation for t
if (t_numer < 0) == denom_is_positive: a2 = g*d-l*e
return False # no collision b2 = n+j*d-i*e
c2 = h*j-m*i
if (s_numer > denom) == denom_is_positive or (t_numer > denom) == denom_is_positive: s = []
return False # no collision if a1 == 0:
s.append(-c1/b1)
return True else:
r1 = b1*b1-4*a1*c1
if r1 >= 0:
def main(): r11 = (-b1+sqrt(r1))/(2*a1)
import sys if -0.0000000001 <= r11 <= 1.0000000001:
s.append(r11)
if len(sys.argv) != 4: r12 = (-b1-sqrt(r1))/(2*a1)
print("usage: %s data.gpx mapwidth paperwidth" % sys.argv[0]) if -0.0000000001 <= r12 <= 1.0000000001:
exit(1) s.append(r12)
t = []
zoom = 10 if a2 == 0:
latmin = math.inf t.append(-c2/b2)
latmax = -1 else:
r2 = b2*b2-4*a2*c2
path = [] if r2 >= 0:
with open(sys.argv[1]) as f: r21 = (-b2+sqrt(r2))/(2*a2)
root = ET.parse(f) if -0.0000000001 <= r21 <= 1.0000000001:
for trkpt in root.findall( t.append(r21)
"./gpx:trk/gpx:trkseg/gpx:trkpt", r22 = (-b2-sqrt(r2))/(2*a2)
{"gpx": "http://www.topografix.com/GPX/1/1"}, if -0.0000000001 <= r22 <= 1.0000000001:
): t.append(r22)
lat = float(trkpt.attrib["lat"]) if not s or not t:
lon = float(trkpt.attrib["lon"]) return [],[]
if lat < latmin: if len(s) == 1 and len(t) == 2:
latmin = lat s = [s[0],s[0]]
if lat > latmax: if len(s) == 2 and len(t) == 1:
latmax = lat t = [t[0],t[0]]
# apply mercator projection return s, t
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 halfwidth = width/2.0
found_smoothing = False tck,u = interpolate.splprep([x,y],s=smoothing)
#for smoothing in [2 ** i for i in range(30)]: unew = np.linspace(0,1.0,subdiv+1)
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) 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 = [] heights = []
offs = [] offs = []
height = 0.0 height = 0.0
@ -251,298 +201,165 @@ def main():
qx.append(out[0][-1]+dx/dl) qx.append(out[0][-1]+dx/dl)
qy.append(out[1][-1]+dy/dl) qy.append(out[1][-1]+dy/dl)
quads = [] quads = []
for (p2x, p2y, p1x, p1y), (p3x, p3y, p0x, p0y) in pairwise(zip(px, py, qx, qy)): patches = []
for (p3x,p3y,p2x,p2y),(p0x,p0y,p1x,p1y) in pairwise(zip(px,py,qx,qy)):
quads.append(((p0x,p0y),(p1x,p1y),(p2x,p2y),(p3x,p3y))) quads.append(((p0x,p0y),(p1x,p1y),(p2x,p2y),(p3x,p3y)))
# check for convex quads (sides intersect) polygon = Polygon(((p0x,p0y),(p1x,p1y),(p2x,p2y),(p3x,p3y)), True)
have_convex = False patches.append(polygon)
for (p0, p1, p2, p3) in quads: containingquad = []
if intersects(p0, p3, p1, p2): for pt in zip(x,y):
have_convex = True # for each point, find the quadrilateral that contains it
break found = []
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): for i,(p0,p1,p2,p3) in enumerate(quads):
p0 = m.view_transform().forward(mapnik.Coord(*p0)) if ptInQuadrilateral(pt,p0,p1,p2,p3):
p1 = m.view_transform().forward(mapnik.Coord(*p1)) found.append(i)
p2 = m.view_transform().forward(mapnik.Coord(*p2)) if found:
p3 = m.view_transform().forward(mapnik.Coord(*p3)) if len(found) > 1:
if i == 1: print "point found in two quads"
ctx.arc(p0.x, p0.y, 10, 0, 2*math.pi) return None
ctx.set_source_rgb(1,0,0) containingquad.append(found[0])
ctx.fill() else:
ctx.arc(p1.x, p1.y, 10, 0, 2*math.pi) containingquad.append(None)
ctx.set_source_rgb(0,1,0) # check if the only points for which no quad could be found are in the
ctx.fill() # beginning or in the end
ctx.arc(p2.x, p2.y, 10, 0, 2*math.pi) # find the first missing ones:
ctx.set_source_rgb(0,0,1) for i,q in enumerate(containingquad):
ctx.fill() if q != None:
ctx.arc(p3.x, p3.y, 10, 0, 2*math.pi) break
ctx.set_source_rgb(1,1,0) # find the last missing ones
ctx.fill() for j,q in izip(xrange(len(containingquad)-1, -1, -1), reversed(containingquad)):
ctx.move_to(p0.x, p0.y) if q != None:
ctx.line_to(p1.x, p1.y) break
ctx.line_to(p2.x, p2.y) # remove the first and last missing ones
ctx.line_to(p3.x, p3.y) if i != 0 or j != len(containingquad)-1:
ctx.close_path() containingquad = containingquad[i:j+1]
ctx.set_source_rgb(1,0,0) x = x[i:j+1]
ctx.stroke() y = y[i:j+1]
surface.finish() # check if there are any remaining missing ones:
if None in containingquad:
vertices = [] print "cannot find quad for point"
triangles = [] 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:
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)): for i,(off,h,(p0,p1,p2,p3)) in enumerate(zip(offs,heights,quads)):
v1 = (p0[0], p0[1], width, off+h) # top right # first, account for the offset of the input image
v2 = (p1[0], p1[1], width, off) # bottom right p0 = p0[0]-bbox[0],p0[1]-bbox[1]
v3 = (p2[0], p2[1], 0, off) # bottom left p1 = p1[0]-bbox[0],p1[1]-bbox[1]
v4 = (p3[0], p3[1], 0, off+h) # top left p2 = p2[0]-bbox[0],p2[1]-bbox[1]
vertices.extend([v3, v2, v1]) p3 = p3[0]-bbox[0],p3[1]-bbox[1]
triangles.append((len(vertices)-3, len(vertices)-2, len(vertices)-1)) # PIL expects coordinates in counter clockwise order
vertices.extend([v3, v1, v4]) p1,p3 = p3,p1
triangles.append((len(vertices)-3, len(vertices)-2, len(vertices)-1)) # 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
with open("/tmp/tinshift.json", "w") as f: if __name__ == '__main__':
print(""" x = []
{ y = []
"file_type": "triangulation_file", import sys
"format_version": "1.0", if len(sys.argv) != 5:
"transformed_components": [ "horizontal" ], print "usage: %s data.csv width smoothing N"%sys.argv[0]
"fallback_strategy": "nearest", print ""
"vertices_columns": [ "source_x", "source_y", "target_x", "target_y" ], print " data.csv whitespace delimited lon/lat pairs of points along the path"
"triangles_columns": [ "idx_vertex1", "idx_vertex2", "idx_vertex3" ], print " width width of the resulting map in degrees"
"vertices": [ %s ], print " smoothing curve smoothing from 0 (exact fit) to higher values (looser fit)"
"triangles": [ %s ] print " N amount of quads to split the path into"
} print ""
""" % (','.join(["[ %f, %f, %f, %f ]"%v for v in vertices]), ','.join(["[%d, %d, %d]"%t for t in triangles])), file=f) print " example usage:"
print " %s Weser-Radweg-Hauptroute.csv 0.286 6 20"%sys.argv[0]
tinshift = mapnik.Projection("+proj=pipeline +step +proj=webmerc +step +proj=tinshift +file=/tmp/tinshift.json") exit(1)
m = mapnik.Map(1600, int(1600*height/width)) with open(sys.argv[1]) as f:
mapnik.load_map_from_string(m, open("/tmp/openstreetmap-carto/mapnik.xml").read(), False, "/tmp/openstreetmap-carto/") for l in f:
gpxstyle = mapnik.Style() a,b = l.split()
gpxrule = mapnik.Rule() # apply mercator projection
gpxsym = mapnik.LineSymbolizer() b = lat2y(float(b))
gpxsym.stroke = mapnik.Color('blue') x.append(float(a))
gpxsym.stroke_width = 5 y.append(b)
gpxsym.stroke_opacity = 0.5 width = float(sys.argv[2])
gpxrule.symbols.append(gpxsym) smoothing = float(sys.argv[3])
gpxstyle.rules.append(gpxrule) N = int(sys.argv[4])
m.append_style('GPXStyle', gpxstyle) main(x,y,width,smoothing,N)
gpxlayer = mapnik.Layer('GPXLayer') #for smoothing in [1,2,4,8,12]:
gpxlayer.datasource = mapnik.Ogr(file = sys.argv[1], layer = 'tracks') # for subdiv in range(10,30):
gpxlayer.styles.append('GPXStyle') # if main(x,y,width,smoothing,subdiv):
m.layers.append(gpxlayer) # print width,smoothing,subdiv
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()