Colegas pesquisadores e entusiastas de SIG e código-livre,

Bem-vindos ao meu blog!

Gostaria de começar com um aviso - posso ser pesquisadora desta área, mas isso não significa que tudo o que faço ou escrevo aqui funcionará para você, em suas próprias configurações de desktop e versões de packages. Me eximo de responsabilidade se você perder dados ou bagunçar sua instalação. Eu também não autorizo nenhum tipo de cópia do meu conteúdo, exceto dos gists, fornecidos com licença do tipo MIT.

Hoje, eu vou apresentar para vocês um gist (fragmento de código) que desenvolvi, para separar as informações de máscaras contidas dentro de uma Fmask em diferentes bandas, com o objetivo de visualizá-las separadamente. Hoje, nós vamos visualizá-las no QGIS, usando o Python Console do QGIS, também conhecido como PyQGIS.

Esta semana, reencontrei em um conjunto de dados de Sensoriamento Remoto, o Fmask. É uma camada de máscara, que geralmente vem como uma das bandas ou rasters fornecidos com o bloco original do produto de Sensoriamento Remoto que você está baixando. Porém, este Geotiff utiliza (de forma muito inteligente) todos os bits de um byte. Cada bit de cada byte em cada pixel é uma classificação diferente, que pode ser posteriormente utilizada, como máscara, por nós Engenheiros ou Pesquisadores de Recursos Hídricos ou de Ciências Ambientais.

Confesso que acho essa forma de guardar as camadas de máscaras muito mais genial do que irritante. Análogo ao que penso sobre o formato de arquivo grib2. E se eu e outras pessoas achamos difícil de ler, bem, nós estamos todos errados, porque esse formato economiza muito espaço! Já imaginou ter que baixar outros 6 rasters para cada tile, só para poder ler as máscaras? Que bom que não precisa, né?

Ok, mas como exatamente são armazenados os dados dentro da F Mask?

Cada informação (não relacionada às outras) é salva em um bit dentro do byte. 0 geralmente significa ausência de tal elemento e 1 geralmente significa presença. Por exemplo, se um determinado bit significa “neve/gelo” e seu valor em um determinado pixel é 1, isso significa que esse pixel tem cobertura de neve ou gelo. Isso é relevante para investigações de sensoriamento remoto, pois, para muitos tipos de análises, esse pixel precisará ser descartado, pois seu valor de refletância será artificialmente muito alto para as bandas utilizadas.

Vejamos, especificamente, o Fmask fornecido com o produto Harmonized Landsat Sentinel, da NASA. O LP DAAC da NASA montou um tutorial sobre o processamento do produto HLS.

Muitos outros produtos da NASA fornecem camadas Fmask ou outras máscaras incorporadas do tipo.

Para o HLS, a Tabela 10 do guia do usuário do HLS nos fornece informações sobre o conteúdo de cada bit na camada Fmask.

Um pouco de teoria

Um byte é formado por 8 bits, numerados de 0 a 7. Cada bit é binário, ou seja, pode ser preenchido apenas com 0 ou 1. Isso permite a representação de 256 valores únicos ( $2^8$ ) por cada byte. Claro, você pode ter diferentes maneiras de armazenar informações numéricas em um byte, como floats e assim por diante, mas esse não é o meu ponto aqui. A questão é que é assim que a maioria dos softwares leria um byte que representa números inteiros, armazenados em um byte.

É assim que o QGIS leria (e lê) seus bytes que compõem os pixels em Fmasks. Por padrão, ele abre o Geotiff do Fmask e mostra o valor composto das máscaras fornecidas, de uma forma que não é necessariamente amigável, a menos que você seja muito bom em fatoração matemática.

Fatoração?

Sim, fatoração. Se você ver que um determinado pixel de interesse tem um valor de 192 na F Mask, o que isso significa? Você saberia, apenas olhando para ele, quais máscaras estão ativas, o que foi sinalizado pelo Fmask e o que não foi? Eu não saberia.

Uma maneira fácil de descobrir quais são as máscaras ativas, seria fatorar o valor 192 em uma base 2. Como aprendemos no colégio:

$2^7=128$

é a maior potência de 2 que pode ser subtraída de 192. Portanto, devemos fazer nosso número menos 128 e refazer esta operação para a próxima maior potência de 2 possível.

$ 192 - 128 = $ 64

Neste caso, é $ 2 ^ 6 $.

$2^6=64$

$ 64 - 64 = 0 $

O resto já é zero, então não precisamos fazer nenhuma operação.

Então, 192 é

$1\cdot128+1\cdot64$

Isso significa que os bits sinalizados são 128 e 64. 128 é $2^7$ e 64 é $2^6$ . Então esses são os bits nas posições 6 e 7. Eles são 1, enquanto os outros bits são 0. Esses dois bits representam, segundo a tabela HLS, a classificação da qualidade do aerossol e, juntos, significam “Alta” Qualidade do Aerossol.

Ok, mas isso não foi difícil. Diga-me novamente por que você fez um código em Python para fazer isso.

Claro! Conseguimos, fatoramos 192 na base 2! Mas vamos imaginar um outro cenário, o raster que você baixou agora parece ter 10 combinações de bits únicas, cada uma resultando em um número diferente que você terá que fatorar. E acontece que você só quer dar uma olhada na máscara que representa cobertura de corpos d’água! Como você pode separar visualmente a máscara de água das outras coisas que estão acontecendo ali?

Se você prefere fazer cálculos na sua cabeça o tempo todo, pode continuar fazendo isso. Neste post, estou tentando fornecer uma alternativa mais simples.

Como funciona o meu código em Python (PyQGIS)?

Ele faz a fatoração na base 2 dentro dele, e então salva cada informação em uma banda separada de um tif raster. Ele foi feito para o produto HLS, mas sinta-se à vontade para atualizá-lo para seu dataset preferido.

Estou aqui pelo código do PyQGIS

Este é o gist, hospedado no GitHub.

Como rodar o código

Escrevi algumas instruções abaixo sobre como executar o gist.

Defina a localização do raster do Fmask

Rasterloc = 'localizacao em disco da Fmask original'

Este é um exemplo de Fmask retirado de um tile HLS.

Defina o nome com o qual deseja salvar seu arquivo final

Saveloc ='localizacao onde o Fmask preprocessado deve ser salvo'

Execute o gist no Console Python do QGIS

Chame a função

New_layer=Fmask_bit_band_separation(rasterloc,saveloc)

Ajuste a simbologia para ver as máscaras separadamente ou para compô-las 3 a 3

Depois de chamar a função no PyQGIS, você notará uma nova camada em seu painel de camadas do QGIS.

Para o nosso exemplo, ela se parece com isso:

Por ser um raster com múltiplas bandas, sua visualização padrão é RGB, cada cor (vermelho, verde, e azul) representando uma das 3 primeiras bandas. Neste caso, as três primeiras bandas são zero.

Podemos alterar as cores representadas pelo RGB para tornar a visualização mais nítida. Neste caso, adicionei as bandas 6 e 7 à visualização.

Às vezes é melhor mudar a simbologia para uma única banda para poder ver uma única máscara, separada das outras. Se quisermos ver a máscara de água, é a banda número 6.

Áreas com valor 1 têm Água, e áreas com valor 0 não. Lembrando que esta máscara é binária e não representa “graus de existência de água”. Nesta imagem, se vê claramente a existência de um córrego que corta o raster lado-a-lado.

Pronto!

Extras

• Existe um tutorial do LP DAAC que ensina como ler dados de uma Fmask em Python. Por que você não está usando o código deles aqui?

Esse é um tutorial muito útil e me ajudou a lidar com a Fmask em um projeto da NASA em que estou colaborando. No entanto, a parte de QA deste tutorial ensina como mascarar áreas de HLS dentro de um código Python, mas não como salvar saídas individuais das diferentes bandas, nem como abrir o raster resultante automaticamente no QGIS. Poderia servir de base, mas optei por fazer os cálculos do meu jeito.

• Por que meu raster da minha Fmask final é muito maior que o raster da original?

Isso ocorre porque o formato original da Fmask é otimizado para usar a quantidade mínima possível de armazenamento. Quando separamos estas máscaras individuais que estavam cada uma em um bit e salvamos cada uma delas em uma banda diferente, estamos salvando cada um destes (originalmente) bits em um byte. Isso é, de certa forma, um desperdício de espaço. Tem um propósito, que é ajudar a nós, pesquisadores, engenheiros, analistas de dados e programadores de SIG a entender os dados. Por outro lado, é sim um desperdício do ponto de vista computacional.