Tempo de leitura: 9 min

Introdução

Há alguns anos publiquei um artigo sobre a conversão de um modelo digital do terreno global para o sistema de coordenadas português e também cortado aos limites de Portugal Continental (http://192.168.1.84:8000/2010/03/mdt-30m-para-portugal/). Depois coloquei-o numa partilha online. Até hoje continua a ser procurado, embora tenha já feito várias referências para novos dados – melhores e mais atuais.

Os dados que recomendo são da Agência Europeia do Ambiente – o EU-DEM – que pode ser obtido aqui: https://land.copernicus.eu/imagery-in-situ/eu-dem/eu-dem-v1.1?tab=metadata. Estes dados são baseados no SRTM, melhorados com uma série de correções que estão documentadas.

Este mdt tem 25m de resolução, com um erro médio quadrático de +-7m!! O que me parece excelente.

Nota técnica – na realidade estes dados não são um mdt, mas sim um mds – modelo digital da superfície, ou seja, não representam a cota do terreno e sim o topo dos objetos na superfície, como árvores e edifícios e outras estruturas. Mas mesmo isto não é bem correto no caso do eu-dem… Aparentemente, estes dados podem representar o topo de objetos, como árvores e estruturas, mas, por outro lado, como sinais radar (usados na missão original) podem penetrar a canópia das árvores, também não há certeza que os dados representem a cota do topo das árvores.

Este artigo é muito simples – pretende disponibilizar uma versão pronta a usar para Portugal. Nada de especial – vamos apenas derivar uma versão que estará projetada para o nosso sistema de coordenadas, e cortada à extensão do nosso país (continente), usando apenas o QGIS. No fim do artigo há um link para descarregar os dados finais preparados para Portugal.

Obter os dados

Os dados são fáceis de obter – estão divididos em quadrículas de 1.000km de lado, em formato tiff, 32bit, e Portugal abrange 2 destas quadrículas. Podemos vê-las usando o visualizador que encontramos no link acima:

Como vemos, temos de descarregar as quadrículas E10N10 e E20N20. Vamos ao separador Download, e escolhemos estas quadrículas. Descarregamos e passamos ao processamento em QGIS.

Visualizar em QGIS

Os dados estão zipados. Descomprimimos e obtemos 1 tif por quadrícula que podemos carregar no QGIS:

Abrimos as propriedades dos 2 ficheiros e vemos a seguinte informação sobre os ficheiros:

• Projeção – IGNF:ETRS89LAEA – ETRS89 Lambert Azimutal Equal Area
• Resolução espacial: 25m
• Largura e Altura: 40.000×40.000 pixeis (1.000kmx1.000km)
• Tipo de pixel: 32bit
• Compressão: LZW
• Pirâmides: Sim
• Valor “nulo” (no data): -3.40282e+38
• Tamanho total em disco: 983.45MB + 956.04 MB

Plano de Ação

Então qual é o nosso plano de ação?

1. Vamos juntar os 2 ficheiros num mosaico virtual, que praticamente não ocupa espaço em disco, mas representa todos os dados como um só conjunto;
2. Vamos cortar este mosaico virtual pelos limites de portugal;
3. E vamos reprojetar para o sistema de coordenadas português. Vamos também escolher uma compressão eficiente e que conserve os dados sem alterações;
4. E no fim, zipamos e colocamos algures online para a malta usar.

Para o corte vamos usar a Carta Administrativa Oficial Portuguese disponível no site da DGT.

1. Unir 2 ficheiros raster num mosaico virtual

Isto é um truque – criar ficheiros virtuais que são apenas pequenos ficheiros de texto que descrevem as operações que queremos fazer, mas sem as fazer realmente. Isto permite criar novas representações dos nossos dados sem ocupar espaço em disco, mas em troca, ocupamos mais o processador e a memória. Já falei várias vezes sobre este fantástico formato VRT, e pode-se ler mais sobre ficheiros VRT aqui.

Para criar um ficheiro VRT a partir de outros ficheiros, usamos a ferramenta “Criar raster virtual” nas ferramentas de processamento do QGIS. Aqui, indicamos os 2 ficheiros que queremos incluir no VRT. A vantagem é que em vez de esperarmos que quase 2GB de dados sejam mastigados para criar um novo ficheiro de mais 2GB, vamos apenas criar muito rapidamente um pequeno ficheiro de texto com 2kb:

O resultado é um ficheiro que se comporta como uma só imagem, que no fundo serve apenas de apontador para 2 ficheiros TIF originais:

De seguida, vamos reprojetar este mosaico virtual para o sistema de coordenadas português.

2. e 3. Cortar e Reprojetar

Bom, reprojetar em QGIS é facílimo. Mas para complicar um bocadinho, vamos cortar e reprojetar num só passo! ahpoisé…

Podemos usar a ferramenta de corte para reprojetar porque tem lá a opção para escolher um sistema de coordenadas diferente para o ficheiro de saída. Assim, é escusado executar depois outra ferramenta só para a reprojeção.

Em primeiro lugar, temos de arranjar o polígono de Portugal continental para cortar o nosso mosaico virtual. A CAOP está disponível para download no site da DGT.

Vê-se bem que dei um toquezinho nos limites de Portugal – o que fiz foi aumentar um pouco a parte litoral para que não deixe nada de fora – é que o eu-dem não tem exatamente os mesmos contornos que a nossa CAOP, como seria de esperar. Por isso, ao aumentar a parte litoral quis evitar deixar de fora algumas partes do eu-dem.

Para fazer o corte do raster com um ficheiro vetorial, usamos a ferramenta “Recortar raster pela camada de máscara”. Aqui vamos usar alguns truques que explico a seguir:

A maioria das opções usadas não têm nada de especial: escolhemos o ficheiro a cortar (mdt), o ficheiro de corte (portugal), uma compressão elevada, e o ficheiro de saída. Mas algumas opções são menos óbvias:

• Na compressão elevada, podemos mudar o parâmetro “PREDICTOR” para o valor 3, que é mais eficiente com valores decimais (que é o nosso caso);
• Nos parâmetros adicionais, vamos incluir “-co tiled=yes” para produzir um ficheiro organizado por blocos comprimidos (isto é importante porque permite melhor compressão e desempenho);
• Nos parâmetros adicionais vamos também aproveitar para dar mais memória ao processamento (é isso que fazem os parâmetros GDAL_CACHEMAX e o WM).

Estas opções em conjunto fazem o processo mais rápido e permitem obter um ficheiro mais pequeno.

O ficheiro resultante tem as seguintes características:

• Projeção – ETRS89 PT TM06 (3763)
• Resolução espacial: 24,99m
• Largura e Altura: 11.633×23.220 pixeis (290kmx580km)
• Tipo de pixel: 32bit
• Compressão: Deflate, Level 9, Predictor 3
• Pirâmides: Não
• Valor “nulo” (no data): -3.40282e+38
• Tamanho total em disco: 346,39 MB

E pronto, o resultado é um mds/mdt para Portugal, já no nosso sistema de coordenadas. Para visualizar, alteramos o sistema de coordenadas do projeto para ETRS89 PT TM06, e vemos que a forma e orientação já é familiar:

Podemos melhorar um pouco o estilo e clareza da visualização alterando a simbologia, e também podemos criar pirâmides/overviews/”vistas gerais” para acelerar a visualização a qualquer escala.

Pirâmides/overviews/”vistas gerais”

As “vistas gerais”, tradução de pirâmides ou overviews, também são referidas neste blog exaustivamente… são versões de resoluções cada vez menores e rápidas (pixel cada vez maior), que tornam a visualização dos dados extremamente rápida – o software escolhe automaticamente que resolução mostra para cada escala.

O costume é usar compressão jpeg, mas neste caso não é possível porque não é compatível com dados 32 bit. Isto foi um pretexto para ir ver as novidades em termos de compressão, e há muitas. Um artigo excelente que compara os novos métodos LERC e ZSTD é o “Guide to GeoTIFF compression and optimization with GDAL” do Koko Alberti – tem no blog mais gemas preciosas que vale a pena ler…

Então, olhando para este artigo vemos que a compressão mais compacta e mais rápida para ficheiros de 32bit é o ZSTD com predictor 3, muitíssimo superior à compressão Deflate e ainda por cima muito mais rápida de visualizar!! Mas mais à frente neste artigo o Koko olha para uma compressão que degrada os dados, a LERC (Limited Error Raster Compression), que rebenta com a escala de compressão… claro que perde dados, e por isso compara-se diretamente com jpeg. Mas tem a vantagem de podermos definir o erro máximo que queremos tolerar. Sendo para pirâmides que servem apenas para visualizar, será uma excelente opção! Qual é o problema?? O GDAL permite criar pirâmides neste formato mas não permite que se defina os parâmetros de compressão, e por isso, não conseguimos definir o erro admissível… para já está fora de ação…

Então, podemos criar as vistas gerais (não me habituo a esta tradução…) usando a compressão ZSTD com predictor=3:

É produzido um ficheiro .ovr com o tamanho de 126MB. Para os dados em causa, esta compressão foi praticamente igual à “High compression” do QGIS… não sei se vale a pena.

E terminámos… o restante é só para nerds… o link para download está no fim do artigo.

Só por carolice, tentei comprimir este ovr com compressão LERC_DEFLATE e um erro de 0,5m, depois adicionei-lhe pirâmides, e obtive um tamanho de 66MB, perto de metade!!!! Como é apenas para visualização, é uma opção incrível quando não conseguimos usar JPEG. (Mas não se esqueçam que ao fazer identify no mdt, vão obter o valor visível, logo se querem ver a cota “verdadeira” têm de fazer zoom até verem a imagem base e não as pirâmides!)

Converter de 32bit para inteiros 16bit

Como sou do século passado e acho que espaço em disco é valioso, vou converter de 32bit para 16bit, porque as elevações no nosso território variam de -5 a 2000 e qualquer coisa metros. Estes valores cabem perfeitamente em 16bit. Mas como 16bit têm de ser inteiros (porque no qgis não há reais de 16bit), também vamos deitar fora as casas decimais… para um erro médio quadrático de 7m, e para os usos em causa, parece-me que não é grave deitar fora ou ganhar até 0,5m de elevação.

Sendo sincero, eu costumo fazer isto na linha de comando, com o GDAL. Mas o desafio aqui era fazer tudo em QGIS, que tem funções próprias e também acesso a todas as ferramentas do GDAL, disponibilizando uma interface de utilizador simpática.

Há várias formas de fazer esta conversão, mas temos de ter cuidado porque não queremos alterar os valores. No nosso caso, queremos preservar os valores de altimetria e passá-los apenas a inteiros de 16bit. Isto porque em 16bit os valores podem variar entre -32768 e +32767 , portanto mais que suficiente para guardar os nossos valores. O único valor que será convertido será o valor nulo que passará automaticamente a -32768. É que há vários algoritmos que recalculam os valores para que fiquem proporcionalmente na mesma “posição” dentro do novo domínio de valores possíveis. Por exemplo, se passarmos de 32bit para 16bit estes algoritmos partem do princípio que existem dados fora do intervalo válido de 16bit e por isso usam uma interpolação linear para os colocar neste intervalo. Assim, um valor de 2147483647 passaria a 32768, e por aí fora.

Então para manter os valores, a ferramenta que usei foi a Converter (Translate): https://docs.qgis.org/3.22/en/docs/user_manual/processing_algs/qgis/rasteranalysis.html#round-raster. Esta ferramenta é apenas uma interface para o comando GDAL – gdal_translate, que permite uma imensidão de parâmetros, mas nós queremos apenas fazer 2 coisas: passar a 16bit e manter a compressão do original:

Produzimos assim um ficheiro de apenas 68MB… Mais as pirâmides com apenas 30MB. Faz alguma diferença… como apenas removemos a parte decimal das cotas, o erro máximo é inferior +-0,5m:

Conclusão

Bem, aquilo que era para ser um simples artigo para disponibilizar um ficheiro de altimetria para Portugal, acabou por ficar um bocadinho extenso 😉

O ficheiro zipado 32bit está aqui: http://192.168.1.84:8000/eu-dem_portugal_32bit.zip.

E para quem quiser algo mais pequeno e não se importar de usar cotas inteiras (sem parte decimal) tem aqui a versão 16bit: http://192.168.1.84:8000/eu-dem-portugal-16bit.zip.

Até breve, e qualquer coisa deixem comentário.