O Côncavo Casco

Criação de um cluster de fronteira usando um K-vizinhos mais próximos abordagem

Um par de meses atrás, Eu escrevi aqui no Meio de um artigo no mapeamento do reino UNIDO acidentes de trânsito (hot spots). Estava preocupado em ilustrar o uso do algoritmo de clustering DBSCAN em dados geográficos. No artigo, utilizei a informação geográfica publicada pelo governo do Reino Unido sobre os acidentes de viação notificados. O meu objectivo era executar um processo de agrupamento baseado na densidade para encontrar as áreas onde os acidentes de trânsito são relatados mais frequentemente. O resultado final foi a criação de um conjunto de cercas que representam estes pontos quentes de acidentes.

ao coletar todos os pontos em um dado conjunto, você pode ter uma idéia de como o conjunto seria em um mapa, mas você vai precisar de uma parte importante de informação: a forma exterior do aglomerado. Neste caso, estamos a falar de um polígono fechado que pode ser representado num mapa como uma geo-cerca. Qualquer ponto dentro da geo-cerca pode ser assumido como pertencendo ao conjunto, o que faz desta forma uma informação interessante: você pode usá-la como uma função discriminadora. Todos os pontos recentemente amostrados que caem dentro do polígono pode ser assumido como pertencendo ao aglomerado correspondente. Como eu insinuei no artigo, você poderia usar tais polígonos para afirmar o seu risco de condução, usando-os para classificar a sua própria localização GPS amostrado.

a questão agora é como criar um polígono significativo a partir de uma nuvem de pontos que compõem um conjunto específico. Minha abordagem no primeiro artigo foi um pouco ingênua e refletiu uma solução que eu já tinha usado em código de produção. Esta solução implicava a colocação de um círculo centrado em cada um dos pontos do aglomerado e, em seguida, a fusão de todos os círculos juntos para formar um polígono em forma de nuvem. O resultado não é muito agradável, nem realista. Além disso, ao usar círculos como a forma base para a construção do polígono final, estes terão mais pontos do que uma forma mais simplificada, aumentando assim os custos de armazenamento e tornando os algoritmos de detecção de inclusão mais lentos de executar.

por outro lado, esta abordagem tem a vantagem de ser computacionalmente simples (pelo menos da perspectiva do desenvolvedor), como ele usa bem Torneadas do cascaded_union função para mesclar todos os círculos juntos. Outra vantagem é que a forma do polígono é implicitamente definida usando todos os pontos do aglomerado.

para abordagens mais sofisticadas, precisamos de alguma forma identificar os pontos de fronteira do aglomerado, os que parecem definir a forma da nuvem do ponto. Curiosamente, com algumas implementações DBSCAN você pode realmente recuperar os pontos de fronteira como um subproduto do processo de agrupamento. Infelizmente, esta informação não está (aparentemente) disponível na implementação do SciKit Learn, por isso temos que fazer.

a primeira abordagem que surgiu na mente foi calcular o casco convexo do conjunto de pontos. Esta é uma bem entendida algoritmo, mas sofre com o problema da não manipulação de formas côncavas, como este:

Desta forma não corretamente capturar a essência subjacente pontos. Se ele fosse usado como um discriminador, alguns pontos seriam incorretamente classificados como estando dentro do aglomerado quando eles não estão. Precisamos de outra abordagem.

a alternativa do casco côncavo

Felizmente, Existem alternativas a este estado de coisas: podemos calcular um casco côncavo. Aqui está o que o côncavo casco parece que quando aplicado ao mesmo conjunto de pontos, como na imagem anterior:

Ou, talvez, um:

Como você pode ver, e ao contrário do casco convexo, não existe uma definição única do que a concavidade do casco de um conjunto de pontos. Com o algoritmo que estou apresentando aqui, a escolha de quão côncavo você quer que seus cascos sejam feitos através de um único parâmetro: k — o número de vizinhos mais próximos considerados durante o cálculo do casco. Vamos ver como isto funciona.

o algoritmo

o algoritmo que estou a apresentar aqui foi descrito há mais de uma década por Adriano Moreira e Maribel Yasmina Santos da Universidade de Minho, Portugal . A partir do resumo:

Este artigo descreve um algoritmo para calcular o envelope de um conjunto de pontos em um plano, que gera convexa, não convexa cascos que representam a área ocupada pelos pontos indicados. O algoritmo proposto é baseado em uma aproximação K-neighbors mais próxima, onde o valor de k, o único parâmetro algoritmo, é usado para controlar a “suavidade” da solução final.

porque vou aplicar este algoritmo à informação geográfica, algumas alterações tiveram de ser feitas, nomeadamente ao calcular ângulos e distâncias . Mas estes não alteram de forma alguma o gist do algoritmo, que pode ser amplamente descrito pelos seguintes passos:

1. Encontre o ponto com a coordenada y (latitude) mais baixa e torne-o o ponto actual.
2. Encontre os pontos k mais próximos do ponto actual.
3. a partir dos pontos k mais próximos, seleccione o que corresponde à maior curva à direita do ângulo anterior. Aqui vamos usar o conceito de rolamento e começar com um ângulo de 270 graus (para Oeste).
4. verifique se, adicionando o novo ponto à cadeia de linhas em crescimento, não se intersecta a si próprio. Se o fizer, seleccione outro ponto do k-Mais Próximo ou reinicie com um valor maior de k.
5. faça do novo ponto o ponto actual e remova-o da lista.
6. após K iterações adicione o primeiro ponto de volta à lista.
7. ciclo para o número 2.

o algoritmo parece ser bastante simples, mas há uma série de detalhes que devem ser atendidos, especialmente porque estamos lidando com Coordenadas geográficas. As distâncias e os ângulos são medidos de uma forma diferente.

o código

aqui estou publicando é uma versão adaptada do Código do artigo anterior. Você ainda vai encontrar o mesmo código de agrupamento e o mesmo gerador de clusters em forma de nuvem. A versão atualizada agora contém um pacote chamado geomath.hullsonde você pode encontrar a classe ConcaveHull. Para criar seus cascos côncavos fazer o seguinte:

No código acima, points é uma matriz de dimensões (N, 2), onde as linhas contêm os pontos observados e as colunas contêm as coordenadas geográficas (longitude, latitude). A matriz resultante tem exatamente a mesma estrutura, mas contém apenas os pontos que pertencem na forma poligonal do aglomerado. Um tipo de filtro.

porque vamos lidar com arrays é natural trazer NumPy para a briga. Todos os cálculos foram devidamente vetorizados, sempre que possível, e foram feitos esforços para melhorar o desempenho ao adicionar e remover itens de arrays (spoiler: eles não são movidos de todo). Uma das melhorias em falta é a paralelização de código. Mas isso pode esperar.

organizei o código em torno do algoritmo como exposto no artigo, embora algumas otimizações foram feitas durante a tradução. O algoritmo é construído em torno de uma série de sub-rotinas que são claramente identificadas pelo jornal, então vamos tirá-las do caminho agora mesmo. Para o seu conforto de leitura, usarei os mesmos nomes usados no papel.

CleanList —limpeza a lista de pontos é realizada no construtor da classe:

Como você pode ver, a lista de pontos é implementado como um NumPy matriz por motivos de desempenho. A limpeza da lista é realizada na linha 10, onde apenas os pontos únicos são mantidos. O conjunto de dados é organizado com observações em linhas e coordenadas geográficas nas duas colunas. Note que eu também estou criando um array booleano na linha 13 que será usado para indexar no array do conjunto de dados principal, facilitando a tarefa de remover itens e, de vez em quando, adicionando-os de volta. Eu vi esta técnica chamada “máscara” na documentação NumPy, e é muito poderosa. Quanto aos números primos, discuti-los-ei mais tarde.

FindMinYPoint — Isso requer uma pequena função:

Esta função é chamada com o conjunto de dados de matriz como argumento e retorna o índice do ponto com a menor latitude. Note que as linhas são codificadas com longitude na primeira coluna e latitude na segunda.

RemovePointindices array. Este array é usado para armazenar os índices ativos na matriz do conjunto de dados principal, assim removendo itens do conjunto de dados é um snap.

embora o algoritmo descrito no Artigo exija a adição de um ponto à matriz que compõe o casco, este é realmente implementado como:

mais tarde, a variável test_hull será atribuída de volta a hull quando a cadeia de linhas é considerada como não-intersecção. Mas estou a adiantar-me ao jogo. Remover um ponto do conjunto de dados é tão simples como:

self.indices = False

adicioná-lo de volta é apenas uma questão de inverter o valor do array no mesmo índice de volta ao true. Mas toda essa conveniência vem com o preço de ter que manter nossas contas nos índices. Mais sobre isto mais tarde.

NearestPoints — as Coisas começam a ficar interessantes aqui porque não estamos lidando com coordenadas planar, então, com Pitágoras e com Haversine:

Note que o segundo e terceiro parâmetros são matrizes no conjunto de dados formato: longitude na primeira coluna e a latitude, no segundo. Como você pode ver, esta função retorna uma matriz de distâncias em metros entre o ponto no segundo argumento e os pontos no terceiro argumento. Assim que tivermos isto, podemos apanhar os vizinhos mais próximos da maneira mais fácil. Mas há uma função especializada para isso e merece algumas explicações:

a função começa criando um array com os índices de base. Estes são os índices dos pontos que não foram removidos do conjunto de dados array. Por exemplo, se em um aglomerado de dez pontos começássemos por remover o primeiro ponto, a matriz de índices de base seria . Em seguida, calculamos as distâncias e ordenamos os índices de matriz resultantes. O primeiro k é extraído e então usado como uma máscara para recuperar os índices de base. Está meio contorcido, mas funciona. Como você pode ver, a função não retorna um array de coordenadas, mas um array de índices no array do conjunto de dados.

SortByAngle — há mais problemas aqui porque não estamos calculando ângulos simples, estamos calculando rolamentos. Estes são medidos como zero graus a norte, com os ângulos a aumentar no sentido horário. Aqui está o núcleo do código que calcula os rolamentos:

A função retorna uma matriz de rolamentos medido a partir do ponto de cujo índice é o primeiro argumento, para os pontos cujos índices estão no terceiro argumento. A ordenação é simples:

neste ponto, a lista de candidatos contém os índices dos pontos k-mais próximos ordenados por ordem decrescente de rolamento.

IntersectQ-em vez de rodar as minhas próprias funções de intersecção de linhas, virei-me de forma para pedir ajuda. Na verdade, ao construir o polígono, estamos essencialmente lidando com uma cadeia de linhas, adicionando segmentos que não se cruzam com os anteriores. Testando para isso é simples: nós pegamos o array do casco em construção, convertemo-lo em um objeto de cadeia de linhas de forma, e testamos se é simples (não auto-intersecção) ou não.

em poucas palavras, uma cadeia de linhas torneadas torna-se complexa se se auto-cruza, de modo que o predicado is_simple se torna falso. Facil.

PointInPolygon-este acabou por ser o mais difícil de implementar. Permita-me explicar olhando para o código que realiza a validação do polígono final do casco (verifica se todos os pontos do aglomerado estão incluídos no polígono):

as funções de Shapely para testar a intersecção e inclusão deveriam ter sido suficientes para verificar se o polígono final do casco sobrepõe todos os pontos do aglomerado, mas não foram. Por quê? A forma é coordenada-agnóstica, então ela irá lidar com Coordenadas geográficas expressas em latitudes e longitudes exatamente da mesma forma que coordenadas em um plano cartesiano. Mas o mundo se comporta de forma diferente quando você vive em uma esfera, e ângulos (ou rolamentos) não são constantes ao longo de um geodésico. O exemplo de referência de uma linha geodésica que liga Bagdá a Osaka ilustra perfeitamente isso. Acontece que em algumas circunstâncias, o algoritmo pode incluir um ponto baseado no critério de rolamento, mas mais tarde, usando algoritmos planares de Shapely, ser considerado sempre um pouco fora do polígono. Isso é o que a pequena correção de distância está fazendo lá.Levei algum tempo a perceber isto. A minha ajuda de depuração foi QGIS, uma grande peça de software livre. Em cada passo dos cálculos suspeitos eu iria enviar os dados em formato WKT em um arquivo CSV para ser lido como uma camada. Um verdadeiro salva-vidas!

finalmente, se o polígono não conseguir cobrir todos os pontos do aglomerado, a única opção é aumentar k e tentar novamente. Aqui eu adicionei um pouco da minha própria intuição.

Prime k

the article suggests that the value of k be increased by one and that the algorithm is executed again from scratch. Os meus primeiros testes com esta opção não foram muito satisfatórios: os tempos de execução seriam lentos em aglomerados problemáticos. Isto foi devido ao aumento lento de k, então eu decidi usar outro esquema de aumento: uma tabela de números primos. O algoritmo já começa com k = 3, Então foi uma extensão fácil para fazê-lo evoluir em uma lista de números primos. Isto é o que você vê acontecendo na chamada recursiva:

eu tenho uma coisa para números primos, você sabe…

Explodir

O côncavo casco polígonos gerados por este algoritmo ainda precisa de algum processamento adicional, porque eles só discriminar pontos dentro do casco, mas não próximo a ele. A solução é adicionar algum enchimento a estes grupos magros. Aqui eu estou usando a mesma técnica usada antes, e aqui está o que parece:

Aqui, eu costumava bem Torneadas do bufferfunção para fazer o truque.

a função aceita um polígono de forma e retorna uma versão inflada de si mesma. O segundo parâmetro é o raio em metros do revestimento adicionado.

executando o código

comece por puxar o código do repositório GitHub para a sua máquina local. O ficheiro que deseja executar está ShowHotSpots.py na pasta principal. Após a primeira execução, o código será lido nos dados de acidentes de tráfego do Reino Unido de 2013 a 2016 e agrupá-lo. Os resultados são então Cache como um arquivo CSV para as corridas subsequentes.

você será então apresentado com dois mapas: o primeiro é gerado usando os clusters em forma de nuvem, enquanto o segundo usa o algoritmo de clustering côncavo discutido aqui. Enquanto o código de geração de polígonos executa, você pode ver que algumas falhas são relatadas. Para ajudar a entender por que o algoritmo falha em criar um casco côncavo, o código escreve os clusters para arquivos CSV para o diretório data/out/failed/. Como de costume, você pode usar o QGIS para importar esses arquivos como camadas.

Essentially this algorithm fails when it does not find enough points to “go around” the shape without self-intersecting. Isto significa que você deve estar pronto para descartar estes aglomerados, ou para aplicar um tratamento diferente para eles (casco convexo ou bolhas fundidas).

Concaveness

It’s a wrap. Neste artigo I apresentou um método de pós-processamento DBSCAN gerou clusters geográficos em formas côncavas. Este método pode fornecer um melhor ajuste fora do polígono para os clusters em comparação com outras alternativas.

Obrigado por ler e desfrutar de remendar com o código!

Kryszkiewicz M., Lasek P. (2010) TI-DBSCAN: Clustering with DBSCAN by Means of the Triangle Inequality. In: Szczuka M., Kryszkiewicz M., Ramanna S., Jensen R., Hu Q. (eds) Rough Sets and Current Trends in Computing. RSCTC 2010. Lecture Notes in Computer Science, vol 6086. Springer, Berlin, Heidelberg

Scikit-learn: Machine Learning in Python, Pedregosa et al., JMLR 12, pp. 2825-2830, 2011

Moreira, A. e Santos, M. Y., 2007, Côncavo do Casco: Uma K-vizinhos mais próximos abordagem para o cálculo da região ocupada por um conjunto de pontos

Calcular a distância, azimute e mais entre a Latitude/Longitude de pontos

repositório no GitHub