PostGIS轨迹分析——横跨180°经线

问题描述

在处理AIS数据中，经常会遇到轨迹线横穿180°经线的情况，这种数据绘制到地图上显示的非常乱，如下图所示：

数据模拟

在geojson.io上模拟一条轨迹线，可以看到轨迹显示的非常好，红框里面的经纬度超过了180°，实际情况不会产生这种数据

使用真实的经纬度仍然还会出现错乱的情况，如下图所示：

解决思路

遍历轨迹上的每一个点，计算当前点和下一个点经度的差值(diff_x)，当diff_x的值大于340°时代表轨迹线是从左到右穿越，小于-340°时则代表轨迹线是从右到左穿越，这样我们就可以编写一个函数来解决此问题：

CREATE OR REPLACE FUNCTION get_lines_from_points(array_geom GEOMETRY[], array_diff_x FLOAT[]) RETURNS SETOF GEOMETRY AS
$$
DECLARE
    array_line  GEOMETRY[];
    line        GEOMETRY;
    array_point GEOMETRY[];
    diff_x      FLOAT;
    i           INT := 1;
BEGIN
    array_point := ARRAY []::GEOMETRY[];
    array_line := ARRAY []::GEOMETRY[];

    FOR diff_x IN SELECT unnest(array_diff_x)
        LOOP
            IF diff_x >= 340 THEN
                array_point := array_point || array_geom[i];
--                raise notice '%',i;
--                raise notice '%',st_astext(array_geom[i]);
                array_line := array_line || ST_MakeLine(array_point);
                array_point := ARRAY []::GEOMETRY[];
            ELSEIF diff_x <= -340 THEN
                array_point := array_point || array_geom[i];
                array_line := array_line || ST_MakeLine(array_point);
                array_point := ARRAY []::GEOMETRY[];
            ELSE
--                raise notice '%',i;
--                raise notice '%',st_astext(array_geom[i]);
                array_point := array_point || array_geom[i];
            END IF;
            i := i + 1;
        END LOOP;
    array_line := array_line || ST_MakeLine(array_point);

    FOR line IN SELECT unnest(array_line)
        LOOP
            RETURN NEXT line;
        END LOOP;
END;
$$
    LANGUAGE plpgsql;

调用函数的sql语句:

with line
         as (select st_geomfromgeojson('{"coordinates":[[-133.56681607072383,38.048432852702376],[-139.5625400372008,26.7466519608887],[-145.80029322035367,25.652134128092058],[-152.40600243048414,24.25006792947383],[-157.12509403666579,24.83850354437604],[-166.03562138914404,25.730242168027047],[-174.60130137567302,25.931369819197315],[-179.039207598898,25.851521506858006],[173.32935654204437,26.040298758510986],[171.01528349485716,29.02363444666598],[170.28599246631518,32.69415629123982],[172.7003594237434,35.35477927943893],[176.05091333553526,37.73333709776331],[-179.3654560402213,39.46966129070009],[-174.16271433310752,39.97444061559179],[-167.33463416628095,39.79967537154252],[-159.09064182365373,39.21533548155793],[-155.6068721801528,38.105767412202454],[-154.62754038722716,36.10105861232415],[-155.91220041908267,33.7318760829549],[-159.52676846180285,31.289767673068425],[-165.0776674179618,30.430887160278573],[-169.0175317306846,30.095326488395187],[-174.15399221361776,29.77057541216257],[-177.41433600934587,30.729135131803034],[178.7977321368836,32.191313894174456],[177.86505030800765,33.66780018177491]],"type":"LineString"}') as geom),
     origin_points as (select (st_dumppoints(line.geom)).path[1] as id, (st_dumppoints(line.geom)).geom as geom
                       from line),
     diff_points as (select origin_points.id                                               as id,
                            geom,
                            st_x(geom) - st_x(lead(geom) over (order by origin_points.id)) as diff_x
                     from origin_points),
     array_points as (select array_agg(geom) as array_geom, array_agg(diff_x) as array_diff_x
                      from diff_points)
select jsonb_build_object(
               'type', 'Feature',
               'geometry', st_asgeojson(get_lines_from_points(array_geom, array_diff_x))::json,
               'properties', jsonb_build_object()
       )
from array_points
;

结果返回了四个线段，知识点：

• st_dumppoints将模拟的线段打散为点
• lead函数获取下一行记录
• array_agg(geom) 和 array_agg(diff_x) 构建了两个用来存储经纬度和经度差的数组，用于传递给get_lines_from_points函数

在geojson.io中的展示效果，可以看到虽然说轨迹线没有之前那么乱了，但是在180°经线断开了

优化函数

当轨迹线穿越180°经线时，要在180°经线上补充一个点，同时下一条轨迹线也要在-180°经线上补充一个点，优化后的函数如下所示：

CREATE OR REPLACE FUNCTION get_lines_from_points(array_geom GEOMETRY[], array_diff_x FLOAT[]) RETURNS SETOF GEOMETRY AS
$$
DECLARE
    array_line  GEOMETRY[];
    line        GEOMETRY;
    array_point GEOMETRY[];
    temp_point  GEOMETRY;
    diff_x      FLOAT;
    middle_y    FLOAT;
    i           INT := 1;
BEGIN
    array_point := ARRAY []::GEOMETRY[];
    array_line := ARRAY []::GEOMETRY[];

    FOR diff_x IN SELECT unnest(array_diff_x)
        LOOP
            IF diff_x >= 340 THEN
                --                 raise notice '%',i;
--                 raise notice '%',st_astext(array_geom[i]);
                array_point := array_point || array_geom[i];
                -- 计算纬度中间值
                middle_y := (st_y(array_geom[i]) + st_y(array_geom[i + 1])) / 2;
                array_point := array_point || ST_SetSRID(ST_MakePoint(180, middle_y), 4326);
                -- 下一条轨迹线的第一个点要补充一个点
                temp_point := ST_SetSRID(ST_MakePoint(-180, middle_y), 4326);

                array_line := array_line || ST_MakeLine(array_point);
                array_point := ARRAY []::GEOMETRY[];
            ELSEIF diff_x <= -340 THEN
                array_point := array_point || array_geom[i];
                middle_y := (st_y(array_geom[i]) + st_y(array_geom[i + 1])) / 2;
                array_point := array_point || ST_SetSRID(ST_MakePoint(-180, middle_y), 4326);
                temp_point := ST_SetSRID(ST_MakePoint(180, middle_y), 4326);

                array_line := array_line || ST_MakeLine(array_point);
                array_point := ARRAY []::GEOMETRY[];
            ELSE
                IF NOT ST_IsEmpty(temp_point) THEN
                    -- 补充的点
                    array_point := array_point || temp_point;
                    temp_point := NULL;
                END IF;
                array_point := array_point || array_geom[i];
            END IF;
            i := i + 1;
        END LOOP;
    array_line := array_line || ST_MakeLine(array_point);

    FOR line IN SELECT unnest(array_line)
        LOOP
            RETURN NEXT line;
        END LOOP;
END;
$$
    LANGUAGE plpgsql;

最终结果

最终轨迹线被分割了若干份，效果如图所示：