Tempo de leitura: 9 min

Neste artigo regresso a um assunto já familiar neste blog – criar mosaicos de ortofotomapas usando o GDAL – (sim eu sei, outra vez?) mas como tenho andado às voltas com as áreas sem informação, que surgem negras nos mosaicos pensei em postar o que acabei por fazer. A solução final é usar máscaras, e não bandas alfa como habitual. Vamos ver como e porquê… assume-se já alguma familiaridade com o GDAL, mas pode sempre saltar as partes teóricas aborrecidas e ver os comandos usados no final do artigo 😉

O Problema

Então a situação é esta: temos uma colecção de 279 ortos que cobrem uma área irregular, não rectangular:

O que acontecerá se criarmos um mosaico? Todas as partes em branco ficarão negras porque não existe informação para essas áreas.

Comprimir

Como cada ortofotomapa tem as 3 bandas usuais, rgb, e ocupam entre 292MB e 465MB cada, temos um total de 122GB. A primeira coisa a fazer, claro, é comprimir…

O comando é sempre o mesmo…

for %I in (directoria_originais\*.tif) do gdal_translate -co photometric=ycbcr -co interleave=pixel -co tiled=yes -co compress=jpeg -co jpeg_quality=75 %I comprimidos\%~nI.tif

Resulta assim uma colecção de ortos comprimidos que usaremos daqui para a frente, e podemos arquivar os originais fora do disco. O espaço ocupado por estes ortos comprimidos passou para apenas 6,39GB (ainda sem pirâmides).

Criar um mosaico virtual

Grande parte do processamento usado neste artigo baseia-se no formato virtual do GDAL, com extensão .vrt. Um ficheiro .vrt é apenas um ficheiro xml que lista as imagens que queremos usar, e permite incorporar alguns processamentos que queremos aplicar, sem efectivamente os realizar adiando-os para outro momento. É por isso muito rápido de criar e testar. Só quando convertemos um ficheiro .vrt para um formato “real”, como tif, é que as operações são efectuadas, ocupando o processador e o espaço em disco necessários. Existem vários artigos que explicam muitas das operações comuns, como este tutorial, existindo também a referência do formato.

O nosso exemplo centra-se na criação de um  mosaico de imagens onde existem espaços sem ortos, e que por isso terá áreas negras nesses espaços…

Para criar o mosaico virtual da forma mais simples basta usar este comando:

gdalbuildvrt mosaico_simples.vrt comprimidos\*.tif

Este comando cria o tal ficheiro mosaico_simples.vrt que apenas lista todos os nossos ortos. A vantagem é que é identificado pelo QGIS e pelo ArcMap como uma imagem única. O aspecto é este:

E cá está o problema das áreas sem informação…

(claro que se tentarem isto sem criar pirâmides vai demorar um bocado!)

Áreas Negras em Mosaicos

Estas áreas resultam de 2 situações…

1. Áreas sem informação relativas a ortos que não foram adquiridos, porque simplesmente não existem ou porque a nossa área de interesse é irregular, não sendo necessário comprar esses ortos (é o caso em mãos). Geralmente, a área de um concelho não constitui um rectângulo perfeito, sendo apenas adquiridos os ortos que cobrem o concelho, resultando nestas áreas sem informação. Nestes casos, em geral, usam-se bandas alfa…que veremos mais adiante.
2. Áreas dentro de um ortofotomapa que não foram cobertas pelo vôo ou que não foram processadas na ortorrectificação ou mesmo até retiradas propositadamente do produto final. Por exemplo, um ortofotomapa junto à fronteira com Espanha terá, tipicamente, áreas negras no território espanhol. O nosso exemplo também tem disto.

Uma solução possível nestes casos é usar o comando nearblack do GDAL, para cada imagem, que consegue detectar as áreas “exteriores” de uma imagem de cor preta ou similar, tornando-as brancas. Note-se que os pixéis continuam a ser desenhados, apenas passam a ser brancos. Como não gosto desta solução, passamos à frente (porque as áreas brancas continuam opacas e é preciso tratar cada imagem individualmente).

A solução da Banda Alfa

Ao criar um mosaico, as áreas onde faltam ortos ficarão negras. A solução passa por informar o GDAL que as áreas sem informação devem ser consideradas transparentes, criando-se usualmente uma 4ª banda “alfa”, onde áreas com informação são 100% opacas, e por isso visíveis, enquanto que as áreas sem informação são 100% transparentes, e por isso invisíveis. Assim, se repetirmos o primeiro comando de criarmos um mosaico virtual com o comando gdalbuildvrt, e adicionar a opção “–addalpha”, temos uma forma espectacularmente simples de resolver a questão:

gdalbuildvrt –addalpha mosaicovirtualalfa.vrt comprimidos\*.tif

Este comando cria o ficheiro mosaicovirtualalfa.vrt que é, novamente, apenas a listagem de todos os ortos que estavam na pasta indicada, mas adicionando também uma nova banda fictícia à imagem que indica ao QGIS/ArcMap que as áreas sem ortos devem ser 100% transparentes. Isto significa que o nosso mosaico passa a ser uma imagem de 4 bandas: rgb+alfa ou “rgba”.

Na imagem abaixo, o espaço entre os 2 ortos estaria preto, mas é transparente, vendo-se o fundo azul-claro definido para o mapa.

Podemos ver que a banda alfa é identificada pelo GDAL usando o comando gdalinfo, que nos  indica para cada banda a expressão ”PER_DATASET_ALPHA” e também a 4ª banda com interpretação de “Alpha”:

No entanto, esta técnica tem 2 problemas: não detecta as áreas negras (sem informação) dentro dos ortos, e tem também problemas em comprimir com jpeg resultando em imagens e pirâmides muito maiores. Nos meus testes, pirâmides para imagens com 4ª banda alfa ocuparam 3x mais espaço em disco…

Comprimir imagens com 4ª banda alfa

Um pequeno parêntesis… A banda alfa é muito útil, mas complica a compressão de imagens de 3 bandas mais comuns – rgb. Isto sucede porque a melhor compressão deste tipo de imagem é conseguida usando um esquema de cor denominado YCbCr, e que não funciona quando existe uma 4ª banda alfa, ou seja, quando temos imagens rgba como o nosso mosaico acima… isto leva a que o mosaico quando comprimido ocupe 2x a 3x mais espaço em disco. Por outro lado, se comprimirmos uma imagem rgba com jpeg, sucede que poderão ocorrer imprecisões nos limites da transparência, o que não é nada desejável.

Por exemplo, se tentarmos este tipo de compressão numa imagem rgba sucede este erro:

gdal_translate -co compress=jpeg -co jpeg_quality=20 -co interleave=pixel -co photometric=ycbcr teste.vrt mosaico.tif

Input file size is 24000, 100000

ERROR 6: PHOTOMETRIC=YCBCR requires a source raster with only 3 bands (RGB)

Assim, o que fazer?? Podemos sempre fazer a compressão sem a opção YCbCr, ocupando 3x mais espaço em disco, mas felizmente temos outras opções…

Ficheiros de Máscara

Existe outra opção para a definição de áreas transparentes, muito usada na área gráfica, e que é usar ficheiros de máscara, ou seja, ficheiros que indicam se um pixel é ou não visível. Cada pixel na máscara tem apenas o valor 0 ou o valor 255, sendo por isso altamente compressível num esquema não-degradante como é a compressão zip, apresentando assim várias vantagens sobre as bandas alfa:

• é um ficheiro separado, não influenciando a compressão da imagem original;
• é altamente compressível num esquema zip (deflate na nomenclatura GDAL), pelo que ocupa muito pouco espaço, menos até que uma banda alfa comprimida em jpeg;
• mantém a precisão na definição das áreas transparentes mesmo quando comprimido;
• e, uma vantagem menos óbvia, é reutilizável, podendo o mesmo ficheiro de máscara ser aplicado a diversas imagens, bastando renomeá-lo.

Então como se criam? Aviso que são vários passos…

Criar um mosaico com máscara usando um shapefile de cutline

Em primeiro lugar, criamos um shapefile que delimita a área dos nossos ortofotmapas. Basta para isso criar um shapefile com a grelha dos ortos do nosso futuro mosaico e depois fazer o dissolve em QGIS/ArcMap.

Para criar o shapefile basta usar o comando gdaltindex. É muito fácil de usar:

gdaltindex quadricula.shp comprimidos\*.tif

Creating new index file…

Depois, faz-se o dissolve para termos apenas um polígono com todas as quadrículas…

Em seguida vamos usar o comando gdalwarp para criar um mosaico e usar o shapefile como o limite de corte, como numa operação Clip, criando a banda alfa neste processo. Como não queremos perder tempo nem espaço em disco, em vez de criar um ficheiro tif vamos criar um mosaico virtual .vrt:

gdalwarp -multi -wm 312 --config GDAL_CACHEMAX 312 -dstalpha -cutline shapes\index_diss.shp -of vrt mosaico_simples.vrt mosaico_cortado_alfa.vrt

A saber:

• -wm 312 –config GDAL_CACHEMAX 312 = parâmetros de uso de memória para acelerar o processo; (aliás, o GDAL_CACHEMAX devia ser obrigatório em todos os comandos!!)
• -multi = usar vários cores do processador se existirem;
• -dstalpha –cutline = usar o shapefile como área de corte, delimitando a área opaca na banda alfa, criando assim um mosaico rgba (4 bandas: rgb+alfa);

Este comando é instantâneo porque adia todo o processamento para quando se efectivamente converter este ficheiro mosaico_cortado_alfa.vrt num ficheiro tif.

Pode-se investigar o que está no ficheiro criado, mosaico_cortado_alfa.vrt usando o notepad. Podemos ver todos os nossos ortos listados, e uma última banda alfa definida pela operação cutline com o nosso shapefile.

O passo seguinte é converter este mosaico virtual de 4 bandas noutro mosaico virtual de 3 bandas, e mais uma máscara que substituirá a banda alfa. Este passo baralha um pouco o processo mas vale a pena… relembro que a shortlist de comandos está no fim.

Para isso, vamos usar o comando gdal_translate. O truque aqui é escolher as bandas que queremos usar… convertendo a 4ª banda alfa numa máscara:

gdal_translate -b 1 -b 2 -b 3 -mask 4 -co alpha=no -co photometric=ycbcr -of vrt mosaico_cortado_alfa.vrt mosaico_mask.vrt

Onde:

• -b 1 –b 2 –b 3 => usar as bandas 1, 2 e 3 como rgb
• -mask 4 => usar a 4ª banda para definir uma máscara
• -noalpha => certificar que deixamos de ter banda alfa
• -co photometric=ycbcr => certificar que o mosaico criado terá a maior taxa possível de compressibilidade

E mais uma vez obtemos um ficheiro .vrt que contém as instruções de processamento para criar uma imagem rgb+máscara.

Resta então o passo final! Converter este .vrt para um tif altamente comprimido, de 3 bandas, e com um ficheiro extra de extensão .msk.

O comando gdal_translate tem a seguinte expressão:

gdal_translate -co alpha=no -co photometric=ycbcr -co interleave=pixel -co tiled=yes -co compress=jpeg -co jpeg_quality=70 mosaico_mask.vrt mosaico_mask.tif

As opções já são conhecidas. E  obtemos o tal tif e o tal msk.

Resultado final

Para visualizarmos o resultado convém criar primeiro as pirâmides com o comando gdaladdo:

gdaladdo -ro -r average --config JPEG_QUALITY_OVERVIEW 70 --config COMPRESS_OVERVIEW JPEG --config INTERLEAVE_OVERVIEW PIXEL mosaico_mask.tif 2 4 8 16 32 64

Podemos agora abrir o mosaico_mask.tif no ArcMap e ver o resultado:

Espectacular! Recapitulando, os ficheiros finais e tamanhos respectivos são os seguintes:

Ficheiros	Tamanho
ortos originais	122 GB
ortos comprimidos	6,39 GB
1 mosaico virtual com máscara	333 kb
Pirâmides do mosaico virtual	1,65 GB
1 ficheiro máscara (.msk)	111 MB
Pirâmides da máscara	28,9 MB
1 mosaico final virtual para QGIS	3 kb

Total dos ficheiros comprimidos: ~8,18 GB, taxa de compressão de 93%!

Falta discutir o último ficheiro na lista… sucede que a máscara é detectada e usada automaticamente pelo ArcMap, mas não pelo QGIS… teremos ainda de fazer um último truque para convencer o QGIS a usar a máscara.

Visualizar máscaras em QGIS

Fazemos o inverso, e criamos uma banda alfa virtual a partir da máscara… não ocupando espaço em disco… eu sei, é confuso 😉

gdal_translate -of vrt -b 1 -b 2 -b 3 -b mask mosaico_final.tif mosaico_final_qgis.vrt

Desta forma, o QGIS já reconhece a transparência, bastando carregar este ficheiro mosaico_final_qgis.vrt.

Truque final: o QGIS detecta bandas alfa automaticamente quando estas estão identificadas no ficheiro correctamente, o que não é o caso no nosso mosaico virtual para QGIS. O comando que usámos acaba por criar a tal 4ª banda alfa mas não a identifica como tal. Para corrigir é fácil: basta abrir o ficheiro mosaico_final_qgis.vrt num editor de texto e adicionar o seguinte texto “<ColorInterp>Alpha</ColorInterp>”:

Conclusões

Com apenas alguns comandos conseguimos obter um mosaico sem áreas negras, usando para isso uma máscara em formato shapefile. O mosaico final tem uma excelente taxa de compressão, e sem ter de escolher uma cor para desaparecer da nossa imagem. E, ainda por cima, compatível com vários programas SIG… o que mais se pode desejar??

Por outro lado, este esquema oferece uma óptima flexibilidade: podemos editar a nossa máscara shapefile para também considerar as áreas negras dentro de ortos periféricos, resolvendo também este problema. E podemos reutilizar a máscara final, bastando copiá-la e renomeá-la…

Lista de Comandos Usados

E, como prometido, cá está a cábula:

1) Criar mosaico VRT a com imagens de uma directoria:
gdalbuildvrt mosaico_simples.vrt tif\*.tif
2) Criar mosaico virtual cortado pelo cutline e com banda alfa
gdalwarp -multi -wm 384 –config GDAL_CACHEMAX 384 -co alpha=yes -dstalpha -cutline cutlineteste.shp -of vrt mosaico_simples.vrt mosaico_cutline_alfa.vrt
3) Criar mosaico virtual rgb com maskline e sem banda alfa, para possibilitar máxima compressão jpeg
gdal_translate -b 1 -b 2 -b 3 -mask 4 -co alpha=no -co photometric=ycbcr -of vrt mosaico_cutline_alfa.vrt mosaico_mask.vrt
4) Criar mosaico tif final comprimido em jpeg qualidade 70
gdal_translate -co alpha=no -co photometric=ycbcr -co interleave=pixel -co tiled=yes -co compress=jpeg -co jpeg_quality=70 mosaico_mask.vrt mosaico_mask.tif
5) adicionar piramides
gdaladdo -ro -r average –config JPEG_QUALITY_OVERVIEW 70 –config COMPRESS_OVERVIEW JPEG –config INTERLEAVE_OVERVIEW PIXEL mosaico_mask.tif 2 4 8 16 32 64
6) Criar mosaico virtual para utilizar em QGIS com banda alfa
gdal_translate -of VRT -b 1 -b 2 -b 3 -b mask mosaico_mask.tif mosaico_mask_qgis.vrt

Voilá: um mosaico rápido, com áreas transparentes, que ocupa pouco espaço em disco, e que funciona em QGIS e ArcMap.