理想情况下,一种算法会给我直线的起点和终点,定义一个“圆”。这些点可以等于GPS修正值之一,不需要插补。是否盘旋。
当我使用带有PostGIS扩展功能的PostgreSQL时,我很好奇是否有更好的方法来解决这个问题。
我已经有一个程序可以计算两个线段:
CREATE OR REPLACE FUNCTION angle_between(
_p1 GEOMETRY(PointZ,4326),
_p2 GEOMETRY(PointZ,4326),
_p3 GEOMETRY(PointZ,4326)
) RETURNS DECIMAL AS $$
DECLARE
az1 FLOAT;
az3 FLOAT;
BEGIN
az1 = st_azimuth(_p2,_p1);
az3 = st_azimuth(_p2,_p3);
IF az3 > az1 THEN
RETURN (
degrees(az3 - az1)::decimal - 180
);
ELSE
RETURN (
degrees(az3 - az1)::decimal + 180
);
END IF;
END;
$$ LANGUAGE plpgsql;
当角度之和大于360度或小于-360度时,应该可以遍历所有线段并进行检查。然后,如果需要,我可以使用st_centroid来检测圆心。
有更好的方法吗?
根据要求,我上传了一个示例航班。
#1 楼
我一直不停地思考这个问题……我能够提出一个存储过程来进行循环计数。示例路径包含109个循环!以下是带有红色质点的循环质心显示的飞行点:
基本上,它穿过依次捕获它们,并在遍历各个点时建立一条线。当我们正在构建的行创建一个循环(使用ST_BuildArea)时,我们计算一个循环并从该点开始再次构建一条线。它的几何形状,它的起点/终点和质心(我也对其进行了一些清理,并命名为更好的变量):
DROP FUNCTION test.find_loop_count(flightid int);
create function test.find_Loop_count(
IN flightid int,
OUT loopnumber int,
OUT loopgeometry geometry,
OUT loopstartend geometry,
OUT loopcentroid geometry
)
RETURNS SETOF record AS
$BODY$
-- s schema 'test' must exist
-- a table 'points' of flight points must exist
-- we are going to iterate through the point path, building a line as we go
-- If the line creates a loop then we count a loop and start over building a new line
-- add the intersection point to the returning recordset
-- add the centroid of the loop to the resulting recordset
-- pass in the flight ID of the flight that you wish to count its loops for example:
-- SELECT * FROM find_loop_count(37);
DECLARE
rPoint RECORD;
gSegment geometry = NULL;
gLastPoint geometry = NULL;
gLoopPolygon geometry = NULL;
gIntersectionPoint geometry = NULL;
gLoopCentroid geometry = NULL;
iLoops integer := 0;
BEGIN
-- for each line segment in Point Path
FOR rPoint IN
WITH
pts as (
SELECT location as geom,datetime,row_number() OVER () as rnum
FROM test.points
WHERE flight_id=flightid
ORDER BY 2)
SELECT ST_AsText(ST_MakeLine(ARRAY[a.geom, b.geom])) AS geom, a.rnum, b.rnum
FROM pts as a, pts as b
WHERE a.rnum = b.rnum-1 AND b.rnum > 1
LOOP
-- if this is the start of a new line then start the segment otherwise add the point to the segment
if gSegment is null then
gSegment=rPoint.geom;
elseif rPoint.geom::geometry=gLastPoint::geometry then
-- do not add this point to the segment because it is at the same location as the last point
else
-- add this point to the line
gSegment=ST_Makeline(gSegment,rPoint.geom);
end if;
-- ST_BuildArea will return true if the line segment is noded and closed
-- we must also flatten the line to 2D
-- lets also make sure that there are more than three points in our line to define a loop
gLoopPolygon=ST_BuildArea(ST_Node(ST_Force2D(gSegment)));
if gLoopPolygon is not NULL and ST_Numpoints(gSegment) > 3 then
-- we found a loop
iLoops:=iLoops+1;
-- get the intersection point (start/end)
gIntersectionPoint=ST_Intersection(gSegment::geometry,rPoint.geom::geometry);
-- get the centroid of the loop
gLoopCentroid=ST_Centroid(gLoopPolygon);
-- start building a new line
gSegment=null;
LOOPNUMBER := iLoops;
LOOPGEOMETRY := gLoopPolygon;
LOOPSTARTEND := gIntersectionPoint;
LOOPCENTROID := gLoopCentroid;
RETURN NEXT;
end if;
-- keep track of last segment
gLastPoint=rPoint.geom;
END LOOP;
RAISE NOTICE 'Total loop count is %.', iLoops;
END;
$BODY$
LANGUAGE plpgsql STABLE
COST 100
ROWS 1000;
>这是一个仅返回循环计数的简单函数:DROP FUNCTION test.find_loop_count(flightid int);
create function test.find_Loop_count(flightid int) RETURNS integer AS $$
-- s schema 'test' must exist
-- a table 'points' of flight points must exist
-- we are going to iterate through the line path, building the line as we go
-- If the line creates a loop then we count a loop and start over building a new line
-- pass in the flight ID of the flight that you wish to count its loops for example:
-- SELECT find_loop_count(37);
DECLARE
segment RECORD;
s geometry = NULL;
lastS geometry = NULL;
b geometry = NULL;
loops integer := 1;
BEGIN
-- for each line segment is Point Path
FOR segment IN
WITH
pts as (
SELECT location as geom,datetime,row_number() OVER () as rnum
FROM test.points
WHERE flight_id=flightid
ORDER BY 2)
SELECT ST_AsText(ST_MakeLine(ARRAY[a.geom, b.geom])) AS geom, a.rnum, b.rnum
FROM pts as a, pts as b
WHERE a.rnum = b.rnum-1 AND b.rnum > 1
LOOP
-- if this is the start of a new line then make s be the segment otherwise add the segment to s
if s is null then
s=segment.geom;
elseif segment.geom::geometry=lastS::geometry then
else
s=ST_Makeline(s,segment.geom);
end if;
-- ST_BuildArea will return true if the line segment is noded and closed
-- we must also flatten the line to 2D
b=ST_BuildArea(st_node(ST_Force2D(s)));
if b is not NULL and st_numpoints(s) > 3 then
RAISE NOTICE 's: %', s;
RAISE NOTICE 'vvvvv %',st_numpoints(s);
RAISE NOTICE 'I found a loop! Loop count is now %', loops;
RAISE NOTICE '^^^^^';
s=null;
loops:=loops +1;
end if;
lastS=segment.geom;
END LOOP;
RAISE NOTICE 'Total loop count is %.', loops-1;
RETURN loops-1;
END;
$$ LANGUAGE plpgsql;
评论
这看起来很有希望。非常感谢。我将需要增强它,因为我对圆圈数不感兴趣,但对起点/终点不感兴趣。但是我想那应该很容易返回。
–pgross
16年9月2日,下午5:18
听起来很聪明。如何处理一个循环与另一个循环相交的情况?还是在找到循环后就跳过初始点?
– PeterHorsbøllMøller
16-09-2在9:11
@PeterHorsbøllMøller它分析线何时循环(ST_BuildArea仅在线创建闭合区域时返回true),而不是寻找交点。
–kttii
16年9月2日在12:32
@pgross糟糕!我有些困惑,忘了起点/终点,但是是的,既然可以区分出循环,那么这很容易确定。
–kttii
16-09-2在12:34
@pgross在我看来,通过定位每个循环的ST_Centroid而不是定位每个循环的开始/结束,您可能会获得更合理的散热位置。你怎么看?当然,该功能可以提供所有三个统计信息。
–kttii
16-09-2在12:37
#2 楼
我注意到gpx文件具有可以利用的时间戳。以下方法也许可以使用。
Make a linesegement with Vi,Vi+1
Make it Polyline
Proceed to Vi+2,Vi+3 check intersection with Polyline
if it intersects
find the point of intersection-Designate this as start/end point of the loop
Make this intersection point as Vi and Vi+1 would next gpx point per time sequence
if the linesegement does not intersect with polyyline then increment 'i'
评论
我发现很难使用ST_Intersects,因为圆圈重叠导致我不得不使用ST_BuildArea。
–kttii
2016年9月1日下午16:21
我给了您赏金,因为您的答案通常在同一轨道上。
–kttii
2016年9月6日15:51
评论
环顾四周,引发了霍夫环变换。这里有一个类似的(尽管带有多边形)postgis用户讨论:lists.osgeo.org/pipermail/postgis-users/2015-February/…谢谢你们俩。我来看看霍夫变换。 osgeo.org上的讨论假设我已经知道了圆的起点和终点,如果我正确理解的话?
您是否看到过:gis.stackexchange.com/questions/37058
@DevdattaTengshe是的,但是还是谢谢。那是一种我必须从外部计算样条曲线和曲率的方法,对吗?外部的,我的意思不是作为过程或直接在数据库上查询。由于航班不会更改,因此一旦将其保存在数据库中,将是一种选择。
您可以将一些示例数据发布为.sql文件吗?