Suchst du nach einem Algorithmus zum Erkennen von Kreisen sowie des Beginns und Endes eines Kreises?

24

Ich habe viele Flugdaten von Segelflugpiloten in Form von GPS-Fixes in einem festgelegten Intervall. Ich möchte die Flugbahn analysieren und den Beginn und das Ende des "Kreisens" erkennen, das der Segelflugpilot machen wird, wenn er Thermik findet.

Im Idealfall würde mir ein Algorithmus einen Start- und einen Endpunkt auf der Linie geben und einen "Kreis" definieren. Diese Punkte können einem der GPS-Fixes entsprechen und müssen nicht interpoliert werden.

Ich könnte einfach den Flugweg entlang gehen, die Wendegeschwindigkeit überprüfen und anhand einiger Kriterien entscheiden, ob der Schirm kreisen soll oder nicht.

Da ich PostgreSQL mit PostGIS-Erweiterung verwende, war ich neugierig, ob es einen besseren Ansatz für dieses Problem gibt. Ich habe bereits ein Verfahren zur Berechnung des Winkels zweier Liniensegmente:

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;

Es sollte möglich sein, alle Liniensegmente zu durchlaufen und zu überprüfen, ob die Summe der Winkel größer als 360 oder kleiner als -360 Grad ist. Dann könnte ich st_centroid verwenden, um bei Bedarf den Mittelpunkt des Kreises zu ermitteln.

Gibt es einen besseren Ansatz?


Wie gewünscht habe ich einen Beispielflug hochgeladen .

Eine Beispielflugbahn

pgross
quelle
1
Hough Circle Transform blickte sich um. Es gibt eine ähnliche (mit einem Polygon versehene) Postgis-Benutzerdiskussion hier: lists.osgeo.org/pipermail/postgis-users/2015-February/…
Barrett
Danke euch beiden. Ich werde einen Blick auf die Hough-Transformation werfen. Bei der Diskussion auf osgeo.org wird davon ausgegangen, dass ich bereits weiß, wo der Kreis beginnt und endet, wenn ich ihn richtig verstanden habe.
Freitag,
Haben Sie dies gesehen: gis.stackexchange.com/questions/37058
Devdatta Tengshe
@ DevdattaTengshe Ja, aber trotzdem danke. Das wäre ein Ansatz, bei dem ich die Splines und die Krümmung von außen berechnen müsste, oder? Mit extern meine ich nicht als Prozedur oder Abfrage direkt in der Datenbank. Da sich die Flüge nicht ändern, sobald sie in der Datenbank vorhanden sind, ist dies eine Option.
Freitag,
Können Sie einige Beispieldaten als .sql-Datei veröffentlichen?
Dbaston

Antworten:

14

Ich konnte nicht aufhören darüber nachzudenken ... Ich konnte mir eine gespeicherte Prozedur ausdenken, um die Schleife zu zählen. Der Beispielpfad enthält 109 Schleifen!

Hier sind die Flugpunkte mit den Schleifenschwerpunkten in Rot dargestellt: Bildbeschreibung hier eingeben

Grundsätzlich werden die Punkte in der Reihenfolge durchlaufen, in der sie erfasst wurden, und es wird eine Linie erstellt, während die Punkte durchlaufen werden. Wenn die gerade erstellte Linie eine Schleife erstellt (mit ST_BuildArea), zählen wir eine Schleife und beginnen von diesem Punkt aus erneut mit dem Erstellen einer Linie.

Diese Funktion gibt eine Datensatzmenge jeder Schleife zurück, die die Schleifennummer, ihre Geometrie, ihren Start- / Endpunkt und ihren Schwerpunkt enthält (ich habe sie auch ein wenig aufgeräumt und bessere Variablennamen erstellt):

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;

Dies ist eine einfache Funktion, um nur die Schleifenzahl zurückzugeben:

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;

kttii
quelle
Das sieht sehr vielversprechend aus. Vielen Dank. Ich muss es verbessern, da mich nicht die Anzahl der Kreise interessiert, sondern die Start- / Endpunkte. Aber das sollte wohl problemlos möglich sein, denke ich.
pgross
Das klingt ziemlich schlau. Wie geht es mit der Situation um, in der eine Schleife eine andere Schleife schneidet? Oder überspringen Sie die Anfangspunkte, sobald Sie eine Schleife gefunden haben?
Peter Horsbøll Møller
@ PeterHorsbøllMøller Analysiert, wann die Linie eine Schleife bildet (ST_BuildArea gibt nur dann true zurück, wenn die Linie einen geschlossenen Bereich erstellt), anstatt nach Schnittpunkten zu suchen.
kttii
@pgross oops true! Ich habe mich ein wenig ablenken lassen und die Start- / Endpunkte vergessen, aber ja, das ist jetzt, da die Schleifen unterschieden werden, leicht genug festzustellen.
kttii
@pgross Es scheint mir, dass Sie wahrscheinlich vernünftigere Positionen der Thermik erhalten würden, wenn Sie den ST_Centroid jeder Schleife lokalisieren, anstatt den Anfang / das Ende jeder Schleife zu lokalisieren. Was denkst du? Natürlich könnte die Funktion alle drei Statistiken liefern.
kttii
3

Mir ist aufgefallen, dass die GPX-Datei einen Zeitstempel hat, der ausgenutzt werden könnte. Vielleicht könnte der folgende Ansatz funktionieren.

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' 
addcolor
quelle
Ich fand es schwierig, ST_Intersects zu verwenden, weil sich Kreise überlappten, was mich dazu veranlasste, ST_BuildArea zu verwenden.
kttii
Ich habe dir das Kopfgeld gegeben, da deine Antwort im Allgemeinen auf dem gleichen Weg ist.
kttii