Rendering curved flight paths from geojson file

Question

I'm working on a Laravel based application where users can create flights and then view the flight routes on maps. Some of these maps consist of up to 10.000 flights. With these numbers Mapbox would sometimes crash the browser or take a really long time to render, especially since we were using arc.js to calculate the Great Circle routes.

We then switched to creating the flight routes as a single geojson file that was then loaded directly by Mapbox. The maps now load really fast (1-2 seconds instead of up to half a minute), but the routes are straight and do not cross the dateline, which is a real problem as that's not how planes fly.

First I looked for some type of setting within Mapbox that would allow me to render lines as great circles, but I couldn't find anything. Next I looked for PHP libraries similar to arc.js to output the routes to the geojson file, but while there are plenty of libraries that will calculate the distance based on a great circle, I didn't find any that actually produce a route. Currently, I'm looking at the database level. We already use PostGIS, so I thought that there might be a way to use it to calculate the route. So far I have this here, cobbled together from various sources, but it's still throwing errors; ST_MakeLine does not exist...:

SELECT ST_Transform(
     ST_Segmentize(
         ST_MakeLine(
           dep.point,
           arr.point
         )::geography,
         100000
     )::geometry,
     3857
   ) as the_geom_webmercator
FROM flights f
LEFT JOIN airports dep ON f.departure_airport_id = dep.id
LEFT JOIN airports arr ON f.arrival_airport_id = arr.id;

I'm kinda hoping that there's a hidden way of displaying curved lines directly in Mapbox, maybe via a source/layer? Failing that, I'd be happy for any pointers for PHP libraries or even PostGIS resources.

Cheers!

Answer 1

So, I ended up going the PHP library route. It might not be as performant as other options, but it was the easiest. In case someone else is wondering how to achieve this, here is my solution using the phpgeo package:

namespace App\Services;

use Location\Bearing\BearingEllipsoidal;
use Location\Coordinate;
use Location\Distance\Vincenty;

class Arc
{
    protected Coordinate $start;
    protected Coordinate $end;

    /**
     * @param  \Location\Coordinate  $start
     * @param  \Location\Coordinate  $end
     */
    public function __construct(Coordinate $start, Coordinate $end)
    {
        $this->start = $start;
        $this->end   = $end;
    }

    /**
     * @param  int  $points
     * @return array
     */
    public function line(int $points = 100): array
    {
        $bearingCalculator  = new BearingEllipsoidal();
        $distanceCalculator = new Vincenty();

        $totalDistance    = $distanceCalculator->getDistance($this->start, $this->end);
        $intervalDistance = $totalDistance / ($points + 1);

        $currentBearing = $bearingCalculator->calculateBearing($this->start, $this->end);
        $currentCoords  = $this->start;

        $polyline = [
            [
                $this->start->getLng(),
                $this->start->getLat(),
            ],
        ];

        for ($i = 1; $i < ($points + 1); $i++) {
            $currentCoords = $bearingCalculator->calculateDestination(
                $currentCoords,
                $currentBearing,
                $intervalDistance
            );

            $point = [
                $currentCoords->getLng(),
                $currentCoords->getLat(),
            ];

            $polyline[$i] = $this->forAntiMeridian(
                $polyline[$i - 1],
                $point
            );

            $currentBearing = $bearingCalculator->calculateBearing(
                $currentCoords,
                $this->end
            );
        }

        $polyline[] = $this->forAntiMeridian(
            $polyline[$i - 1],
            [
                $this->end->getLng(),
                $this->end->getLat(),
            ]
        );

        return array_values($polyline);
    }

    /**
     * @param  array  $start
     * @param  array  $end
     * @return array
     */
    protected function forAntiMeridian(array $start, array $end): array
    {
        $startLng = $start[0];
        $endLng   = $end[0];

        if ($endLng - $startLng > 180) {
            $end[0] -= 360;
        } elseif ($startLng - $endLng > 180) {
            $end[0] += 360;
        }

        return $end;
    }
}

I found the contents of the line method, which does the bulk of the work, somewhere else, but can't seem to find it again. I just freshened it up a bit and added calculations for when a line crosses the anti meridian. Mapbox allows coordinates that are out of bounds for just this reason.

Answer 2

In addition to your php based solution:

If you're getting the error ST_MakeLine does not exist you are either passing the wrong parameters to the function or the PostGIS extension isn't installed at all. That being said, to get the curved lines you're looking for, you can transform your points to a spherical projected reference system before ST_Segmentize. After that, in order to get a single GeoJSON String, you might wanna ST_Collect the rows and then serialise it with ST_AsGeoJSON.

WITH airports (point) AS (
  VALUES 
    ('SRID=4326;POINT(-120.26 35.38)'::geometry),
    ('SRID=4326;POINT(-98.24 25.74)'::geometry),
    ('SRID=4326;POINT(-36.01 -8.23)'::geometry),
    ('SRID=4326;POINT(131.13 -26.11)'::geometry),
    ('SRID=4326;POINT(64.03 76.25)'::geometry),
    ('SRID=4326;POINT(30.99 27.46)'::geometry)
)
SELECT 
 ST_AsGeoJSON(
  ST_Collect(
    ST_Transform(
      ST_Segmentize(
        ST_MakeLine(
          -- Transforming from 4326 to 53027
          ST_Transform('SRID=4326;POINT(12.62 52.95)'::geometry, 53027),
          ST_Transform(arr.point, 53027)),
        100000
     ),4326))) -- Here you can transform it to the SRS of your choice, e.g. WGS84
FROM airports arr;

Demo: db<>fiddle