从 geojson 文件渲染弯曲的飞行路径

Rendering curved flight paths from geojson file

我正在开发一个基于 Laravel 的应用程序,用户可以在其中创建航班,然后在地图上查看飞行路线。其中一些地图包含多达 10.000 个航班。使用这些数字,Mapbox 有时会使浏览器崩溃或需要很长时间才能呈现,尤其是因为我们使用 arc.js 来计算大圆路线。

然后我们转而将飞行路线创建为单个 geojson 文件,然后由 Mapbox 直接加载。现在地图加载速度非常快(1-2 秒,而不是最多半分钟),但路线是笔直的,不会越过日期变更线,这是一个真正的问题,因为飞机不是这样飞行的。

首先,我在 Mapbox 中寻找某种类型的设置,允许我将线条渲染为大圆圈,但我找不到任何东西。接下来我寻找 PHP 类似于 arc.js 的库来将路线输出到 geojson 文件,但是虽然有很多库会根据大圆来计算距离,但我没有找到任何库实际上产生了一条路线。目前,我正在研究数据库级别。我们已经在使用 PostGIS,所以我想可能有一种方法可以用它来计算路线。到目前为止,我这里有这个,从各种来源拼凑而成,但它仍然会抛出错误; ST_MakeLine 不存在...:[=​​12=]

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;

我有点希望有一种隐藏的方式可以直接在 Mapbox 中显示曲线,也许是通过 source/layer?如果做不到这一点,我很乐意为 PHP 库甚至 PostGIS 资源提供任何指示。

干杯!

所以,我最终选择了 PHP 图书馆路线。它的性能可能不如其他选项,但它是最简单的。如果有人想知道如何实现这一点,这里是我使用 phpgeo 包的解决方案:

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;
    }
}

我在其他地方找到了 line 方法的内容,它完成了大部分工作,但似乎无法再次找到它。我只是稍微更新了一下,并添加了一条线何时穿过反子午线的计算。出于这个原因,Mapbox 允许坐标越界。

除了基于 php 的解决方案:

如果您遇到错误 ST_MakeLine does not exist,您要么向函数传递了错误的参数,要么根本没有安装 PostGIS 扩展。也就是说,要获得您正在寻找的曲线,您可以将您的点转换为 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 .

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;

演示:db<>fiddle