Mirrors/MeshCore-Solo
1
0
Fork 0
mirror of https://github.com/MarekZegare4/MeshCore-Solo.git synced 2026-10-09 19:26:40 +00:00
Code Issues Packages Projects Releases Wiki Activity

feat(ui-lvgl): contour lines on the vector map

Browse Source 
osm_vector.py --dem adds contours from the Terrain Tiles on AWS Open Data
(SRTM, ~25 m; tools/maps/dem.py, standard library only: PNG decode,
smoothing, marching squares, line joining). Every 100 m from zoom 12,
every 20 m from zoom 15; the credit names the elevation source.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
This commit is contained in:
Jakub
2026-09-27 11:35:21 +02:00
co-authored by Claude Opus 5.5
parent 117b7513bc
commit 16144bd48b
4 changed files with 235 additions and 8 deletions
Download Patch File Download Diff File Expand all files Collapse all files
tools/maps/dem.py
+201
View File
@@ -0,0 +1,201 @@
"""Contour lines from a DEM, for osm_vector.py (standard library only).
Elevation comes from the Terrain Tiles on AWS Open Data ("terrarium" PNGs:
height = R * 256 + G + B / 256 - 32768 m; sources SRTM, GMTED, ETOPO1 and
others, see https://github.com/tilezen/joerd/blob/master/docs/attribution.md),
fetched once into a cache folder. The grid is smoothed a little, cut by
marching squares at every `step` metres, and the pieces joined into lines in
world coordinates (Web Mercator 0..1), each with its height.
"""
import math, os, struct, sys, urllib.request, zlib
DEM_ZOOM = 12 # ~25 m a pixel at 49 N: what SRTM has
URL = 'https://s3.amazonaws.com/elevation-tiles-prod/terrarium/{z}/{x}/{y}.png'
ATTRIBUTION = 'elevation: Terrain Tiles (Mapzen / AWS Open Data; SRTM and others)'
def read_png(data):
"""8-bit RGB / RGBA, not interlaced -> (w, h, channels, bytes)."""
assert data[:8] == b'\x89PNG\r\n\x1a\n', 'not a PNG'
pos, idat, w = 8, bytearray(), 0
while pos < len(data):
n, kind = struct.unpack('>I4s', data[pos:pos + 8])
body = data[pos + 8:pos + 8 + n]
if kind == b'IHDR':
w, h, depth, ctype, _, _, interlace = struct.unpack('>IIBBBBB', body)
assert depth == 8 and ctype in (2, 6) and not interlace, 'unsupported PNG'
ch = 3 if ctype == 2 else 4
elif kind == b'IDAT':
idat += body
pos += 12 + n
raw = zlib.decompress(bytes(idat))
stride = w * ch
out = bytearray(stride * h)
prev = bytearray(stride)
for y in range(h):
f = raw[y * (stride + 1)]
line = bytearray(raw[y * (stride + 1) + 1:(y + 1) * (stride + 1)])
if f == 1:
for i in range(ch, stride):
line[i] = (line[i] + line[i - ch]) & 255
elif f == 2:
for i in range(stride):
line[i] = (line[i] + prev[i]) & 255
elif f == 3:
for i in range(stride):
line[i] = (line[i] + ((line[i - ch] if i >= ch else 0) + prev[i]) // 2) & 255
elif f == 4:
for i in range(stride):
a = line[i - ch] if i >= ch else 0
b = prev[i]
c = prev[i - ch] if i >= ch else 0
p = a + b - c
pa, pb, pc = abs(p - a), abs(p - b), abs(p - c)
line[i] = (line[i] + (a if pa <= pb and pa <= pc else b if pb <= pc else c)) & 255
out[y * stride:(y + 1) * stride] = line
prev = line
return w, h, ch, out
def tile(cache, x, y):
path = os.path.join(cache, str(DEM_ZOOM), str(x), f'{y}.png')
if not os.path.exists(path):
os.makedirs(os.path.dirname(path), exist_ok=True)
req = urllib.request.Request(URL.format(z=DEM_ZOOM, x=x, y=y), headers={'User-Agent': 'meshcore-maps/1'})
with urllib.request.urlopen(req, timeout=60) as r, open(path + '.tmp', 'wb') as f:
f.write(r.read())
os.replace(path + '.tmp', path)
w, h, ch, px = read_png(open(path, 'rb').read())
return [[px[(j * w + i) * ch] * 256 + px[(j * w + i) * ch + 1] + px[(j * w + i) * ch + 2] / 256 - 32768
for i in range(w)] for j in range(h)]
def grid(cache, x0, y0, x1, y1):
"""Heights of the tiles x0..x1, y0..y1 stitched: rows of floats."""
rows = []
for ty in range(y0, y1 + 1):
band = [tile(cache, tx, ty) for tx in range(x0, x1 + 1)]
for j in range(256):
rows.append([v for t in band for v in t[j]])
return rows
def smooth(g):
"""3x3 box blur: the 1 px noise of the DEM makes wiggly lines."""
h, w = len(g), len(g[0])
out = [row[:] for row in g]
for j in range(1, h - 1):
a, b, c = g[j - 1], g[j], g[j + 1]
o = out[j]
for i in range(1, w - 1):
o[i] = (a[i - 1] + a[i] + a[i + 1] + b[i - 1] + b[i] + b[i + 1] + c[i - 1] + c[i] + c[i + 1]) / 9
return out
def contours(g, step):
"""Marching squares -> {height: [polyline of (col, row) grid points]}."""
h, w = len(g), len(g[0])
segs = {} # height -> list of (edge key, point, edge key, point)
for j in range(h - 1):
r0, r1 = g[j], g[j + 1]
for i in range(w - 1):
a, b, c, d = r0[i], r0[i + 1], r1[i + 1], r1[i] # corners clockwise from top-left
lo, hi = min(a, b, c, d), max(a, b, c, d)
k0, k1 = math.floor(lo / step) + 1, math.floor(hi / step)
if k0 > k1:
continue
for k in range(k0, k1 + 1):
v = k * step
# Crossings on the cell's edges: top, right, bottom, left.
cr = []
if (a < v) != (b < v):
cr.append((('h', i, j), (i + (v - a) / (b - a), j)))
if (b < v) != (c < v):
cr.append((('v', i + 1, j), (i + 1, j + (v - b) / (c - b))))
if (d < v) != (c < v):
cr.append((('h', i, j + 1), (i + (v - d) / (c - d), j + 1)))
if (a < v) != (d < v):
cr.append((('v', i, j), (i, j + (v - a) / (d - a))))
s = segs.setdefault(v, [])
if len(cr) == 2:
s.append((cr[0], cr[1]))
elif len(cr) == 4: # a saddle: pair by the centre's side
centre = (a + b + c + d) / 4
if (centre < v) == (a < v):
s.append((cr[0], cr[1])); s.append((cr[2], cr[3]))
else:
s.append((cr[0], cr[3])); s.append((cr[1], cr[2]))
out = {}
for v, sl in segs.items():
at = {} # edge key -> segments touching it
for n, (p, q) in enumerate(sl):
at.setdefault(p[0], []).append(n)
at.setdefault(q[0], []).append(n)
used = [False] * len(sl)
lines = []
for n in range(len(sl)):
if used[n]:
continue
used[n] = True
p, q = sl[n]
line = [p, q]
for grow_end in (True, False): # extend from the end, then from the start
while True:
key = line[-1][0] if grow_end else line[0][0]
nxt = next((m for m in at[key] if not used[m]), None)
if nxt is None:
break
used[nxt] = True
a2, b2 = sl[nxt]
far = b2 if a2[0] == key else a2
if grow_end:
line.append(far)
else:
line.insert(0, far)
lines.append([pt for _, pt in line])
out[v] = lines
return out
def chaikin(pts):
"""One round of corner cutting: softer lines than the grid's."""
if len(pts) < 3:
return pts
closed = pts[0] == pts[-1]
out = [] if closed else [pts[0]]
for (x0, y0), (x1, y1) in zip(pts, pts[1:]):
out.append((0.75 * x0 + 0.25 * x1, 0.75 * y0 + 0.25 * y1))
out.append((0.25 * x0 + 0.75 * x1, 0.25 * y0 + 0.75 * y1))
if closed:
out.append(out[0])
else:
out.append(pts[-1])
return out
def contour_lines(cache, lon0, lat0, lon1, lat1, step):
"""[(height, polyline in world coords)] for a box (degrees)."""
n = 1 << DEM_ZOOM
def wxy(lon, lat):
s = math.sin(math.radians(lat))
return (lon + 180.0) / 360.0 * n, (0.5 - math.log((1 + s) / (1 - s)) / (4 * math.pi)) * n
ax, ay = wxy(lon0, lat1)
bx, by = wxy(lon1, lat0)
x0, y0, x1, y1 = int(ax), int(ay), int(bx), int(by)
print(f'DEM: {(x1 - x0 + 1) * (y1 - y0 + 1)} tiles at z{DEM_ZOOM}', file=sys.stderr)
g = smooth(grid(cache, x0, y0, x1, y1))
# Only the box (the tiles reach past it): its pixels, plus one.
c0, r0 = max(0, int((ax - x0) * 256) - 1), max(0, int((ay - y0) * 256) - 1)
c1, r1 = min(len(g[0]), int((bx - x0) * 256) + 2), min(len(g), int((by - y0) * 256) + 2)
g = [row[c0:c1] for row in g[r0:r1]]
size = 256 * n
out = []
for v, lines in contours(g, step).items():
for line in lines:
if len(line) < 4:
continue
line = chaikin(line)
out.append((v, [((x0 * 256 + c0 + c + 0.5) / size, (y0 * 256 + r0 + r + 0.5) / size) for c, r in line]))
return out
tools/maps/osm_vector.py
+23 -3
View File
@@ -7,7 +7,10 @@ marked hiking routes in their waymark colours, water, forest, meadows, rock,
buildings; and named points for the labels (places, peaks, huts, springs...).
No street names (the raster map has them).
tools/maps/osm_vector.py area.json [points.json] --out vmap/
tools/maps/osm_vector.py area.json [points.json] --out vmap/ [--dem dem-cache/]
With --dem, contour lines too (tools/maps/dem.py: Terrain Tiles fetched into
that folder): every 20 m (from zoom 15) with a darker one every 100 m.
Three data zooms: 10 (drawn at z10-11), 12 (z12-13), 14 (z14-18); the device picks the
data tile covering the tile it draws and scales it. Copy the output folder to
@@ -56,6 +59,8 @@ and the points (a second file):
Map data (c) OpenStreetMap contributors, ODbL.
"""
import argparse, json, math, os, struct, sys
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
import dem
from collections import defaultdict
EXTENT = 4096
@@ -69,6 +74,7 @@ L_PATH, L_TRACK, L_SERVICE, L_MINOR, L_TERTIARY, L_SECONDARY, L_PRIMARY, L_TRUNK
L_PATH_HARD = 29 # a path of demanding / alpine difficulty (sac_scale)
L_ROUTE = 50 # (old: one route, its RGB565 colour)
L_ROUTES = 51 # the routes along a stretch
L_CONTOUR, L_CONTOUR_IDX = 15, 16 # contour lines, every 20 m / 100 m
HIGHWAY = {
'path': L_PATH, 'footway': L_PATH, 'steps': L_PATH, 'bridleway': L_PATH, 'cycleway': L_PATH,
@@ -81,7 +87,8 @@ HIGHWAY = {
'trunk': L_TRUNK, 'trunk_link': L_TRUNK, 'motorway': L_TRUNK, 'motorway_link': L_TRUNK,
}
# Lowest data zoom a class goes in (smaller ones would be clutter / weight).
MIN_DZ = {P_BUILDING: 14, L_SERVICE: 14, L_PATH: 12, L_PATH_HARD: 12, L_TRACK: 12, L_STREAM: 12, L_MINOR: 12}
MIN_DZ = {P_BUILDING: 14, L_SERVICE: 14, L_PATH: 12, L_PATH_HARD: 12, L_TRACK: 12, L_STREAM: 12, L_MINOR: 12,
L_CONTOUR: 14, L_CONTOUR_IDX: 12}
HARD_SAC = ('demanding_mountain_hiking', 'alpine_hiking', 'demanding_alpine_hiking', 'difficult_alpine_hiking')
# Points (labels): class, lowest data zoom.
@@ -147,7 +154,7 @@ def short_name(name):
# Natural land cover (not roads, buildings, water): its edges are vague
# anyway, so simplified harder -- most of the points are here.
NATURAL = (P_MEADOW, P_SCRUB, P_FOREST, P_ROCK)
NATURAL = (P_MEADOW, P_SCRUB, P_FOREST, P_ROCK, L_CONTOUR, L_CONTOUR_IDX) # (and the contours: from a 25 m grid)
NATURAL_TOL = 3
WAYMARK = { # osmc:symbol / colour -> RGB888
@@ -390,6 +397,9 @@ def main():
ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
ap.add_argument('json', nargs='+')
ap.add_argument('--out', default='vmap')
ap.add_argument('--dem', metavar='CACHE', help='add contour lines; DEM tiles are kept in this folder')
ap.add_argument('--contour-step', type=int, default=20)
ap.add_argument('--index-step', type=int, default=100)
a = ap.parse_args()
elements = []
for fn in a.json:
@@ -446,6 +456,14 @@ def main():
for cols, ways in bundles.items():
packed = sum(c << (4 * i) for i, c in enumerate(cols))
feats.append((L_ROUTES, packed, 'line', join_lines(ways)))
if a.dem:
lons = [p['lon'] for e in elements for p in (e.get('geometry') or []) if p]
lats = [p['lat'] for e in elements for p in (e.get('geometry') or []) if p]
nc = 0
for v, line in dem.contour_lines(a.dem, min(lons), min(lats), max(lons), max(lats), a.contour_step):
feats.append((L_CONTOUR_IDX if round(v) % a.index_step == 0 else L_CONTOUR, 0, 'line', [line]))
nc += 1
print(f'{nc} contour lines', file=sys.stderr)
print(f'{len(feats)} features, {len(way_routes)} route ways', file=sys.stderr)
total_bytes = total_tiles = 0
@@ -539,6 +557,8 @@ def main():
print(f'{len(points)} points', file=sys.stderr)
with open(os.path.join(a.out, 'attribution.txt'), 'w') as f:
f.write('© OpenStreetMap contributors (ODbL)\n')
if a.dem:
f.write(dem.ATTRIBUTION + '\n')
print(f'{total_tiles} tiles, {total_bytes / 1024:.0f} KB', file=sys.stderr)