|
| 1 | +"""Search the OpenStreetMap public GPS trace archive, and turn one trace into KML. |
| 2 | +
|
| 3 | + python3 gpstraces.py 47.815,7.940,47.832,7.960 # what has been ridden here |
| 4 | + python3 gpstraces.py 47.815,7.940,47.832,7.960 --near 7.95,47.82 # sort by distance |
| 5 | + python3 gpstraces.py --trace 12004364 # profile one trace |
| 6 | + python3 gpstraces.py --trace 12004364 --emit # <coordinates> block |
| 7 | +
|
| 8 | +This is the fallback for when OSM has no mapped trail geometry: real riders' |
| 9 | +GPS logs. It is a different database from the one `trails.py` queries - traces |
| 10 | +are uploaded files, not tagged map objects, and nobody has curated them. |
| 11 | +
|
| 12 | +Two endpoints, and the difference between them matters: |
| 13 | +
|
| 14 | + api/0.6/trackpoints?bbox=W,S,E,N&page=N (the listing form, used here) |
| 15 | + Public, no auth. 5000 points per page, so page until a page repeats. |
| 16 | + Returns lat/lon/time only - the elevation is stripped. Traces whose |
| 17 | + owner chose "identifiable" or "trackable" carry <name> and <url>; |
| 18 | + the rest arrive as anonymous <trk> blocks with no owner and no id. |
| 19 | +
|
| 20 | + /trace/<id>/data (the download form, used by --trace) |
| 21 | + Also public and unauthenticated, on www.openstreetmap.org rather than |
| 22 | + api.openstreetmap.org. Returns the ORIGINAL uploaded file: full |
| 23 | + precision, full point count, and <ele> on every point. This is the one |
| 24 | + you actually build a line from. |
| 25 | +
|
| 26 | + api/0.6/gpx/<id>/data is the documented equivalent and returns 401 |
| 27 | + without OAuth. Use the /trace/ URL. |
| 28 | +
|
| 29 | +Bounding boxes are given here as S,W,N,E, like the other scripts. The API |
| 30 | +wants W,S,E,N; this converts. The API caps a bbox at 0.25 square degrees. |
| 31 | +
|
| 32 | +Licence: traces are ODbL like the rest of OSM. Credit them as |
| 33 | +"Geometry simplified from OpenStreetMap GPS trace <id> (ODbL)". |
| 34 | +""" |
| 35 | + |
| 36 | +import bz2 |
| 37 | +import gzip |
| 38 | +import io |
| 39 | +import math |
| 40 | +import sys |
| 41 | +import urllib.request |
| 42 | +import xml.etree.ElementTree as ET |
| 43 | +import zipfile |
| 44 | + |
| 45 | +from trails import length, simplify |
| 46 | + |
| 47 | +UA = "alpine-mtb-map/1.0 (+https://github.com/vemonet/alpine-mtb-map)" |
| 48 | +BBOX_URL = "https://api.openstreetmap.org/api/0.6/trackpoints?bbox={w},{s},{e},{n}&page={p}" |
| 49 | +TRACE_URL = "https://www.openstreetmap.org/trace/{id}/data" |
| 50 | +GPX0 = "{http://www.topografix.com/GPX/1/0}" |
| 51 | + |
| 52 | + |
| 53 | +def get(url): |
| 54 | + req = urllib.request.Request(url, headers={"User-Agent": UA}) |
| 55 | + with urllib.request.urlopen(req, timeout=60) as r: |
| 56 | + return r.read() |
| 57 | + |
| 58 | + |
| 59 | +def scan(bbox, max_pages=20): |
| 60 | + """Page the bbox endpoint until it stops giving new data. |
| 61 | +
|
| 62 | + Each <trk> element is kept separately rather than merged by name: one |
| 63 | + upload can cross the bbox several times, and merging the visits draws a |
| 64 | + straight bar between them. Anonymous traces have no url and cannot be |
| 65 | + downloaded in full - their bbox points are all you will ever get. |
| 66 | + """ |
| 67 | + s, w, n, e = (float(x) for x in bbox.split(",")) |
| 68 | + out, seen = [], set() |
| 69 | + for p in range(max_pages): |
| 70 | + body = get(BBOX_URL.format(w=w, s=s, e=e, n=n, p=p)) |
| 71 | + if body in seen: |
| 72 | + break |
| 73 | + seen.add(body) |
| 74 | + trks = ET.fromstring(body).findall(GPX0 + "trk") |
| 75 | + if not trks: |
| 76 | + break |
| 77 | + for trk in trks: |
| 78 | + pts = [ |
| 79 | + {"lat": float(t.get("lat")), "lon": float(t.get("lon"))} |
| 80 | + for t in trk.iter(GPX0 + "trkpt") |
| 81 | + ] |
| 82 | + if len(pts) < 2: |
| 83 | + continue |
| 84 | + url = trk.findtext(GPX0 + "url") or "" |
| 85 | + # No timestamps means this is the unordered pool, not a path. The |
| 86 | + # API dumps every trace marked merely "public" into anonymous |
| 87 | + # blocks of 5000 points in no order at all, at the tail of the |
| 88 | + # pagination. Chain them and you get a 1800 km scribble. |
| 89 | + ordered = trk.find(f".//{GPX0}trkpt/{GPX0}time") is not None |
| 90 | + out.append( |
| 91 | + { |
| 92 | + "name": trk.findtext(GPX0 + "name") or "(anonymous)", |
| 93 | + "id": url.rsplit("/", 1)[-1] if url else "", |
| 94 | + "user": url.split("/user/")[-1].split("/")[0] if url else "", |
| 95 | + "pts": pts, |
| 96 | + "ordered": ordered, |
| 97 | + "km": length(pts) / 1000 if ordered else 0.0, |
| 98 | + } |
| 99 | + ) |
| 100 | + return out |
| 101 | + |
| 102 | + |
| 103 | +def decompress(body): |
| 104 | + """/trace/<id>/data serves the file exactly as uploaded, not as GPX. |
| 105 | +
|
| 106 | + Most of the archive predates the web upload form, so a large share of it is |
| 107 | + .gpx.bz2, .gpx.gz or .zip. At Les Gets 9 of 16 downloads were bzip2. Feed |
| 108 | + those to an XML parser and you get "not well-formed: line 1, column 7", |
| 109 | + which reads like a corrupt trace rather than a compressed one. |
| 110 | + """ |
| 111 | + if body[:3] == b"BZh": |
| 112 | + return bz2.decompress(body) |
| 113 | + if body[:2] == b"\x1f\x8b": |
| 114 | + return gzip.decompress(body) |
| 115 | + if body[:2] == b"PK": |
| 116 | + with zipfile.ZipFile(io.BytesIO(body)) as z: |
| 117 | + names = [n for n in z.namelist() if n.lower().endswith(".gpx")] |
| 118 | + return z.read(names[0] if names else z.namelist()[0]) |
| 119 | + return body |
| 120 | + |
| 121 | + |
| 122 | +def parse_full(body): |
| 123 | + """Read a downloaded trace. GPX 1.0 and 1.1 differ only in namespace.""" |
| 124 | + root = ET.fromstring(decompress(body).lstrip(b"\xef\xbb\xbf")) |
| 125 | + ns = root.tag[: root.tag.index("}") + 1] |
| 126 | + pts = [] |
| 127 | + for t in root.iter(ns + "trkpt"): |
| 128 | + ele = t.findtext(ns + "ele") |
| 129 | + pts.append( |
| 130 | + { |
| 131 | + "lat": float(t.get("lat")), |
| 132 | + "lon": float(t.get("lon")), |
| 133 | + "ele": float(ele) if ele else None, |
| 134 | + } |
| 135 | + ) |
| 136 | + return pts |
| 137 | + |
| 138 | + |
| 139 | +def main(argv): |
| 140 | + if "--trace" in argv: |
| 141 | + tid = argv[argv.index("--trace") + 1] |
| 142 | + pts = parse_full(get(TRACE_URL.format(id=tid))) |
| 143 | + ele = [p["ele"] for p in pts if p["ele"] is not None] |
| 144 | + km = length(pts) / 1000 |
| 145 | + if ele: |
| 146 | + drop = ele[0] - ele[-1] |
| 147 | + print( |
| 148 | + f"trace {tid}: {km:.1f} km, {len(pts)} pts, " |
| 149 | + f"{max(ele):.0f} -> {min(ele):.0f} m, net drop {drop:.0f} m" |
| 150 | + f"{' (STORED UPHILL - reverse it)' if drop < 0 else ''}" |
| 151 | + ) |
| 152 | + else: |
| 153 | + print(f"trace {tid}: {km:.1f} km, {len(pts)} pts, no elevation in the file") |
| 154 | + if "--emit" in argv: |
| 155 | + simp = simplify(pts) |
| 156 | + print( |
| 157 | + f"\n<!-- GPS trace {tid}: {km:.1f} km, {len(pts)} pts -> {len(simp)}. " |
| 158 | + f"Source: OSM GPS trace {tid} (ODbL) -->" |
| 159 | + ) |
| 160 | + coords = " ".join(f"{p['lon']:.6f},{p['lat']:.6f},0" for p in simp) |
| 161 | + print(f"<coordinates>{coords}</coordinates>") |
| 162 | + return |
| 163 | + |
| 164 | + bbox = argv[1] |
| 165 | + near = None |
| 166 | + if "--near" in argv: |
| 167 | + lon, lat = (float(x) for x in argv[argv.index("--near") + 1].split(",")) |
| 168 | + near = {"lon": lon, "lat": lat} |
| 169 | + trks = scan(bbox) |
| 170 | + if not trks: |
| 171 | + print("no public GPS traces in this bbox") |
| 172 | + return |
| 173 | + |
| 174 | + def hav(a, b): |
| 175 | + r = 6371000 |
| 176 | + p1, p2 = math.radians(a["lat"]), math.radians(b["lat"]) |
| 177 | + h = ( |
| 178 | + math.sin((p2 - p1) / 2) ** 2 |
| 179 | + + math.cos(p1) * math.cos(p2) * math.sin(math.radians(b["lon"] - a["lon"]) / 2) ** 2 |
| 180 | + ) |
| 181 | + return 2 * r * math.asin(math.sqrt(h)) |
| 182 | + |
| 183 | + if near: |
| 184 | + for t in trks: |
| 185 | + t["d"] = min(hav(near, p) for p in t["pts"]) |
| 186 | + trks.sort(key=lambda t: t["d"]) |
| 187 | + else: |
| 188 | + trks.sort(key=lambda t: -t["km"]) |
| 189 | + |
| 190 | + named = sum(1 for t in trks if t["id"]) |
| 191 | + pool = sum(1 for t in trks if not t["ordered"]) |
| 192 | + print( |
| 193 | + f"{len(trks)} track segments, {named} downloadable, " |
| 194 | + f"{len(trks) - named - pool} anonymous but ordered, {pool} unordered pool (unusable)\n" |
| 195 | + ) |
| 196 | + for t in trks[:40]: |
| 197 | + if not t["ordered"]: |
| 198 | + continue |
| 199 | + d = f"{t['d']:5.0f} m " if near else "" |
| 200 | + who = f"{t['id']:>10} {t['user'][:16]:16s}" if t["id"] else f"{'-':>10} {'':16s}" |
| 201 | + print(f"{d}{t['km']:7.2f} km {len(t['pts']):6d} pts {who} {t['name'][:44]}") |
| 202 | + |
| 203 | + |
| 204 | +if __name__ == "__main__": |
| 205 | + if len(sys.argv) < 2: |
| 206 | + sys.exit(__doc__) |
| 207 | + main(sys.argv) |
0 commit comments