using mapnik

main
parent d3692ff587
commit 845a99eec4
Signed by: josch
GPG Key ID: F2CBA5C78FBD83E1

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

Loading…
Cancel
Save