NavDataExplorer/geojson.js

import fs from "fs";
import geo from "geojson";
import geolib from "geolib";

const π = Math.PI;

Number.prototype.toRadians = function () {
  return (this * Math.PI) / 180;
};
Number.prototype.toDegrees = function () {
  return (this * 180) / Math.PI;
};

const nmiToMetre = (nmi) => nmi * 1852.0;

const intersection = (p1, brng1, p2, brng2) => {
  if (isNaN(brng1)) throw new TypeError(`invalid brng1 ${brng1}`);
  if (isNaN(brng2)) throw new TypeError(`invalid brng2 ${brng2}`);

  // see www.edwilliams.org/avform.htm#Intersection

  const φ1 = p1.lat.toRadians(),
    λ1 = p1.lon.toRadians();
  const φ2 = p2.lat.toRadians(),
    λ2 = p2.lon.toRadians();
  const θ13 = Number(brng1).toRadians(),
    θ23 = Number(brng2).toRadians();
  const Δφ = φ2 - φ1,
    Δλ = λ2 - λ1;

  // angular distance p1-p2
  const δ12 =
    2 *
    Math.asin(
      Math.sqrt(
        Math.sin(Δφ / 2) * Math.sin(Δφ / 2) +
          Math.cos(φ1) * Math.cos(φ2) * Math.sin(Δλ / 2) * Math.sin(Δλ / 2)
      )
    );
  if (Math.abs(δ12) < Number.EPSILON) return p1; // coincident points

  // initial/final bearings between points
  const cosθa =
    (Math.sin(φ2) - Math.sin(φ1) * Math.cos(δ12)) /
    (Math.sin(δ12) * Math.cos(φ1));
  const cosθb =
    (Math.sin(φ1) - Math.sin(φ2) * Math.cos(δ12)) /
    (Math.sin(δ12) * Math.cos(φ2));
  const θa = Math.acos(Math.min(Math.max(cosθa, -1), 1)); // protect against rounding errors
  const θb = Math.acos(Math.min(Math.max(cosθb, -1), 1)); // protect against rounding errors

  const θ12 = Math.sin(λ2 - λ1) > 0 ? θa : 2 * π - θa;
  const θ21 = Math.sin(λ2 - λ1) > 0 ? 2 * π - θb : θb;

  const α1 = θ13 - θ12; // angle 2-1-3
  const α2 = θ21 - θ23; // angle 1-2-3

  if (Math.sin(α1) == 0 && Math.sin(α2) == 0) return null; // infinite intersections
  if (Math.sin(α1) * Math.sin(α2) < 0) return null; // ambiguous intersection (antipodal/360°)

  const cosα3 =
    -Math.cos(α1) * Math.cos(α2) + Math.sin(α1) * Math.sin(α2) * Math.cos(δ12);

  const δ13 = Math.atan2(
    Math.sin(δ12) * Math.sin(α1) * Math.sin(α2),
    Math.cos(α2) + Math.cos(α1) * cosα3
  );

  const φ3 = Math.asin(
    Math.min(
      Math.max(
        Math.sin(φ1) * Math.cos(δ13) +
          Math.cos(φ1) * Math.sin(δ13) * Math.cos(θ13),
        -1
      ),
      1
    )
  );

  const Δλ13 = Math.atan2(
    Math.sin(θ13) * Math.sin(δ13) * Math.cos(φ1),
    Math.cos(δ13) - Math.sin(φ1) * Math.sin(φ3)
  );
  const λ3 = λ1 + Δλ13;

  const lat = φ3.toDegrees();
  const lon = λ3.toDegrees();

  return { lat, lon, name: "INTC" };
};

const SPEED = 4 / 60; //nmi/s
const TURN = 3; //deg/s
const VS = 1; //fps

// ADJ to location, (NG), E is neg.
const MAGVAR = -1;

const PROCEDURE_ID = 10394;

const waypoints = JSON.parse(fs.readFileSync("Waypoints.json"));

const terminal = JSON.parse(fs.readFileSync("Terminals.json")).filter(
  ({ ID }) => ID === PROCEDURE_ID
)[0];
const runway = JSON.parse(fs.readFileSync("Runways.json")).filter(
  ({ ID }) => ID === terminal.RwyID
)[0];

const procedure = JSON.parse(fs.readFileSync(`TermID_${PROCEDURE_ID}.json`));

const points = [{ lat: runway.Latitude, lon: runway.Longitude }];
const lines = [];

procedure.forEach((leg) => {
  const waypoint = waypoints.filter(({ ID }) => ID === leg.WptID)[0];

  switch (leg.TrackCode) {
    case "AF":
    case "CA":
    case "CD":
      break;
    case "CF": {
      points.push({
        lat: leg.WptLat,
        lon: leg.WptLon,
        name: waypoint?.Ident ?? undefined,
      });

      if (leg.TurnDir) {
        const line = [
          [points[points.length - 2].lon, points[points.length - 2].lat],
        ];

        let inbdCrs = Math.round(
          geolib.getGreatCircleBearing(
            {
              latitude: points[points.length - 3].lat,
              longitude: points[points.length - 3].lon,
            },
            {
              latitude: points[points.length - 2].lat,
              longitude: points[points.length - 2].lon,
            }
          )
        );
        let outbCrs = Math.round(leg.Course - MAGVAR);
        outbCrs =
          outbCrs >= 360
            ? outbCrs - 360
            : outbCrs < 0
            ? outbCrs + 360
            : outbCrs;

        const delta = Math.abs(inbdCrs - 180 - outbCrs);

        let perpCrs;
        if (leg.TurnDir === "R") {
          perpCrs = inbdCrs + 90 + delta / 2;
          perpCrs = perpCrs >= 360 ? perpCrs - 360 : perpCrs;
        } else {
          perpCrs = inbdCrs - 90 - delta / 2;
          perpCrs = perpCrs < 0 ? perpCrs + 360 : perpCrs;
        }

        const perpFix = intersection(
          points[points.length - 2],
          Math.round(perpCrs - MAGVAR),
          points[points.length - 1],
          Math.round(leg.Course + 180 - MAGVAR)
        );

        const midFix = geolib.getCenter([
          {
            latitude: points[points.length - 2].lat,
            longitude: points[points.length - 2].lon,
          },
          {
            latitude: perpFix.lat,
            longitude: perpFix.lon,
          },
        ]);
        const dist = geolib.getDistance(
          {
            latitude: points[points.length - 2].lat,
            longitude: points[points.length - 2].lon,
          },
          {
            latitude: perpFix.lat,
            longitude: perpFix.lon,
          }
        );
        const mid1 = geolib.computeDestinationPoint(midFix, dist / 2, inbdCrs);

        let invCrs = outbCrs + 180;
        invCrs = invCrs >= 360 ? invCrs - 360 : invCrs;
        const mid2 = geolib.computeDestinationPoint(midFix, dist / 2, invCrs);

        line.push([mid1.longitude, mid1.latitude]);
        line.push([mid2.longitude, mid2.latitude]);
        line.push([perpFix.lon, perpFix.lat]);
        line.push([
          points[points.length - 1].lon,
          points[points.length - 1].lat,
        ]);

        lines.push({ line });
      } else {
        lines.push({
          line: [
            [points[points.length - 2].lon, points[points.length - 2].lat],
            [points[points.length - 1].lon, points[points.length - 1].lat],
          ],
        });
      }
      break;
    }
    case "CI":
      break;
    case "CR": {
      const intc = intersection(
        points[points.length - 1],
        Math.round(leg.Course - MAGVAR),
        { lat: leg.NavLat, lon: leg.NavLon },
        Math.round(leg.NavBear - MAGVAR)
      );

      points.push(intc);
      lines.push({
        line: [
          [points[points.length - 2].lon, points[points.length - 2].lat],
          [points[points.length - 1].lon, points[points.length - 1].lat],
        ],
      });
    }
    case "DF":
    case "FA":
    case "FC":
    case "FD":
    case "FM":
    case "HA":
    case "HF":
    case "HM":
    case "IF":
    case "PI":
    case "RF":
    case "TF":
    case "VA":
    case "VD":
    case "VI":
      break;
    case "VM": {
      const end = geolib.computeDestinationPoint(
        {
          latitude: points[points.length - 1].lat,
          longitude: points[points.length - 1].lon,
        },
        nmiToMetre(10),
        leg.Course - MAGVAR
      );
      lines.push({
        line: [
          [points[points.length - 1].lon, points[points.length - 1].lat],
          [end.longitude, end.latitude],
        ],
      });
    }
    case "VR":
    case "AF":
    default:
      break;
  }
});

const output = geo.parse([...points, ...lines], {
  Point: ["lat", "lon"],
  LineString: "line",
});

console.log(JSON.stringify(output, null, 2));