stammtisch-hersbruck/frontend/src/lib/gpx.ts

/**
 * GPX-Verarbeitung im Browser.
 *
 * Das Backend bleibt reines PocketBase ohne Hooks: Beim Upload liest das
 * Frontend die Datei, rechnet Länge, Höhenmeter und Bounding-Box aus und
 * schickt Originaldatei plus abgeleitete Werte gemeinsam an die API.
 */

export type TrackPoint = { lat: number; lng: number; ele: number | null }

export type ParsedGpx = {
    points: TrackPoint[]
    geojson: { type: 'LineString'; coordinates: [number, number][] }
    /** Kumulative Distanz je Koordinate in `geojson.coordinates`, in Metern. */
    coordDistances: number[]
    elevation: { d: number; ele: number }[]
    bounds: [[number, number], [number, number]]
    distance_m: number
    ascent_m: number
}

/** Erdradius in Metern */
const R = 6_371_000

/**
 * Anstiege unterhalb dieser Schwelle gelten als GPS-Rauschen. Ohne sie
 * summieren sich Messfehler zu absurden Höhenmetern — ein flacher Waldweg
 * käme leicht auf mehrere hundert Meter.
 */
const ELE_NOISE_M = 3

/** Distanz zweier Punkte in Metern (Haversine). */
export function haversine(a: TrackPoint, b: TrackPoint): number {
    const toRad = (d: number) => (d * Math.PI) / 180
    const dLat = toRad(b.lat - a.lat)
    const dLng = toRad(b.lng - a.lng)
    const lat1 = toRad(a.lat)
    const lat2 = toRad(b.lat)

    const h =
        Math.sin(dLat / 2) ** 2 +
        Math.sin(dLng / 2) ** 2 * Math.cos(lat1) * Math.cos(lat2)

    return 2 * R * Math.asin(Math.sqrt(h))
}

/**
 * Douglas-Peucker: entfernt Punkte, die kaum von der Verbindungslinie ihrer
 * Nachbarn abweichen. Die Toleranz wird verdoppelt, bis das Ergebnis unter
 * maxPoints liegt — so bleibt die Form erhalten, während die Datenmenge
 * beherrschbar wird.
 */
export function simplify(points: TrackPoint[], maxPoints = 2000): TrackPoint[] {
    return simplifyIndices(points, maxPoints).map((i) => points[i])
}

/**
 * Wie `simplify`, liefert aber die Indizes der erhaltenen Punkte statt der
 * Punkte selbst. Damit lässt sich zu jeder vereinfachten Koordinate die
 * ursprüngliche Distanz nachschlagen.
 */
export function simplifyIndices(points: TrackPoint[], maxPoints = 2000): number[] {
    if (points.length <= maxPoints) {
        return points.map((_, i) => i)
    }

    let tolerance = 0.00001
    let out = points.map((_, i) => i)

    for (let i = 0; i < 30 && out.length > maxPoints; i++) {
        out = douglasPeuckerIndices(points, points.map((_, k) => k), tolerance)
        tolerance *= 2
    }

    return out
}

function douglasPeuckerIndices(all: TrackPoint[], idx: number[], tolerance: number): number[] {
    if (idx.length < 3) return idx

    let maxDist = 0
    let split = 0

    for (let i = 1; i < idx.length - 1; i++) {
        const d = perpendicularDistance(all[idx[i]], all[idx[0]], all[idx[idx.length - 1]])
        if (d > maxDist) {
            maxDist = d
            split = i
        }
    }

    if (maxDist <= tolerance) return [idx[0], idx[idx.length - 1]]

    const left = douglasPeuckerIndices(all, idx.slice(0, split + 1), tolerance)
    const right = douglasPeuckerIndices(all, idx.slice(split), tolerance)
    return [...left.slice(0, -1), ...right]
}

/**
 * Abstand eines Punktes zur Geraden durch start und end — in Grad, nicht in
 * Metern. Für den Vergleich innerhalb eines Tracks genügt das und spart die
 * teure Projektion.
 */
function perpendicularDistance(p: TrackPoint, start: TrackPoint, end: TrackPoint): number {
    const dx = end.lng - start.lng
    const dy = end.lat - start.lat

    if (dx === 0 && dy === 0) {
        return Math.hypot(p.lng - start.lng, p.lat - start.lat)
    }

    const t =
        ((p.lng - start.lng) * dx + (p.lat - start.lat) * dy) / (dx * dx + dy * dy)
    const clamped = Math.max(0, Math.min(1, t))

    return Math.hypot(
        p.lng - (start.lng + clamped * dx),
        p.lat - (start.lat + clamped * dy),
    )
}

/**
 * Liest eine GPX-Datei. Wirft mit einer für Nutzer verständlichen Meldung,
 * wenn die Datei kein gültiges XML ist oder keine Punkte enthält.
 */
export function parseGpx(xml: string): ParsedGpx {
    const doc = new DOMParser().parseFromString(xml, 'application/xml')

    if (doc.querySelector('parsererror')) {
        throw new Error('Die Datei ist kein gültiges XML.')
    }

    // Manche Programme exportieren Routen (rtept) statt Tracks (trkpt).
    let nodes = Array.from(doc.getElementsByTagName('trkpt'))
    if (nodes.length === 0) {
        nodes = Array.from(doc.getElementsByTagName('rtept'))
    }

    if (nodes.length === 0) {
        throw new Error('Die GPX-Datei enthält keine Punkte.')
    }

    const points: TrackPoint[] = []
    for (const n of nodes) {
        const latAttr = n.getAttribute('lat')
        const lngAttr = n.getAttribute('lon')
        // Number(null) ist 0, nicht NaN — ohne diese Prüfung landen Punkte
        // ohne Koordinaten bei 0/0 im Golf von Guinea und verfälschen
        // Streckenlänge und Bounding-Box.
        if (latAttr === null || lngAttr === null) continue

        const lat = Number(latAttr)
        const lng = Number(lngAttr)
        if (!Number.isFinite(lat) || !Number.isFinite(lng)) continue

        const eleText = n.getElementsByTagName('ele')[0]?.textContent
        const ele = eleText != null && eleText.trim() !== '' ? Number(eleText) : null

        points.push({ lat, lng, ele: Number.isFinite(ele as number) ? ele : null })
    }

    if (points.length === 0) {
        throw new Error('Die GPX-Datei enthält keine Punkte.')
    }

    // Länge, Höhenprofil und Höhenmeter in einem Durchlauf
    let distance = 0
    let ascent = 0
    const elevation: { d: number; ele: number }[] = []
    let lastCountedEle: number | null = points[0].ele

    // Kumulative Distanz je Originalpunkt — Brücke zwischen elevation (ein
    // Eintrag pro Originalpunkt) und der vereinfachten Geometrie.
    const distAt: number[] = new Array(points.length)
    distAt[0] = 0

    if (points[0].ele != null) elevation.push({ d: 0, ele: points[0].ele })

    for (let i = 1; i < points.length; i++) {
        distance += haversine(points[i - 1], points[i])
        distAt[i] = distance

        const ele = points[i].ele
        if (ele == null) continue

        elevation.push({ d: Math.round(distance), ele })

        if (lastCountedEle == null) {
            lastCountedEle = ele
        } else {
            const delta = ele - lastCountedEle
            if (delta >= ELE_NOISE_M) {
                ascent += delta
                lastCountedEle = ele
            } else if (delta <= -ELE_NOISE_M) {
                // Abstieg über Schwelle: neuer Bezugspunkt, aber nicht gezählt
                lastCountedEle = ele
            }
            // Kleine Schwankungen (unter 3 m) ignorieren, Bezugspunkt bleibt
        }
    }

    const lats = points.map((p) => p.lat)
    const lngs = points.map((p) => p.lng)

    const keep = simplifyIndices(points)
    const simplified = keep.map((i) => points[i])

    return {
        points: simplified,
        geojson: {
            type: 'LineString',
            coordinates: simplified.map((p) => [p.lng, p.lat] as [number, number]),
        },
        coordDistances: keep.map((i) => Math.round(distAt[i])),
        elevation,
        bounds: [
            [Math.min(...lngs), Math.min(...lats)],
            [Math.max(...lngs), Math.max(...lats)],
        ],
        distance_m: Math.round(distance),
        ascent_m: Math.round(ascent),
    }
}

/** Ein auf die Trail-Linie projizierter Punkt. */
export type LinePosition = {
    lng: number
    lat: number
    /** Index des Segmentanfangs in `coordinates` */
    index: number
    /** Lage im Segment zwischen 0 und 1 */
    t: number
}

/**
 * Projiziert einen beliebigen Punkt auf die nächstgelegene Stelle der Linie.
 * Damit landet ein Marker exakt auf dem Track, egal wie genau geklickt wurde.
 *
 * Gerechnet wird in Grad, die Längengrade aber mit cos(lat) gestaucht — sonst
 * wäre ein Grad Ost-West bei uns fast doppelt so „lang" wie einer Nord-Süd und
 * die Projektion würde spürbar daneben liegen. Für die paar Meter Suchradius
 * genügt das; eine echte Projektion wäre hier Aufwand ohne Gewinn.
 */
export function nearestOnLine(
    coords: [number, number][],
    lng: number,
    lat: number,
): LinePosition | null {
    if (!coords.length) return null
    if (coords.length === 1) {
        return { lng: coords[0][0], lat: coords[0][1], index: 0, t: 0 }
    }

    const kx = Math.cos((lat * Math.PI) / 180) || 1

    let best: LinePosition | null = null
    let bestDist = Infinity

    for (let i = 0; i < coords.length - 1; i++) {
        const [ax, ay] = coords[i]
        const [bx, by] = coords[i + 1]

        const dx = (bx - ax) * kx
        const dy = by - ay
        const px = (lng - ax) * kx
        const py = lat - ay

        const len = dx * dx + dy * dy
        // Doppelter Punkt im Track: das Segment hat keine Richtung.
        const t = len === 0 ? 0 : Math.max(0, Math.min(1, (px * dx + py * dy) / len))

        const ox = px - t * dx
        const oy = py - t * dy
        const dist = ox * ox + oy * oy

        if (dist < bestDist) {
            bestDist = dist
            best = { lng: ax + (bx - ax) * t, lat: ay + (by - ay) * t, index: i, t }
        }
    }

    return best
}

/**
 * Streckenkilometer einer projizierten Stelle. `coordDistances` hat einen
 * Eintrag je Koordinate; zwischen zwei Stützpunkten wird linear interpoliert.
 * Ohne Distanztabelle (Versionen von vor dem Feld) gibt es kein Ergebnis.
 */
export function distanceAlong(pos: LinePosition, coordDistances: number[]): number | null {
    if (coordDistances.length <= pos.index) return null

    const a = coordDistances[pos.index]
    const b = coordDistances[pos.index + 1]
    if (b == null) return a

    return a + (b - a) * pos.t
}

/**
 * Gegenrichtung zu `distanceAlong`: die Koordinate zu einem Streckenkilometer.
 * Ebenfalls interpoliert, damit ein Punkt beim Überfahren gleichmäßig wandert
 * statt von Stützpunkt zu Stützpunkt zu springen.
 */
export function positionAtDistance(
    coords: [number, number][],
    coordDistances: number[],
    distance: number,
): [number, number] | null {
    if (!coords.length || coordDistances.length < 2) return null

    const last = Math.min(coords.length, coordDistances.length) - 1
    if (distance <= coordDistances[0]) return coords[0]
    if (distance >= coordDistances[last]) return coords[last]

    // Binäre Suche nach dem Segment, in dem die Distanz liegt.
    let lo = 0
    let hi = last
    while (hi - lo > 1) {
        const mid = (lo + hi) >> 1
        if (coordDistances[mid] <= distance) lo = mid
        else hi = mid
    }

    const span = coordDistances[hi] - coordDistances[lo]
    const t = span === 0 ? 0 : (distance - coordDistances[lo]) / span

    return [
        coords[lo][0] + (coords[hi][0] - coords[lo][0]) * t,
        coords[lo][1] + (coords[hi][1] - coords[lo][1]) * t,
    ]
}