Cálculo de um centróide de polígono esférico

11

Eu gostaria de uma maneira geral de calcular centróides para polígonos em uma esfera.

Até agora, a melhor referência online parece ser:

Ferramentas para gráficos e formas, por Jeff Jenness.

O método descrito sugere a decomposição do polígono em múltiplos triângulos esféricos e o cálculo da média dos centróides esféricos do triângulo, ponderados pela área do triângulo esférico.

Eu sei que existem várias maneiras de definir um centróide de polígono esférico, mas estou procurando algo análogo às seguintes definições de pontos e polilinhas:

  • Pontos : média aritmética dos vetores cartesianos que representam os pontos.
  • Polilinhas : média ponderada de vetores cartesianos representando pontos médios de cada segmento de linha, ponderados pelo comprimento (esférico) de cada segmento.

Parece uma continuação razoável ter centróides poligonais definidos como a média ponderada da decomposição triangular, ponderada por área.

Minha pergunta é se o método na referência acima funcionará independentemente da decomposição do triângulo usada. Em particular, ele menciona a decomposição em triângulos em relação a um ponto arbitrário, mesmo externo ao polígono, de modo que alguns triângulos terão áreas negativas que contribuem com um peso negativo.

Relacionado: Como encontrar o centro da geometria de um objeto?

Jason Davies
fonte

Respostas:

9

Ele não funcionará de maneira consistente, mesmo quando você executar todas as triangulações em relação a um único ponto fixo . O problema é que os cálculos esféricos e euclidianos estão sendo misturados sem considerar o que eles significam.

Uma maneira de tornar isso óbvio é considerar um triângulo extremo, como quase a metade de um hemisfério. Por exemplo, começando em (lon, lat) = (-179, 0), corra ao longo do equador até (0, 0), depois suba até o polo norte em (0, 90) e volte ao início em (- 179, 0). Este é um triângulo 90-179-90 que compreende a maior parte da metade norte do hemisfério ocidental. O problema é que seus pontos finais (mostrados como pontos brancos na figura) estão praticamente em um plano: um está no poste e os outros dois estão quase em lados opostos. Assim, sua média, projetada de volta à esfera (o ponto vermelho), está quase no mastro - mas isso é o mais longe possível de qualquer centro razoável:

Triângulo esférico grande

Como outro exemplo, vamos triangular um polígono representando o hemisfério superior em relação ao seu centro, o Polo Norte. Sempre dividiremos o hemisfério ocidental em duas metades iguais, cada uma com um triângulo 90-90-90 (evitando assim qualquer problema com enormes triângulos que se estendem pelo hemisfério). O hemisfério oriental, no entanto, será dividido em n semilunas iguais. Os vértices do lune k ( k = 1, 2, ..., n ) têm coordenadas (lon, lat)

((k-1) * 180/n, 0),  (k * 180/n, 0),  (k * 180/n, 90).

Lunes para k = 8

Esta figura mostra a configuração para k = 8. Os pontos vermelhos são os "centros" individuais do triângulo calculados de acordo com o documento "Ferramentas para gráficos e formas", páginas 65-67.

Fazendo os cálculos, acho que com k = 2, o centro ponderado por área está realmente no Polo Norte (como seria indicado por considerações de simetria), mas à medida que n aumenta, o resultado muda rapidamente para o hemisfério ocidental e, no limite, aproxima-se de uma latitude de 89.556 graus ao longo da longitude -90 graus. Isso fica a aproximadamente 50 quilômetros ao sul do próprio Polo Norte.

É certo que um erro de +/- 50 quilômetros para um polígono de 20.000 quilômetros é pequeno; a quantidade total de variação arbitrária devido a diferentes triangulações neste caso é de apenas 0,5%. Evidentemente, os erros relativos podem ser arbitrariamente grandes, incluindo triângulos negativos (basta adicionar e subtrair alguns triângulos realmente grandes em relação a um triângulo pequeno). Independentemente disso, qualquer um que se esforce para fazer cálculos esféricos evidentemente está tentando evitar erros de projeção, portanto está procurando alta precisão. Este método de triangulação não pode ser recomendado.

whuber
fonte
Você demonstrou que erros podem acumular-se para n grande, mas não está claro que a abordagem seja necessariamente falha. Qual valor de n você usou para atingir o valor limite?
perfil completo de Jason Davies
Além disso, muito obrigado por fazer os cálculos e analisar isso em profundidade. Eu só quero um pouco mais de esclarecimento antes que eu possa resolver o problema. :)
Jason Davies
Jason, adicionei um exemplo preliminar para lhe dar alguma intuição. O próprio limite é rapidamente abordado; algumas dúzias de lunas oferecem vários dígitos significativos. Mas o novo exemplo deve repousar qualquer dúvida persistente de que essa triangulação ponderada faça algo razoável - exceto em pequenos triângulos, nos quais é muito melhor você fazer os cálculos nas coordenadas projetadas. A única razão para fazer cálculos esféricos é quando sua área de análise é verdadeiramente global, pois todas as projeções introduzem muita distorção.
whuber
1
Fantástico, obrigado. Portanto, se eu entendi corretamente, simplesmente calcular a média dos vetores cartesianos não resulta em um centróide razoável para um triângulo esférico (particularmente grandes como o seu primeiro exemplo). Vou investigar métodos melhores, por exemplo, encontrar a interseção de medianas do grande círculo.
Jason Davies
BTW, ainda estou esperançoso de que um centróide esférico de área semelhante ao anterior funcione. Imagine cada polígono recebendo um volume 3D adicionando um vértice na origem da esfera. Em seguida, suspenda a esfera por uma corda invisível conectada à sua origem e encontre um equilíbrio estável. O centróide é o ponto mais baixo (é a projeção do centro de massa na superfície esférica). Isso deve funcionar além de alguns casos ambíguos, por exemplo, uma faixa ao redor do equador, onde posso apenas escolher um ponto sensato. É um prazer discutir em uma nova pergunta se você acha que vale a pena.
21412 Jason Jason
3

É uma boa idéia enumerar propriedades que o centróide de um polígono deve ter. Aqui estão os meus critérios:

(a) É uma propriedade do interior do polígono (em vez dos vértices ou arestas). Assim, dividir uma aresta em duas inserindo um vértice adicional não deve alterar a posição do centróide. Observe que a definição de centróide de Jenness falha nesse critério, pois a posição do centróide depende de como um polígono é dividido em triângulos.

(b) Perturbar um pouco a forma do polígono deve mover um pouco o centróide. Aqui é necessário impor uma restrição à extensão geral do polígono (por exemplo, a um único hemisfério). Sem essa restrição, é fácil construir casos em que o centróide suba de repente para o lado oposto da terra com um leve movimento de um vértice. Essa condição exclui métodos que exigem que o centróide esteja dentro do polígono.

(c) Deve reduzir à definição planar de centróide para pequenos polígonos.

Aqui estão duas abordagens atendem a esses critérios:

(1) Calcule o centróide para o polígono elipsóide em três dimensões e projete-o de volta à superfície elipsóide (ao longo de um normal para o elipsóide). Grande vantagem: o centróide pode ser calculado dividindo o polígono em formas mais simples.

(2) O centróide é o ponto com uma distância geodésica mínima de RMS para todos os pontos no interior do polígono. Veja Buss e Fillmore, "Médias esféricas e aplicações para splines esféricas e interpolação", ACM Transactions on Graphics 20 , 95–126 (2001). Grande vantagem: o ponto resultante não depende de como a superfície é incorporada no R 3 .

Infelizmente, nenhuma dessas definições é simples de colocar em prática. No entanto , o primeiro método pode ser realizado simplesmente para uma esfera. A melhor área "elementar" a ser usada é o quadrilátero delimitado por uma aresta do polígono, dois meridianos através dos pontos finais da aresta e o equador. O resultado para todo o polígono implica somar as contribuições sobre as arestas. (É necessário executar etapas adicionais se o polígono envolver um poste.)

Suponha que os pontos finais da aresta sejam (φ 1 , λ 1 ) e (φ 2 , λ 2 ). Deixe os azimutes da aresta e os pontos finais em α 1 e α 2 . Assumindo que o raio da esfera é 1, a área do quadrilátero é

  A = α 2 - α 1
      = 2 tan −1 [tan ½ (λ 2 - λ 1 ) sen ½ (φ 2 + φ 1 ) / cos ½ (φ 2 + φ 1 )]

(Esta fórmula para a área, devido a Bessel, é substancialmente melhor comportada numericamente do que a fórmula de L'Huilier comumente usada para a área de um triângulo.)

Os componentes do centróide para este quadrilátero são dados por

  2 Umx ⟩ = φ 2 sin (λ 2 - λ 0 ) - φ 1 sin (λ 1 - λ 0 )
  2 Umay ⟩ = cos α 02 - σ 1 ) - (φ 2 cos (λ 2 - λ 0 ) - φ 1 cos (λ 1 - λ 0 ))
  2 Umaz ⟩ = (λ 2 - λ 1 ) - sin α 02 - σ1 )

onde σ 2 - σ 1 é o comprimento da aresta, e λ 0 e α 0 são a longitude e o azimute da aresta em que ela atravessa o equador, e os eixos x e y são orientados de modo que o cruzamento do equador esteja em x = 1, y = 0. ( z é o eixo através do polo, é claro).

cffk
fonte
Você pode explicar por que a posição dos centróides de Jenness dependerá de como um polígono é dividido em triângulos? Sei pelo exemplo da @ whuber que o cálculo do centróide de Jenness é errado para triângulos esféricos, mas e se um centróide baseado em medianas do triângulo esférico for usado? Isso ainda falhará?
Jason Davies
Jenness efetivamente substitui o polígono esférico por um conjunto de triângulos planares e calcula seu centróide. Claramente (?), O resultado dependerá do particionamento. Fazer o cálculo que descrevi usando os centróides de triângulos esféricos é bom. Você pode encontrar a fórmula para o centróide em JE Brock, O tensor de inércia para um triângulo esférico, J. Applied Mechanics 42, 239 (1975) dx.doi.org/10.1115/1.3423535
cffk
Dei outra olhada no jornal de Brock. Sua fórmula para o centro de massa de um triângulo esférico envolve uma soma sobre as bordas do triângulo. Portanto, pode ser trivialmente generalizado para aplicar a um polígono (sem a necessidade de dividi-lo em triângulos).
Cffk 28/12/12
Você também se importa em fornecer uma referência para o cálculo da área devido a Bessel? Não consigo encontrá-lo em lugar nenhum, e estou interessado em escrever uma rotina de área de polígono esférica rápida (e precisa). Obrigado!
Jason Davies
Encontrei e percebi que você a traduziu para o inglês, então, obrigado. :)
Jason Davies