diff --git a/docs/development/l2-roadmap.md b/docs/development/l2-roadmap.md index a488c720..69bb0026 100644 --- a/docs/development/l2-roadmap.md +++ b/docs/development/l2-roadmap.md @@ -264,8 +264,11 @@ Ideas to come back to (2026-09-26): white band, less simplified; demanding / alpine paths dotted; labels of named points (places, peaks with height, huts, passes, springs...) as a layer over the tiles, placed by priority. No street names (raster). -- [ ] Vector maps next: contours, a PBF pipeline for regions + download on - the device, styling. +- [x] Contours on the vector map: osm_vector.py --dem (tools/maps/dem.py, + standard library: Terrain Tiles from AWS Open Data, SRTM ~25 m, + smoothed, marching squares); 100 m lines from z12, 20 m from z15. +- [ ] Vector maps next: a PBF pipeline for regions + download on the + device, styling (anti-aliasing, dark theme). User list, second batch (2026-09-26): diff --git a/examples/companion_radio/ui-lvgl/map/VectorTileProvider.h b/examples/companion_radio/ui-lvgl/map/VectorTileProvider.h index 335b3a3a..c9cccdd2 100644 --- a/examples/companion_radio/ui-lvgl/map/VectorTileProvider.h +++ b/examples/companion_radio/ui-lvgl/map/VectorTileProvider.h @@ -47,7 +47,7 @@ public: return true; } - const char* attribution() const override { return _have ? "\xC2\xA9 OpenStreetMap contributors (ODbL)" : _fb.attribution(); } + const char* attribution() const override { return _have ? "\xC2\xA9 OpenStreetMap contributors (ODbL); contours: Terrain Tiles (SRTM, AWS Open Data)" : _fb.attribution(); } bool hasData() const { return _have; } uint32_t lastMs() const { return _last_ms; } // The last tile's time split: reading the data, areas, lines (ms). @@ -131,7 +131,7 @@ private: } // Styles. Widths in px at zoom 14, scaled with the zoom. - enum : uint8_t { L_STREAM = 20, L_RIVER = 21, L_PATH_HARD = 29, L_PATH = 30, L_TRACK = 31, L_SERVICE = 32, L_TRUNK = 37, + enum : uint8_t { L_CONTOUR = 15, L_CONTOUR_IDX = 16, L_STREAM = 20, L_RIVER = 21, L_PATH_HARD = 29, L_PATH = 30, L_TRACK = 31, L_SERVICE = 32, L_TRUNK = 37, L_ROUTE = 50, L_ROUTES = 51 }; // Waymark colours of L_ROUTES (index 1.., osm_vector.py PALETTE). static uint16_t routeColour(int i) { @@ -152,6 +152,8 @@ private: } static bool lineStyle(uint8_t c, uint16_t& col, float& w) { switch (c) { + case L_CONTOUR: col = rgb(0xB89668); w = 1.0f; return true; + case L_CONTOUR_IDX: col = rgb(0x9A7448); w = 1.0f; return true; case L_STREAM: col = rgb(0x86B6D8); w = 1.0f; return true; case L_RIVER: col = rgb(0x86B6D8); w = 3.0f; return true; case L_PATH: col = rgb(0xA8502A); w = 1.2f; return true; @@ -196,7 +198,8 @@ private: bool road = cls >= L_SERVICE && cls <= L_TRUNK; bool want = pass == 0 ? (cls < L_SERVICE) : pass == 3 ? (cls == L_ROUTE || cls == L_ROUTES) : road; if (!want) continue; - if (_z < 13 && (cls == L_PATH || cls == L_PATH_HARD || cls == 7)) continue; // paths from z13, buildings from z14 (data) + if (_z < 13 && (cls == L_PATH || cls == L_PATH_HARD || cls == 7)) continue; + if (_z < 15 && cls == L_CONTOUR) continue; // every 20 m: too dense further out // paths from z13, buildings from z14 (data) // Off the drawn tile (with a margin for the widest line): not even read. const int32_t M = 12 * 16, S = TILE_PX * 16; if (sx(bb[2]) < -M || sx(bb[0]) > S + M || sy(bb[3]) < -M || sy(bb[1]) > S + M) continue; diff --git a/tools/maps/dem.py b/tools/maps/dem.py new file mode 100644 index 00000000..93ecce48 --- /dev/null +++ b/tools/maps/dem.py @@ -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 diff --git a/tools/maps/osm_vector.py b/tools/maps/osm_vector.py index e41e2dc2..9e2d0348 100644 --- a/tools/maps/osm_vector.py +++ b/tools/maps/osm_vector.py @@ -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)