Buscando algoritmo para detectar circulando e começo e fim de círculo?

24

Eu tenho muitos dados de vôo de pilotos de planador na forma de correções de GPS em um intervalo fixo. Gostaria de analisar a trajetória de vôo e detectar o início e o final do 'círculo' que o piloto de planador fará quando encontrar as térmicas.

Idealmente, um algoritmo me daria um ponto inicial e final na linha, definindo um "círculo". Esses pontos podem ser iguais a uma das correções de GPS e não precisam ser interpolados.

Eu simplesmente poderia caminhar ao longo da trajetória de vôo, verificar a taxa de curva e ter alguns critérios para decidir se o planador está circulando ou não.

Como estou usando o PostgreSQL com extensão PostGIS, fiquei curioso para saber se existe uma abordagem melhor para esse problema. Eu já tenho um procedimento para calcular o ângulo de dois segmentos de linha:

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;

Deve ser possível fazer um loop em todos os segmentos de linha e verificar se a soma dos ângulos é maior que 360 ​​ou menor que -360 graus. Então eu poderia usar o st_centroid para detectar o centro do círculo, se necessário.

Existe uma abordagem melhor?


Conforme solicitado, enviei um exemplo de voo .

Um exemplo de trajetória de vôo

pgross
fonte
11
Olhando em volta, chutou o Hough Circle Transform. Há uma semelhante (com um polígono embora) PostGIS discussão utilizador aqui: lists.osgeo.org/pipermail/postgis-users/2015-February/...
Barrett
Obrigado a ambos. Vou dar uma olhada na Hough Transform. A discussão no osgeo.org supõe que eu já saiba onde o círculo começa e termina, se o entendi corretamente?
Pgross
Você já viu isso: gis.stackexchange.com/questions/37058
Devdatta Tengshe
@DevdattaTengshe Sim, mas obrigado de qualquer maneira. Essa seria uma abordagem em que eu teria que calcular os splines e a curvatura externamente, certo? Por externo, não quero dizer como um procedimento ou consulta diretamente no banco de dados. Como os voos não mudam, uma vez que estejam no banco de dados, isso seria uma opção.
Pbruta
Você pode postar alguns dados de amostra como um arquivo .sql?
Dbaston

Respostas:

14

Eu não conseguia parar de pensar nisso ... Consegui criar um Procedimento Armazenado para fazer a contagem do loop. O caminho do exemplo contém 109 loops!

Aqui estão os pontos de vôo mostrados com os centróides do loop em vermelho: insira a descrição da imagem aqui

Basicamente, ele percorre os pontos na ordem em que foram capturados e constrói uma linha conforme itera pelos pontos. Quando a linha que estamos construindo cria um loop (usando ST_BuildArea), contamos um loop e começamos a construir uma linha novamente a partir desse ponto.

Essa função retorna um conjunto de registros de cada loop que contém o número do loop, sua geometria, seu ponto inicial / final e seu centróide (eu também o limpei um pouco e criei melhores nomes de variáveis):

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;

Esta é uma função simples para retornar apenas a contagem de loop:

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
fonte
Isso parece muito promissor. Muito obrigado. Vou precisar aprimorá-lo, pois não estou interessado no número de círculos, mas nos pontos de partida / chegada. Mas isso deve ser facilmente possível retornar, eu acho.
pgross 02/09
Isso parece muito inteligente. Como ele lida com a situação em que um loop cruza outro loop? Ou você está pulando os pontos iniciais depois de localizar um loop?
Peter Horsbøll Møller 02/09/16
@ PeterHorsbøllMøller Analisa quando a linha faz um loop (ST_BuildArea retornará true somente quando a linha criar uma área fechada) em vez de procurar interseções.
kttii
@pgross oops verdade! Fiquei um pouco desviado e esqueci os pontos de início / fim, mas sim, essa é uma determinação fácil o suficiente agora que os loops são diferenciados.
kttii
@pgross Parece-me que você provavelmente obteria locais mais razoáveis ​​das térmicas localizando o ST_Centroid de cada loop em vez de localizar o início / fim de cada loop. O que você acha? Obviamente, a função poderia fornecer todas as três estatísticas.
kttii
3

Notei que o arquivo gpx tem carimbo de data e hora que pode ser explorado. Talvez a abordagem abaixo possa funcionar.

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
fonte
Achei difícil usar o ST_Intersects por causa da sobreposição de círculos, o que me levou a usar o ST_BuildArea.
kttii 1/09/16
Eu lhe dei a recompensa, já que sua resposta geralmente está na mesma trilha.
kttii