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.383,38.0376],[-139.08,26.87],[-145.367,25.2058],[-152.414,24.383],[-157.579,24.604],[-166.004,25.7047],[-174.302,25.7315],[-179.0,25.8006],[173.437,26.0986],[171.016,29.098],[170.518,32.982],[172.34,35.893],[176.026,37.331],[-179.13,39.009],[-174.752,39.179],[-167.095,39.252],[-159.073,39.793],[-155.28,38.2454],[-154.716,36.415],[-155.267,33.49],[-159.285,31.8425],[-165.08,30.8573],[-169.06,30.0187],[-174.776,29.257],[-177.587,30.3034],[178.36,32.4456],[177.765,33.491]],"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;

最终结果

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