# Seeking Algorithm to detect circling and beginning and end of circle?

I've got lots of flight data from glider pilots in the form of gps fixes in a fixed interval. I would like to analyze the flight path and detect the start and the end of the 'circling' the glider pilot will do when he finds thermals.

Ideally, an algorithm would give me a start and an end point on the line, defining one "circle". These points could be equal to one of the gps fixes and do not need to be interpolated.

I simply could walk along the flight path, check the turn-rate and have some criteria to decide if the glider is circling or not.

As I am using PostgreSQL with PostGIS extension, I was curious if there is a better approach to this issue. I already have an procedure to calculate the angle of two line segments:

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

It should be possible to loop over all line segments and check, when the sum of the angles is greater than 360 or less than -360 degrees. Then I could use st_centroid to detect the center of the circle, if needed.

Is there a better approach?

As requested, I uploaded an example flight. • Looking around kicked up Hough Circle Transform. There's a similar (with a polygon though) postgis user discussion here: lists.osgeo.org/pipermail/postgis-users/2015-February/… – Barrett Aug 12 '16 at 15:18
• Thank you both. I'll have a look at the Hough Transform. The discussion at osgeo.org assumes that I already know where the circle starts and ends, if I understood it correctly? – pgross Aug 12 '16 at 15:28
• Have you seen this: gis.stackexchange.com/questions/37058 – Devdatta Tengshe Aug 12 '16 at 15:59
• @DevdattaTengshe Yes, but thanks anyway. That would be an approach where I would have to calculate the splines and the curvature externally, right? By external, I mean not as a procedure or query directly on the database. As the flights do not change, once they are in the database, this would be an option. – pgross Aug 12 '16 at 16:06
• Can you post some sample data as a .sql file? – dbaston Aug 12 '16 at 16:49

I couldn't stop thinking about this... I was able to come up with a Stored Procedure to do the loop counting. The example path contains 109 loops!

Here are the flight points shown with the loop centroids in red: Basically, it runs through the points in the order they were captured and builds a line as it iterates through the points. When the line we are building creates a loop (using ST_BuildArea) then we count a loop and start building a line again from that point.

This function returns a recordset of each loop which contains the loop number, its geometry, its start/end point and its centroid (I also cleaned it up a bit and made better variable names):

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

This is a simple function to return only the loop count:

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

• This looks very promising. Thank a lot. I will need to enhance it, as I am not interested in the number of circles but the start/end points. But that should be easily possible to return I guess. – pgross Sep 2 '16 at 5:18
• That sounds pretty clever. How does it handle the situation where one loop intersects another loop? Or are you skipping the initial points once you locate a loop? – Peter Horsbøll Møller Sep 2 '16 at 9:11
• @PeterHorsbøllMøller It analyzes when the line makes a loop (ST_BuildArea will only return true when the line creates a closed area) rather than looking for intersections. – kttii Sep 2 '16 at 12:32
• @pgross oops true! I got a bit sidetracked and forgot about the start/end points but yes, that's an easy enough determination now that loops are distinguished. – kttii Sep 2 '16 at 12:34
• @pgross It seems to me that you would probably get more reasonable locations of the thermals by locating the ST_Centroid of each loop rather than locating the start/end of each loop. What do you think? Of course, the function could provide all three statistics. – kttii Sep 2 '16 at 12:37

I noticed that the gpx file has time stamp which could be exploited. Perhaps the below approach could work.

``````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'
``````
• I found it difficult to use ST_Intersects because of overlapping circles which led me to use ST_BuildArea. – kttii Sep 1 '16 at 16:21
• I gave you the bounty since your answer is generally on the same track. – kttii Sep 6 '16 at 15:51