24

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.

A sample flight path

  • 1
    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
14

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: enter image description here

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
3
+200

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.