As tecnologias LiDAR, do acrónimo em inglês Light Detection And Ranging são comummente usadas como pulsos de luz laser emitidas (e recebidas) por um plataforma aerotransportadora (avião, helicóptero, drone), por plataforma móvel (comboio, automóvel) ou em suporte fixo no solo.
Estes levantamentos tanto são feitos por varrimentos lineares como circulares, entre outros. Não deixando de ser uma tecnologia de Deteção Remota, terá sempre erros associados, que cabe ao engenheiro minimizar e acima de tudo quantificar.
Não sendo o contexto teórico o motivo desta página, partilharei alguns métodos de tratamento de nuvens de pontos, recorrendo a softwares como o Lastools e CloudCompare.
Para os exemplos usados, recorrerei aos levantamentos do projeto áGIL do ICNF, possíveis de descarregar aqui. Estes levantamentos fazem parte duma fase experimental para apoio ao mundo técnico e científico, estando prevista uma grande empreitada para o país inteiro, disponibilizando EVENTUALMENTE de modo gratuito esses dados a todos os interessados.
Como filho das Beiras, usarei a Serra da Lousã como cobaia.
Primeiro, um mapa do enquadramento da Serra da Lousã com os seus tiles do levantamento LiDAR:
Cada tile destes contem uma nuvem de pontos que deverá ser processada com o ultimate end de produzir um mega e altamente detalhado Modelo Digital de Terreno (MDT) da área toda, onde 13,1 GB de informação XYZ ("+ IVA") deverá dar um único ficheiro TIFF.
Como será do conhecimento do público técnico, existem variadíssimos MDTs disponíveis, que na verdade alguns se tratam de MDT e outros de Modelos Digitais de Superfície (MDS), abundando estes últimos. Falamos duns globais, de outros americanos e outros europeus, por exemplo do SRTM, ASTER, EU-DEM, ALOS-PALSAR, entre outros.
Para nos certificarmos que estamos a trabalhar só com pontos de solo (bareground) requer classificar os pontos, pelo que usar-se-á a linha de comando para correr os processos do Lastools.
Primeiro passo: simplesmente escolhemos um tile para tratar, que irá ser o que contém o belíssimo Castelo da Lousã, embora que possa aparecer cortado por se encontrar no extremo Norte do tile. Por consulta ao QGis:
Abrindo a nuvem de pontos no CloudCompare, em falsa cor, facilmente se reconhece o castelo e a ermida na parte superior da nuvem. Obviamente que se tratando dum levantamento aéreo, não se irá verificar, ou muito dificilmente se terão os pontos das paredes, aparecendo só os telhados ou superfícies lisas a flutuar no espaço tridimensional.
Rapidamente também nos apercebemos pela orientação das "linhas" o sentido do varrimento, o que facilmente subentende que são varrimentos paralelos no mesmo sentido.
Trabalhar-se-á este tile inteiro diretamente na linha de comandos, assim:
Tratando-se de mais de 3 milhões e 600 mil de pontos, recorrendo ao CloudCompare mais uma vez para ler a nuvem classificada, apresenta-se da seguinte maneira:
Todos os pontos azuis são os pontos objeto - aqueles que influenciariam a criação de um MDS - e a verde os pontos bareground, que são os que interessam para a realização de um MDT - ou MDE, como Modelo Digital de Elevação.
Tendo uma densidade de cerca de 10 pulsos/m2, perfaz que se possa esperar uma resolução de cerca de 30 centímetros por pixel no MDT.
Outro exercício que se poderia fazer por análise em detalhe (cada nuvem é um caso isolado) seria separar em diferentes classes o que é vegetação baixa, da vegetação alta, dos edifícios, das linhas elétricas (caso as existam), ou mesmo dos pontos aberrantes, o chamado "lixo" que não queremos. Sendo de interesse só os pontos de chão, passamos à frente.
Esse procedimento poderia ser feito pela ferramenta lasgrid do Lastools, mas depressa nos aperceberíamos que nesta área tão densamente povoada por árvores e, em geral, objetos, por muito que se lhe definisse um pixel (step) de 30 cm, expandindo cada ponto também com a mesma dimensão ou mesmo que maior, iria gerar um raster com "buracos", correspondendo a pixels NoData, o que não nos interessa. Como o objetivo é ter uma resolução de 30 cm, este método não funcionará de todo.
Resultado dos pontos bareground usando a ferramenta lasgrid
Esta ferramenta tem aplicabilidade correta se estivéssemos a trabalhar com os pontos-objeto para modelar, por exemplo, um modelo de copas das árvores, comummente conhecido por Canopy Height Model (CHM).
Usando a ferramenta las2dem, permite interpolar os pixels ditos NoData e preencher esses buracos com pixels interpolados entre os vizinhos mais próximos. Continua-se a falar de pontos bareground para estes procedimentos.
O mesmo procedimento, feito usando a ferramenta las2dem.
Como se pode verificar, nos extremos do tile encontra-se partes inacabadas, cantos cortados e edges incompletos. Como se pretende construir um MDT de áreas maiores, usar-se-á a mesma ferramenta para fazer merge a mais tiles de modo a corrigir essas imperfeições. Eis um exemplo da nuvem de pontos ao canto superior esquerdo deste tile:
Como se pode verificar, o las2dem não interpolou para lá dos pontos bare ground (verdes) na diagonal até ao canto superior esquerdo da nuvem de pontos. Fazendo um seccionamento àquele canto e ponto em vista de perfil, verifica-se que os pontos bareground ora bem identificados, justificam a dificuldade em encontrar forma de interpolar naquele canto, no extremo esquerdo do perfil, quando se avista arvoredo que deverá ultrapassar os 15-20 metros de altura.
Para determinar as alturas do arvoredo e modelar variados processos como o já referido CHM, usarei outro post para pôr isso em detalhe, podendo ser feito também através de uma normalização da nuvem de pontos.
Abrindo a área de tiles para 2x2 contendo o já anterior definido, irei fazer o merge das 4 nuvens de pontos para fazer um MDT dos 4 e, assim esperando, colmatar as falhas dos cantos e dos edges. Se usasse uma dimensão de pixel maior, muito provavelmente estas deformações seriam mitigadas, mas eventualmente nunca ultrapassadas.
Usando mais uma vez a linha de comandos para processar as outras 3 nuvens de pontos, obtém-se os 4 ficheiros classificados necessários para o passo seguinte.
Rapidamente se consegue observar quais são os pontos bareground neste universo de pontos de 4 nuvens de pontos postas lado a lado, que à escala de visualização se destacam logo as estradas e caminhos florestais, e uma lavrada para plantação de alguma espécie de árvores (👀).
Correndo o las2dem para as 4 nuvens:
Obtém-se um MDT de alta resolução para os 4 tiles, continuando a haver apenas imperfeições nos seus 4 cantos e nalguns edges, facto que irá ser ultrapassado à medida que mais área se vai acrescentando.
Obtém-se um MDT de alta resolução, sem pixéis calcados nas zonas de união de tiles, perfazendo um total de 40 MB, com que podemos trabalhar para os mais diversos fins, desde extrair curvas de nível, calcular declives, obter volumes, ou até comparar altitudes com outros MDT, verdadeiros modelos de terreno, e não de superfície.
----------------------------------------------------------------------------------
CONTINUAÇÃO
DESENVOLVIMENTO DO MDT PARA AS NUVENS DE PONTOS TODAS
Para os próximos passos, tratou-se de correr o processo de classificação de todas as nuvens de pontos. Desenvolveu-se um código Python que permitisse ler os nomes das nuvens originais e que lhes desse um output - numa lista - de novos nomes, que poderiam eventualmente ser usados no lasground, lasheight e las2dem como uma massiva lista corrida de inputs e outputs, que seria sucedida dos parâmetros mas rapidamente se explorou a documentação e se verificou que há forma mais simples de o fazer por parâmetros das próprias ferramentas. Assim, usando um comando da linha de comandos, exporta-se um simples ficheiro .txt que pode ser chamado à linha de comandos quando se correr qualquer uma das ferramentas.
Usando este método, evita-se ter que recorrer a este método, que comprova-se que o lasground não consegue processar mais que um input de cada vez, e sucessivamente mais que um output. Este exemplo serviria para renomear de antemão os ficheiros .laz já classificados.
Assim, recorrendo à linha de comandos de novo de modo a correr com a lista de ficheiros .laz de todas as nuvens de pontos, corre-se o seguinte:
Como se pode verificar nos parâmetros, para alem dos outputs irem ter todos o mesmo nome do original, mas não deixarão de ser em formato .laz (caso omitisse iria dar ficheiros .las, muito mais pesados, lento de produzir e com a mesma informação), acrescentou-se outros parâmetros que permitem que o programa refine a classificação, e não acontecer casos como telhados, terraços e outras superfícies lisas serem confundidas com pontos bare ground.
Este processo demorou 2h41m35s a correr para os 845 tiles, durando entre cerca de 8 a 14 segundos por tile, mais tempo de escrita em disco [cerca de 2 segundos: criação (cópia original) + modificação com o processo].
Depois do processo corrido, foi principal preocupação analisar a classificação um a um nos tiles dos cantos, alguns bem pequeninos e nada percetíveis na primeira figura. Essa análise começou por ser feita começando pelos avisos de não conseguir classificar, e pelos tamanhos de ficheiros no explorador. Após devida triagem, excluíram-se os seguintes:
Os ficheiros das nuvens classificadas, agora uma amostra de 821 nuvens, enunciam-se pelo sufixo "_1" no final do nome, igual ao original.
Estas nuvens, antes de avançar para o las2dem, devem, a bem, ser novamente analisadas, de modo a limpar eventuais spikes de pontos aberrantes, e pô-los à parte, em classificação própria. Este passo podia ser só ignorado, mas pelo bem dum MDT de qualidade, mais vale prevenir que remediar.
Do mesmo jeito que no anterior, corre-se o comando de criar uma lista em .txt com os ficheiros a ser usados, e (re)classificam-se os pontos de terreno entre os -20 cm e os 10 cm de altura e todos os valores abaixo dos -20 cm e acima dos 35 metros como sendo pontos aberrantes. Este valor de 35 metros é para efeitos de teste, pois em nuvens particulares deveria ser estendido para alturas maiores, uma vez que este valor diz respeito à altura em Z em cada nuvem, e em terrenos mais escarpados os pontos mais altos dos objetos ultrapassam esta altura. Existe algumas nuvens que contêm muito "lixo" proveniente do levantamento LiDAR original, como manchas a flutuar ou colunas inteiras de pontos a trespassar o terreno, para alturas muito altas ou muito baixas. Boa política seria antes, analisar nuvem a nuvem, e eliminar manualmente esses pontos demasiado aberrantes. São 845 tiles, irei deixar para fins de teste assim como está, e na eventualidade de declives ou MDT com valores "astronómicos", refinarei.
Exemplo de comparação entre a mesma nuvem de pontos antes e depois do lasheight:
Na primeira nuvem temos o processo original de classificação, estando a verde os pontos bareground e na segunda os pontos classificados como bareground (a azul), objetos (branco) e resíduos (laranja), que são para ignorar no caso da realização do MDT.
Este processo de classificação das 821 nuvens demorou 1h25m24s a concluir, demorando entre 2 e 12 segundos por nuvem (algumas nuvens, frações de segundos), e criou outro sufixo "_1" à frente do já existente. Estas nuvens são as que irão ser chamadas, por fim, para a criação do MDT final com 30 cm de resolução.
Como não se pretende classificar ao detalhe o que são árvores, casas (telhados) e outros objetos (como as torres eólicas), não irá ser feito um refinamento a esses objetos usando o lasclassify. Poderei fazê-lo em post próprio.
CRIAÇÃO DO MDT
UTILIZANDO O PLUGIN BLAST2DEM
Como anunciado anteriormente, para a criação do MDT terá que se fazer uma análise a cada nuvem de pontos classificada, colocá-las todas em cache e só depois começará o processamento. Das 821 nuvens de pontos, rapidamente se concluiu por tentativa e erro que há nuvens que induzem problemas à criação do raster final, tendo sido excluídas mais 4 nuvens, perfazendo 817 nuvens a ser computadas.
A ferramenta las2dem apresenta limitações ao número de pontos que consegue processar, logo exige que se use o plugin do Lastools, blast2dem, que faz exatamente a mesma coisa, mas consegue processar até 2 biliões de pontos; os pontos todos destas nuvens, ultrapassam em muito as capacidades do las2dem.
Correndo a ferramenta blast2dem, definindo-lhe o sistema de coordenadas, um step de 30 cm, usando triangulação para distâncias até 200 metros e usando todas as 817 nuvens, temos como resultado final o seguinte MDT:
[CONTINUA]
[Próximo passo avaliar a exatidão do levantamento LiDAR recorrendo às cotas de base dos vértices geodésicos, usando a Raíz Quadrada do Erro Médio - RMSE]
[Ainda próximo passo: criar as curvas de nível metro a metro e de 10 em 10 metros. Irei disponibilizar as curvas de 25 metros para consulta por parte dos interessados que assim se manifestem]
[esta página irá ser continuada à medida que mais desenvolvimentos surgirem]