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

问题描述

我正在开发基于 Laravel 的应用程序，用户可以在其中创建航班，然后在地图上查看航班路线。其中一些地图包含多达 10.000 个航班。有了这些数字，MapBox 有时会导致浏览器崩溃或需要很长时间来渲染，特别是因为我们使用 arc.js 来计算 Great Circle 路线。

然后我们切换到将飞行路线创建为单个 geojson 文件，然后由 MapBox 直接加载。现在地图加载速度非常快（1-2 秒而不是半分钟），但路线是笔直的，不会越过日期变更线，这是一个真正的问题，因为这不是飞机飞行的方式。

首先，我在 MapBox 中寻找某种类型的设置，可以让我将线条渲染为大圆圈，但我找不到任何东西。接下来，我寻找类似于 arc.js 的 PHP 库来将路线输出到 geojson 文件，但是虽然有很多库可以根据大圆计算距离，但我没有找到任何实际生成路线的库。目前，我正在查看数据库级别。我们已经在使用PostGIS，所以我认为可能有一种方法可以使用它来计算路线。到目前为止，我在这里有这个，从各种来源拼凑而成，但它仍然抛出错误； ST_MakeLine 不存在...：

$ curl -D /dev/stdout -d mydata$'\n' 127.0.0.1:3005
HTTP/1.1 200 OK
Date: Wed Mar 31 2021 15:14:39 GMT-0400 (Eastern Daylight Time)
Connection: close
Content-Type: text/plain
Access-Control-Allow-Origin: *

POST / HTTP/1.1
Host: 127.0.0.1:3005
User-Agent: curl/7.74.0
Accept: */*
Content-Length: 7
Content-Type: application/x-www-form-urlencoded

mydata

我有点希望有一种隐藏的方式可以直接在 MapBox 中显示曲线，也许是通过源/图层？如果做不到这一点，我很乐意为 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 扩展。话虽如此，为了获得您正在寻找的曲线，您可以在 ST_Segmentize 之前将您的点转换为 spherical 投影参考系统。之后，为了获得单个 GeoJSON 字符串，您可能需要 ST_Collect 行，然后使用 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;

演示：db<>fiddle