Tag Archives: Mosaicos

Raster Overviews

Overviews GDAL em modo Turbo

Tempo de leitura: 7 min

TLDR: Neste post discutimos formas de acelerar o processo de criação de overviews, e no fim usamos um script que reduz o tempo de processamento em 20%-50%. O script é apresentado abaixo e está no github.

Na visualização de rasters é obrigatório construir as overviews ou pirâmides, para conseguirmos uma visualização rápida.

As overviews são uma série de cópias do nosso raster com resoluções cada vez menores (pixeis maiores), e geralmente cada nível aplica uma redução de 50% na resolução. Por exemplo, numa resolução original de 0,30m/pixel, as overviews são imagens com resoluções de 0,60 – 1,20 – 2,40m/pixel e assim sucessivamente até que não faz sentido reduzir mais a imagem.

Overviews ou pirâmides permitem uma visualização rápida de rasters, através de imagens de resolução reduzida. (Obitdo em: https://eurogeographics.org/wp-content/uploads/2018/04/WCS-NLSS.pdf.)

Em geral, a construção destas pirâmides é feito com o comando gdaladdo, e é o processo mais moroso quando processamos grandes áreas. Nem a conversão com compressão, nem a união de muitos rasters leva tanto tempo.

Actualmente, com discos SSD rapidíssimos e memória super-abundante, e processadores multi-core, o comando gdaladdo que constrói overviews continua a usar apenas 1 core… por outro lado, é mais lento que outros comandos, como o gdal_translate.

Recentemente processei novos mosaicos para o Alentejo, desta vez com ortofotomapas com 0,30m de resolução, rgb+nir. E, claro, construir overviews foi uma tortura… mais de 8h para cada metade (dividi a área em 2 blocos este/oeste). O processador nunca passou dos 17% (i7 de 4 cores/8threads), e o disco SSD nunca passou de uns miseráveis 5MB/s (quando o disco é capaz de 1000MB/s). Muito frustrante…

O processo que uso consiste sempre em manter os ficheiros independentes, e criar um mosaico .vrt. Por hábito não crio mosaicos tif enormes. Este processo é descrito em artigos anteriores.

Depois de pesquisar online, vi 3 sugestões para melhorar o tempo de criar overviews:

Isto ensinou-me uma série de coisas novas:

  1. Os ficheiros .ovr são na verdade ficheiros TIFF multi-página (herança do tempo dos faxes!), onde um tiff é “colado” a outro dentro do mesmo ficheiro. Eu não sabia isto sobre os .ovr. Ou seja, cada resolução é um tiff, dentro do ovr, que é também um tiff (matrioska?).
  2. É possível juntar vários tiff num só tiff multi-página usando o comando “tiffcp tiff1 tiff2 tiffunido”.
  3. O OSGEO4W inclui uma versão “geo-activada” dos comandos tiff, mantendo as características SIG dos ficheiros.
  4. Podemos ter overviews de overviews, juntando a extensão .ovr ao ficheiro .ovr anterior, numa sucessão que funciona em gdal, qgis, e arcgis. Deve funcionar nos restantes programas, como geoserver, mapserver, etc.

Teste

Vamos fazer um teste com uma série mais pequena de ortofotomapas, para vermos qual é a melhoria no tempo de criação de overviews: vamos usar 3 processos simultânos de gdal_translate, onde cada processo constrói um resample diferente (x2, x4, x8), e renomeando-os para terem extensão .ovr acrescida.

A nossa coleção de ortos de testes é constituída por:

  • 29 ficheiros 3 bandas x 8 bit, num total de 276MB, já comprimidos em tiff/jpeg, com 5km de lado, e 0,30m de resolução.
  • Um mosaico virtual .vrt com todos os 29 ortos, “teste_script.vrt”, com dimensão de 166.667 x 100.000 pixeis.
Quadrículas dos 29 ortofotomapas do nosso teste.

O método consiste em executar 3 comandos em simultâneo:

  • Processo 1: resample para 0,6m/pixel
  • Processo 2: resample para 1,2m/pixel
  • Processo 3: resample para 2,4m/pixel, mais construção de pirâmides para este resample apenas

Assim, no processo 1 teremos este comando:

gdal_translate -of gtiff -tr 0.6 0.6 -ro -r average --config GDAL_CACHEMAX 1024 -co photometric=ycbcr -co interleave=pixel -co tiled=yes -co compress=jpeg teste_script.vrt teste_script.vrt.ovr

Ou seja, construimos uma cópia do mosaico, em formato tiff, com 2x o tamanho do pixel original (0,6m/pixel) e damos o nome certo para que seja automaticamente reconhecido como overviews do original -> teste_script.vrt.ovr.

No processo 2, construímos um resample com 4x o tamanho do pixel (1,2m/pixel), e damos o nome que o faz ser reconhecido como overviews do 1º processo:

gdal_translate -of gtiff -tr 1.2 1.2 -ro -r average --config GDAL_CACHEMAX 1024 -co photometric=ycbcr -co interleave=pixel -co tiled=yes -co compress=jpeg teste_script.vrt teste_script.vrt.ovr.ovr  

No processo 3, construímos o 3º nível, com 8x a resolução (2,4m/pixel), e com um nome que o marque como as overviews do 2º nível:

gdal_translate -of gtiff -tr 2.4 2.4 -ro -r average --config GDAL_CACHEMAX 1024 -co photometric=ycbcr -co interleave=pixel -co tiled=yes -co compress=jpeg teste_script.vrt teste_script.vrt.ovr.ovr.ovr

Já que sabemos que este resample é muitíssimo mais rápido que o 1º, terminando por isso muito cedo, podemos aproveitar para criar pirâmides para este 3º nível. Isto permitirá termos a série completa de overviews no final:

gdaladdo -ro -r average --config GDAL_CACHEMAX 1024 --config COMPRESS_OVERVIEW JPEG --config PHOTOMETRIC_OVERVIEW YCBCR --config INTERLEAVE_OVERVIEW PIXEL  teste_script.vrt.ovr.ovr.ovr

O script faz uma série de correções aos nomes dos ficheiros caso detecte uma máscara externa, que é o caso do nosso teste (ficheiro .msk). Ficamos assim com os seguintes ficheiros:

      50 194 teste_script.vrt
  12 434 020 teste_script.vrt.msk
   6 003 920 teste_script.vrt.msk.ovr
   1 525 126 teste_script.vrt.msk.ovr.ovr
     395 087 teste_script.vrt.msk.ovr.ovr.ovr
     267 858 teste_script.vrt.msk.ovr.ovr.ovr.ovr
 466 889 237 teste_script.vrt.ovr
 116 024 726 teste_script.vrt.ovr.ovr
  32 529 152 teste_script.vrt.ovr.ovr.ovr
  11 135 760 teste_script.vrt.ovr.ovr.ovr.ovr
               10 File(s)    647 255 080 bytes

O tempo de execução foi de 06:47,4 min. E podemos ver a ocupação de CPU, disco e memória muito mais altos:

E funciona? Vamos a ver…

gdalinfo teste_script.vrt
 Driver: VRT/Virtual Raster
 Files: teste_script.vrt
        teste_script.vrt.ovr
        teste_script.vrt.ovr.ovr
        teste_script.vrt.ovr.ovr.ovr
        teste_script.vrt.ovr.ovr.ovr.ovr
        teste_script.vrt.msk
        230_060_irg.tif
        230_065_irg.tif
...  ...  ...  ...  ...
Size is 166667, 100000
 Coordinate System is:
... ... ... ... ...
Band 1 Block=128x128 Type=Byte, ColorInterp=Red
   Min=0.000 Max=255.000
   Minimum=0.000, Maximum=255.000, Mean=109.688, StdDev=100.399
   Overviews: 83334x50000, 41667x25000, 20833x12500, 10417x6250, 5209x3125, 2605x1563, 1303x782, 652x391, 326x196, 163x98
   Mask Flags: PER_DATASET
 ...  ...  ...  ...  ... 

O gdalinfo reconhe todas as pirâmides. E o QGIS?

Identificação das overviews pelo QGIS.

Pequeno à parte: Já em artigos anteriores referi que o QGIS tem de ser “convencido” a ler máscaras externas. Isto não causa problemas ao processo. Aparentemente, o GDAL tem um comportamento diferente com máscaras externas, em que as expõe com valores 0/1, em vez de valores 0/255 como acontece com máscaras internas. Sem esta correção a máscara não é detectadas correctamente pelo QGIS, e temos de a ignorar, aparecendo as zonas pretas sem dados. Se corrigirmos editando o vrt, tudo aparece correctamente. Mas em qualquer dos casos as overviews funcionam:

QGIS e overviews em ação, velocidade real.

Com o gdaladdo “normal”

Para compararmos, vamos criar overviews com o processo normal:

timing "gdaladdo --config GDAL_CACHEMAX 1024 --config COMPRESS_OVERVIEW JPEG --config PHOTOMETRIC_OVERVIEW YCBCR --config INTERLEAVE_OVERVIEW PIXEL teste_gdaladdo.vrt"
14:52:45,10
a executar o comando indicado: "gdaladdo --config GDAL_CACHEMAX 1024 --config COMPRESS_OVERVIEW JPEG --config PHOTOMETRIC_OVERVIEW YCBCR --config INTERLEAVE_OVERVIEW PIXEL teste_gdaladdo.vrt"
0…10…20…30…40…50…60…70…80…90…100 - done.
0…10…20…30…40…50…60…70…80…90…100 - done.
15:01:17,19

Assim, o processo normal demorou 08:32,1 min.

E os ficheiros deste processo normal são:

      50 194 teste_gdaladdo.vrt
  23 195 175 teste_gdaladdo.vrt.msk
 827 911 591 teste_gdaladdo.vrt.ovr
  3 File(s)    851 156 960 bytes

Nota: este vrt tem uma máscara externa .msk. Como me esqueci do parâmetro -ro (readonly), as overviews da máscara foram adicionadas ao próprio msk. Também me esqueci do método resample average, que seria mais lento…

Comparação

O processo de construir overviews em paralelo, divido em 3 processos simultâneos, é 20% mais rápido, e ainda com um bónus de ocupar menos 23% em disco! (não sei porquê)

MétodoTempoTamanhoDiscoCPURAM
gdaladdo normal 08:32,1 850MB5MB/s15%1GB
gdal_translate x3 06:47,4 647 MB19MB/s45%3GB

Ou seja, conseguimos subir o uso do CPU para 45%, e o disco para 19MB/s. Nada mau. A memória ocupada pelo processo depende do uso que fizermos da flag –config GDAL_CACHEMAX. No nosso caso, definimo-la como 1GB. Logo 3 processos ocupam obviamente 3x esta quantidade.

O ganho de velocidade resulta do processamento em simultâneo – enquanto se processa o 1º nível, processam-se logo os restantes e as máscaras também desses níveis caso exista máscara no raster original.

Script

Numa tentativa de automatizar o processo, criei um script bat para windows. Pode ser obtido aqui: https://github.com/dncpax/Turbo_GDAL_Overviews .

Algumas notas interessantes sobre o bat:

  • É possível imitar uma execução em background usando “start /b” dentro do bat.
  • A shell DOS só faz aritmética de inteiros, por isso temos de indicar as 3 resoluções que queremos – x2, x4 e x8 – porque em geral não são inteiros e não os conseguimos calcular no bat.
  • Usamos DELAYEDEXPANSION porque precisamos de mostrar o tempo de execução.
  • Temos de renomear os ficheiros resultantes da máscara externa (.msk) porque ficam com nomes que impedem o seu reconhecimento. O script trata disso com uma série de renames.

A migração para Linux deve ser fácil, porque 80% do script é só validação de argumentos. O que interessa são comandos gdal. Voluntários procuram-se…

Para executar indicamos o raster e as 3 resoluções iniciais das overviews:

turbo_overviews.bat teste_script.vrt 0.6 1.2 2.4  
Inicio em: 13:29:24,69  
Input file size is 166667, 100000  
Input file size is 166667, 100000  00
Input file size is 166667, 100000  
0…10…20…30…40…50….60…70…..80…..9010….100 - done.  
..0…10…20…30….40…50…60…70…80…90.20..100 - done.  
0…10…20…30…40…50…60…70…80…90….100 - done.  
10..30…..40…20..50…..6030…..70…40.80…..90…50.100 - done.  
…60…70…80…90…100 - done.  
Fim em: 13:36:12,07

Conclusões

Se calhar este post é optimista: só fiz 1 teste sério… pode ser o caso de não funcionar mais vez nenhuma ;)…

Para mim é realmente estranho a falta de processamento multi-core no GDAL. Talvez seja uma questão de tempo, mas já se sente a falta. O que existe é muito incipiente e apenas funciona em cenários que não me são aplicáveis (e.g. compressão DEFLATE).

Há mais alternativas a este processo, mas nos testes que fiz não tive tão bons resultados.

Há outros programas que podem fazer resample de imagens e com processamento multi-core, como o imagemagick. Isto obrigará a copiar a georeferenciação para as imagens resultantes porque estes programas não reconhecem a componente SIG das imagens. Mas pode ser interessante.

De qualquer forma, por agora, este processo parece funcionar bem. Falta testar com um mosaico “à séria”. Pode ser que os 20% de maior rapidez se confirmem!

Até breve!

Adenda

Teste com um mosaico um pouco maior…

Este mosaico é similar mas maior: tem uma dimensão de 233.334 x 383.334 pixeis, em 258 ficheiros, num total de 10,3GB. Demorou 41:38,2 min, em vez de 01:35 h com o GDAL v.3.0.4 (na v2.3.0 tinha demorado 08:42:33 h) do gdaladdo… Ou seja, um ganho de 56%!

Observámos uma ocupação do disco interessante: mais de 60MB/s…

Mosaico maior em construção…

Os tamanhos foram similares: 4.4GB vs 4.9GB.

E correctamente visualizado em QGIS:

Mosaico de 10GB com overviews criadas em 40min…

Neste mosaico mais encorpado, a melhoria é enorme… Curioso em ver a aplicação em mosaicos muito maiores.

Mosaicos com áreas transparentes

Tempo de leitura: 9 minNeste 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 😉 Ler artigo completo

Mosaicos de imagens em MapServer, com GDAL

Tempo de leitura: 12 min

Quando fica incumbido de publicar informação geográfica na web é quase sempre confrontado com esta questão – como publicar aquela colecção de ortos, composta por uma directoria cheia de ficheiros tif? (ou jpeg, ou ecw, ou sid…)

Este artigo nasceu num desses momentos, e da necessidade de determinar quais as opções existentes no MapServer, e saber qual delas é a melhor.

AVISO – artigo longo, com resultados de uma série de testes!!! Conclusões no final para os mais stressados…

Cenário inicial

Imaginemos um pequeno concelho que tem uma colecção de 45 ortofotomapas, em formato TIFF, não comprimido. Cada imagem ocupa 234MB, o que totaliza 10GB de imagens.

O objectivo é publicar estas imagens como uma cobertura que abrange o concelho, numa aplicação webgis, servida pelo MapServer.

A melhor opção

Para escolher a melhor forma de publicar as imagens é preciso saber o que é mais importante para nós: velocidade ou espaço em disco?

A verdade é que a compressão das imagens impõe uma penalização na velocidade com que os programas conseguem mostrar essas imagens. Essa penalização pode ser muito pequena ou muito grande, consoante o tipo de compressão e o tipo de imagem criada durante o processo de compressão. Geralmente, os ficheiros ECW e SID são extremamente comprimidos e relativamente rápidos. Em geral também, o formato JPEG2000 comprime muito as imagens mas é lento. O formato mais antigo JPEG é um compromisso entre os 2 grupos anteriores. Outra regra: assume-se que o MapServer é mais rápido com o formato TIFF, sem compressão, do que com formatos comprimidos, mas estas imagens ocupam muito mais espaço em disco.

Portanto, se a sua única preocupação é velocidade, e tem espaço em disco que não se importa de gastar com os seus ortofotomapas, em principio já sabe a sua resposta – use TIFFs sem compressão. Mas nada como verificar se esta “verdade” realmente se aplica aos seus dados em particular. Ahh, e não se esqueça de considerar o tráfego que vai gerar na rede se quiser usar os ortos em programas SIG desktop…

Temos ainda de saber como criar um layer em MapServer que use todas as imagens como um conjunto único. Afinal não quer que o utilizador seja obrigado a ligar 45 imagens uma-por-uma.

Para iniciar os testes converteram-se todas as imagens para os diversos formatos. Os tamanhos totais em cada formato ficaram assim:

Formato Dimensão
TIF+overheads 13,5 GB
JPG+overheads 915 MB
ECW 605 MB
JP2 509 MB

Os tamanhos indicados incluem imagens de resolução reduzida, chamadas overheads ou pirâmides. Mais sobre isto adiante…

Para ver os comandos de conversão do GDAL para cada formato pode consultar este artigo anterior – GDAL: Formatos comprimidos.

Usar mosaicos em MapServer

Um mosaico de imagens em MapServer é um layer do tipo RASTER que aponta para um conjunto de imagens em vez de uma só.

A forma tradicional de referenciar um mosaico é criando um shapefile que contém um quadriculado com as extensões das imagens. Este shapefile é criado com um comando do GDAL chamado gdaltindex.

Mas há uma nova opção. Podemos usar o formato VRT, um formato virtual que é constituído por um ficheiro de texto que indica a fonte dos dados e a forma como são organizados. Um ficheiro VRT pode listar um conjunto de imagens de forma a que o MapServer e outras aplicações o leiam como uma imagem única, quando na verdade será composto por todas as nossas imagens. O comando para criar um mosaico VRT é o gdalbuildvrt.

Por último, temos a opção de força bruta: juntar todas as imagens numa só, criando uma imagem enorme em disco que abrange toda a nossa área de interesse. Ou seja, no nosso caso em estudo esta imagem teria 10GB. Para criar este mosaico monolítico usamos o comando gdal_merge.

Vamos ao longo do artigo analisar cada opção e no final medir os tempos que o MapServer demora a visualizar cada tipo de mosaico.

Mosaico VRT

O GDAL suporta um formato virtual, VRT, que lista ficheiros já existentes e os descreve como um todo. Podemos agregar várias imagens numa só, definir um novo sistema de coordenadas, ou no caso de ficheiros vectoriais definir filtros de atributos, e renomear atributos, tudo sem alterar os ficheiros originais. Os ficheiros VRT são ficheiros de texto em formato XML e podem ser editados manualmente, ou criados com o gdal_translate ou com o comando gdalbuildvrt. Mais info aqui.

Para construir um mosaico a partir das imagens numa directoria, basta usar o seguinte comando:

gdalbuildvrt mosaico_tif.vrt directoria\*.tif

O ficheiro “mosaico_tif.vrt” é criado com o seguinte conteúdo:

<VRTDataset rasterXSize="40000" rasterYSize="20000">
  <SRS>LOCAL_CS[&quot;unnamed&quot;,UNIT[&quot;metre&quot;,1,AUTHORITY[&quot;EPSG&quot;,&quot;9001&quot;]]]</SRS>
  <GeoTransform>-1.6000000000000000e+004, 5.0000000000000000e-001, 0.0000000000000000e+000,-7.0000000000000000e+004, 0.0000000000000000e+000,-5.0000000000000000e-001</GeoTransform>
  <VRTRasterBand dataType="Byte" band="1">
    <ColorInterp>Red</ColorInterp>
    <SimpleSource>
      <SourceFilename relativeToVRT="1">tif\003941Argbx.tif</SourceFilename>
      <SourceBand>1</SourceBand>
      <SourceProperties RasterXSize="8000" RasterYSize="10000" DataType="Byte" BlockXSize="8000" BlockYSize="1"/>
      <SrcRect xOff="0" yOff="0" xSize="8000" ySize="10000"/>
      <DstRect xOff="0" yOff="0" xSize="8000" ySize="10000"/>
    </SimpleSource>

… seguem-se as restantes imagens em sucessivos SimpleSource até terminar a lista de imagens, para a 1ª banda. Depois inicia-se a 2ª banda:

</VRTRasterBand>
  <VRTRasterBand dataType="Byte" band="2">
    <ColorInterp>Green</ColorInterp>
    <SimpleSource>
      <SourceFilename relativeToVRT="1">tif\003941Argbx.tif</SourceFilename>
      <SourceBand>2</SourceBand>
      <SourceProperties RasterXSize="8000" RasterYSize="10000" DataType="Byte" BlockXSize="8000" BlockYSize="1"/>
      <SrcRect xOff="0" yOff="0" xSize="8000" ySize="10000"/>
      <DstRect xOff="0" yOff="0" xSize="8000" ySize="10000"/>
    </SimpleSource>

… seguem depois as imagens para compor a 3ª banda…

</VRTRasterBand>
  <VRTRasterBand dataType="Byte" band="3">
    <ColorInterp>Blue</ColorInterp>
    <SimpleSource>
      <SourceFilename relativeToVRT="1">tif\003941Argbx.tif</SourceFilename>
      <SourceBand>3</SourceBand>
      <SourceProperties RasterXSize="8000" RasterYSize="10000" DataType="Byte" BlockXSize="8000" BlockYSize="1"/>
      <SrcRect xOff="0" yOff="0" xSize="8000" ySize="10000"/>
      <DstRect xOff="0" yOff="0" xSize="8000" ySize="10000"/>
    </SimpleSource>

Finalizando-se o ficheiro assim:

   </VRTRasterBand>
</VRTDataset>

O QGIS 1.4.0 consegue ler este tipo de ficheiro, mostrando-o como uma só imagem:

Mosaico VRT visualizado no QGIS 1.4.0
Mosaico VRT visualizado no Quantum GIS 1.4.0.

Mas a performance é um problema. Para visualizar este mosaico, o QGIS tem de abrir as 45 imagens, ler todos os pixels, e mostrá-los. Ainda por cima, nesta escala a maior parte da informação é inútil porque não se distinguem todos os pixels. Para resolver este problema, criam-se overheads…

Criar Overheads

Os ficheiro TIFF originais têm 234MB cada um, e para visualizar todos os ortos o QGIS tem de ler todas as imagens, num total de 9GB, e mostrá-las numa escala onde o todo pormenor se perde. Este problema aplica-se a todos os programas SIG, incluindo o MapServer.

Para evitar este desperdício de tempo, criam-se overheads ou pirâmides, que são imagens de resolução decrescente gravadas num só ficheiro com extensão OVR. Assim, para uma escala pequena, como a do exemplo anterior, o QGIS tem apenas de ler e mostrar a imagem de pior resolução, já que se ajusta bem ao pormenor visível a essa escala. Consoante o utilizador faz “zoom in”, o QGIS selecciona a imagem com a resolução apropriada a essa escala, até que a partir de uma dada escala começa a mostrar a imagem original. Mas neste momento já só é necessário ler uma pequena porção da imagem para cobrir a área visível. E assim se acelera a visualização das imagens.

Para criar overheads usamos o comando do GDAL, gdaladdo:

gdaladdo -ro <ficheiro_de_imagem> 2 4 8 16 32 64 128

O parâmetro “-ro” indica que a imagem original não deve ser alterada, e portanto as overheads serão criadas num ficheiro separado, que terá a extensão OVR.

Os números depois do nome da imagem são os níveis a criar com resolução reduzida. O n.º 2 indica metade da resolução original, 4 é 1/4 e assim sucessivamente. Até que nível de redução calculamos depende da extensão geográfica das imagens e da sua resolução original. Uma forma fácil de determinar até que redução devemos chegar é consultar o QGIS…

O QGIS consegue agora criar overheads – acedendo às propriedades de um tema de imagem, na secção de Pyramids existe um botão para criar pirâmides (embora seja mais lento que o comando GDAL). Aqui o QGIS mostra a lista de resoluções que aconselha – é só contar quantos níveis são aconselhados pelo QGIS. Na imagem seguinte, podemos ver o QGIS a criar pirâmides para uma das imagens:

Construir pirâmides/overheads no Quantum GIS 1.4.0.
Construir pirâmides/overheads no Quantum GIS 1.4.0.

Como o QGIS indicava 7 níveis de pirâmides para as minhas imagens, foi o que usei no comando: 2, 4, 8, 16, 32, 64, 128 = 7 níveis.

Criei então overheads para o mosaico VRT. O resultado foi um ficheiro OVR com 776MB, ou seja, 33% do tamanho original. Esta é a taxa comum na criação de overheads e temos de prever este consumo adicional de disco.

Para testar as overheads, carreguei de novo o mosaico VRT no QGIS. Agora a imagem total foi mostrada quase instantaneamente. Mas ao fazer zoom in a velocidade degradou-se. Significa que o QGIS deixa de usar as pirâmides do VRT e passa a mostrar as imagens originais… penso que poderá ser um bug.

Portanto a solução passa por criar overheads para cada uma das imagens originais, e ver se assim o QGIS utiliza sempre as overheads.

Para isso usei um comando de linha que percorre todas as imagens TIFF numa directoria e executa o comando gdaladdo:

for %I in (directoria_originais\*.tif) do gdaladdo -ro %I 2 4 8 16 32 64 128

O resultado foi a criação de um ficheiro OVR para cada TIFF, cada um com 82MB (35% dos originais). Voltando a carregar o mosaico VRT no QGIS, o desempenho foi excelente em todos os níveis de zoom, provando que o QGIS usou sempre as overheads.

NOTA: por curiosidade abri o mosaico VRT no ArcGIS, e tudo funcionou perfeitamente, incluindo as overheads! Sabe-se há muito que o ArcGIS incluía o GDAL, não se sabendo bem para o quê. Pelos vistos, servirá também para ler VRT’s e overheads.

O caso ECW

O formato ECW agrada-me muito. Tem compressões equivalentes ao MrSID, é gratuito para comprimir imagens até 500MB, é rápido a visualizar, e já inclui pirâmides no próprio ficheiro. Poupa assim espaço em disco 2x: na compressão e nas pirâmides.

A questão que se coloca é saber até que ponto vai a penalização na performance do MapServer ao ter que descomprimir os ECW para os visualizar a cada Zoom e a cada Pan.

Para criar um mosaico VRT com ECW’s tivemos de converter todos os TIFF originais, com o seguinte comando:

for %I in (directoria_originais\*.tif) do gdal_translate -of ECW -co TARGET=90 %I directoria_destino\%~nI.ecw

As imagens originais passaram de 13GB (contando com as overheads) para 605MB – poupando 95% do espaço! Podiamos ainda comprimir mais, mas neste estudo preferi manter os defaults…

Criámos um novo VRT com as imagens ECW, e testámos no QGIS. A visualização foi muito rápida e sem degradação ao aumentar a escala.

O caso JPEG

Este formato é muito apreciado por ser familiar, bastante rápido, e apresentar taxas de compressão bastante atractivas (embora não consiga acompanhar a compressão dos novos formatos como MrSID e ECW). A questão que se coloca é saber se a multa devida pela descompressão se nota ou não em MapServer. Não esquecer que será necessário criar overheads. Neste teste optei por criar as pirâmides também em formato JPEG, por forma a ser mais fiel ao formato. Se este formato for suficientemente rápido poderá ser uma boa alternativa ao TIFF. Para criar pirâmides em formato JPEG pode-se usar o seguinte comando (está tudo documentado na página do gdaladdo):

gdaladdo --config COMPRESS_OVERVIEW JPEG --config PHOTOMETRIC_OVERVIEW YCBCR --config INTERLEAVE_OVERVIEW PIXEL imagem 2 4 8 16 32 64 128

O resultado foram imagens JPEG e OVR ocupando 915MB (93% de compressão). Em QGIS o resultado foi mais uma vez excelente… visualização rápida de inicio, sem degradação com os zoom in… Além da velocidade, fiquei surpreendido com a compressão atingida e com a óptima qualidade das imagens ao ver todo o mosaico.

O caso JPEG2000

Este formato é uma alternativa aos formatos ECW e MrSID, mas mais aberto. No entanto, tem-se revelado sistematicamente lento demais para o poder utilizar a não ser para fins de arquivo.

Decidi mesmo assim inclui-lo no teste. Converti todas as imagens TIFF para JP2, e criei um mosaico VRT. No teste com QGIS, a velocidade foi muito rápida de início, apenas demorando um pouco mais com o aumento da escala de visualização. No entanto, a sua qualidade de imagem ao ver todo o mosaico foi excelente.

Mosaicos Shapefile (TILEINDEX)

A forma tradicional de usar mosaicos no MapServer é criar um shapefile com polígonos que cobrem a área de cada ortofotomapa. Este shapefile é depois referenciado no mapfile em vez do ficheiro VRT.

A questão que se coloca é saber se a utilização de um VRT implica uma redução de performance, quando comparada com a utilização de um shapefile.

Criei por isso mosaicos usando shapefiles, e medi os tempos de visualização em MapServer. Afinal, nem sempre novos métodos se revelam melhores que os antigos. O comando que usei foi o seguinte:

gdaltindex mosaico.shp directoria\*.tif

Não descobri forma de visualizar este tipo de mosaico em QGIS…

Mosaicos monolíticos

A última variante de mosaicos conhecida pela humanidade – juntar todas as imagens numa única, e gigante, imagem.

Para criar este mosaico vamos usar a ferramenta gdal_merge, que é um script python incluído no GDAL (pelo menos na versão FWTools). O plano é juntar todos os tiffs originais, criando um novo tiff. Mas nesta altura, o espaço em disco começava a escassear…

Para poupar algum espaço em disco, optei por usar compressão JPEG. Como o ficheiro poderá ser maior que 4GB usei a opção “bigtiff=yes” (por limitações do formato tiff). E como o ficheiro cobrirá uma grande extensão geográfica, para acelerar o acesso a pequenas áreas, usei também a opção “tiled=yes”. Assim, o comando para criar o mosaico foi o seguinte:

gdal_merge -o mosaico.tif -of gtiff -co compress=jpeg -co tiled=yes -co tfw=yes -co bigtiff=yes -v imagens_originais\*.tif

Isto produziu um ficheiro com 3,5GB, dada a compressão JPEG dentro do ficheiro TIFF. Confuso? Sucede que o formato TIFF suporta uma variedade de processos de compressão. Um deles é o JPEG. Claro que esta compressão poderá comprometer a velocidade do mosaico, mas não tinha mais espaço em disco para criar um mosaico sem compressão.

O passo seguinte foi criar overheads/pirâmides para este mosaico. Também aqui optei por criar pirâmides comprimidas em JPEG, usando o seguinte comando:

gdaladdo --config COMPRESS_OVERVIEW JPEG --config PHOTOMETRIC_OVERVIEW YCBCR --config INTERLEAVE_OVERVIEW PIXEL -ro mosaico.tif 2 4 8 16 32 64 128

Este comando criou um ficheiro OVR com 320MB, sendo apenas 9% do original devido à compressão usada. No total, a compressão foi de 72%. Nada mau…

Em QGIS este mosaico foi extremamente rápido, mostrando o ficheiro numa fracção de segundo, e zooms e pans foram também muito rápidos. Definitivamente, foi das melhores prestações senão a melhor.

NOTA: seria interessante criar um mosaico enorme ECW, mas não foi possível por causa do limite de 500MB imposto pela licença gratuita do compressor ECW (incluído no GDAL).

E agora… com MapServer

Todos os testes foram efectuados num portátil,  com Windows 7, 32-bit, 2 GB de memória, e MS4W 2.2.7 (inclui o MapServer 5.0.2). Portanto, podemos esperar melhores resultados num PC ou servidor, onde o disco rígido será em princípio bastante mais rápido. Os testes foram feitos com uma pequena aplicação OpenLayers, que inicia com um mosaico visível em toda a extensão, sendo depois feitos 4 zooms no centro do mapa, manualmente.

Os tempos foram obtidos no log produzido pelo MapServer. Para criar este log, inseriram-se as seguintes linhas no mapfile:

MAP
    CONFIG  "MS_ERRORFILE" "/ms4w/mapserver.log"
    DEBUG 5

O mapa foi configurado para gerar imagens em JPEG, através do AGG, usando esta secção no mapfile:

  OUTPUTFORMAT
    NAME jpg
    DRIVER AGG/JPEG
    IMAGEMODE RGB
  END

Para definir um layer com um mosaico VRT usou-se a seguinte sintaxe no mapfile, indicando o ficheiro VRT no tag DATA:

   LAYER
    NAME 'mosaicoTotal_tif'
    TYPE RASTER
    DUMP true
    TEMPLATE fooOnlyForWMSGetFeatureInfo
    EXTENT -16500.000000 -83122.641509 4500.000000 -66877.358491
    DATA 'mosaicoTotal_tifoverhds.vrt'
    METADATA
      'ows_title' 'mosaicoTotal_tif'
    END
    STATUS OFF
    TRANSPARENCY 100
    PROJECTION
    'init=EPSG:27492'
    END
  END

Note-se que foram criados vários layers deste tipo, havendo um VRT para cada tipo de imagem: TIFF, ECW, JPG, JP2.

Para os mosaicos baseados em shapefile usou-se a seguinte sintaxe (TILEINDEX e TILEITEM no lugar de DATA):

   LAYER     NAME 'mosaicoTotalshp_tif'     TYPE RASTER     DUMP true     TEMPLATE fooOnlyForWMSGetFeatureInfo     EXTENT -16500.000000 -83122.641509 4500.000000 -66877.358491     TILEINDEX 'mosaicoshp_tif.shp'     TILEITEM 'location'     METADATA       'ows_title' 'mosaicoTotalshp_tif'     END     STATUS OFF     TRANSPARENCY 100     PROJECTION     'init=EPSG:27492'     END   END

Resultados do MapServer

Usando mosaicos VRT obtiveram-se os seguintes tempos:

Zoom TIFF+overhds ECW JP2 JPG TIFmonolítico
Inicial 0.298 0.744 1.108 0.81 0.094
Z1 0.354 0.609 3.205 0.945 0.259
Z2 0.419 0.761 3.173 1.018 0.291
Z3 0.415 0.709 3.053 1 0.284
Z4 0.65 0.658 2.533 0.966 0.276
Totais 2.136 3.481 13.072 4.739 1.204
Espaço disco 13,5GB 605MB 509MB 915MB 3,85GB*

*usando compressão JPEG, e incluindo overheads.

E em gráfico (a última categoria “Totais” representa a soma de todos os zooms, evidenciando os ganhos acumulados nos formatos mais rápidos):

Desempenho de mosaicos VRT em MapServer.
Desempenho de mosaicos VRT em MapServer.

E o vencedor é… o mosaico monolítico, destacadíssimo. Dos mosaicos VRT, o mais rápido foi o mosaico de TIFF’s não comprimidos, como se supunha de início, com o formato ECW em 2º lugar, bastante melhor que o JPEG em 3º. O formato JPEG2000 é completamente desaconselhado.

Usando mosaicos baseados em Shapefile, obtiveram-se os seguintes resultados:

Zoom TIFF+overhds ECW JP2 JPG TIFmonolítico
Inicial 0.287 2.374 4.206 2.366 0.094
Z1 0.337 2.355 6.121 2.355 0.259
Z2 0.247 1.082 3.583 1.139 0.291
Z3 0.186 0.656 3.069 0.76 0.284
Z4 0.176 0.657 2.412 0.926 0.276
Totais 1.233 7.124 19.391 7.546 1.204
Espaço disco 13,5GB 605MB 509MB 915MB 3,85GB*

*usando compressão JPEG, e incluindo overheads.

E em gráfico:

Desempenho de mosaicos SHP em MapServer Desempenho de mosaicos SHP em MapServer

 

E o vencedor é… de novo o TIFF monolítico, mas havendo várias alterações. Só o mosaico de TIFF’s não comprimidos melhorou o desempenho em quase 100%, reduzindo de 2,1s para 1,2s. Todos os outros pioraram bastante, por vezes duplicando o tempo de visualização. Não encontrei explicação para esta penalização. O shapefile continha apenas 45 polígonos, sendo difícil acreditar que possa causar tão grande penalização, mesmo considerando que o MapServer tem de determinar quais os polígonos visíveis num momento, determinar as áreas de cada um, e abrir as imagens correspondentes… parece-me pouco trabalho para explicar diferenças de 1,6s como no caso dos ECW… Mas factos são factos.

Conclusões

Se não existirem impedimentos, deverá unir as suas imagens num mosaico único, em formato TIFF, comprimido em JPEG, e criar pirâmides também com compressão JPEG. Obtém o maior desempenho possível e com o extra de poupar 70% do espaço em disco. Esta abordagem foi simplesmente a melhor… mais vale guardar os seus originais num arquivo.

Quanto aos mosaicos VRT mostraram ser uma boa opção. Quase todos os formatos tiveram melhor desempenho quando usados através deste formato virtual, do que usando um mosaico baseado em shapefile. A excepção foi o mosaico de TIFF’s não comprimidos. Aqui, usar o mosaico shapefile foi melhor, muito melhor.

Nota Final

O formato VRT não serve apenas para criar mosaicos sem alterar as imagens originais, e por isso poderá ser uma opção atractiva em certas situações. Podemos, por exemplo, definir uma deslocação de +200km, +300km sem andar a editar ficheiros de coordenadas – alguém se lembra de passar por isto?? Parece-me que sim.

Cábula de Comandos

Todos os comandos usados no artigo…

Criar mosaico VRT a com imagens de uma directoria:
gdalbuildvrt mosaico_tif.vrt directoria\*.tif
Criar overheads/pirâmides de uma imagem, num ficheiro separado
gdaladdo -ro imagem.tif 2 4 8 16 32 64 128
Criar overheads/pirâmides para todas as imagens TIFF numa directoria
for %I in (directoria_originais\*.tif) do gdaladdo -ro %I 2 4 8 16 32 64 128
Converter todos os TIFF numa directoria para ECW noutra directoria
for %I in (directoria_originais\*.tif) do gdal_translate -of ECW -co TARGET=90 %I directoria_destino\%~nI.ecw
Criar overheads/pirâmides com compressão JPEG
gdaladdo –config COMPRESS_OVERVIEW JPEG –config PHOTOMETRIC_OVERVIEW YCBCR –config INTERLEAVE_OVERVIEW PIXEL imagem.tif 2 4 8 16 32 64 128
Criar mosaico baseado em shapefile das imagens TIFF numa directoria
gdaltindex mosaico.shp directoria\*.tif
Criar um mosaico monolítico tif, com compressão jpeg, das imagens tiff numa directoria
gdal_merge -o mosaico.tif -of gtiff -co compress=jpeg -co tiled=yes -co tfw=yes -co bigtiff=yes -v imagens_originais\*.tif
MapServer: configurar um mapfile para gerar um ficheiro log com tempos
MAP    CONFIG  “MS_ERRORFILE” “/ms4w/mapserver.log”    DEBUG 5