Versão abreviada…
Podem aceder aqui a um motor de mosaicos de imagens Sentinel-2 RGB e IRG para Portugal, com serviço WMS, com suporte temporal… O serviço WMS está funcional, mas para usar em QGIS é preciso algumas definições (ler abaixo, muuuito abaixo).
Alguns avisos: isto é um projecto pessoal, de carolice, tem muitos defeitos, eu sei, que podem ou não vir a ser resolvidos… Estou muito interessado em ouvir sugestões, e para isso nada melhor que o twitter ou os comentários aqui no blog.
Se a carga for demasiada no servidor, os pedidos são “desacelerados”, por forma a manter o servidor equilibrado. Por favor, não usem scripts de download… Pretendo incluir a função de download em breve. Entretanto, se precisarem de alguma imagem é só dizer, eu farei os possíveis para responder atempadamente.
E pronto, agora quem tem curiosidade e paciência pode continuar a ler…
Introdução
Acho as imagens Sentinel um prodígio, a sério. Temos imagens de 10m de resolução, 7 em 7 dias, para grande parte do planeta, gratuitas! É espantoso…
Tenho usado imagens que pesquiso e descarrego a partir dos sites de distribuição da ESA. Mas é um processo moroso obter as imagens – cada uma é 1GB num zip com 13 ficheiros jpeg2000, e mais uma dezena de outros ficheiros de metadados. Ainda mais moroso se quiser juntar algumas imagens. Mas é tudo gratuito.
Também consulto visualizadores web com imagens processadas, com equilibrio de cores, e sem nuvens, e com várias combinações de bandas. É só pesquisar… também é gratuito.
Agora, se subirmos um pouco o grau de exigência e quisermos sobrepor a nossa informação às imagens, podemos usar um serviço WMS. Mas aqui já é pago, e começa nos 20€/mês por utilizador particular ou académico. Se for empresa o preço já sobe. E é completamente justificado. Mas também é um pouco contraditório em relação ao que se pretende do programa Copernicus, que seria disseminar o mais possível os resultados do programa na sociedade civil… 20€/mês parece-me um pouco contraditório.
Na mais recente incursão pelos sites de descarga das imagens, comecei a fazer alguns scripts, muito simples. Mas a coisa foi-se alongando, e acabei por automatizar grande parte do que seria um servidor wms, com dados processados periodicamente. Foi assim que nasceu este projecto – WMS Time Machine!
Isto não é bem uma invenção… Os nossos vizinhos aqui ao lado do Institut Cartogràfic i Geològic de Catalunya têm um servidor WMS-T de imagens Sentinel *do caraças!*…
http://www.icgc.cat/en/Public-Administration-and-Enterprises/Services/Online-services-Geoservices/WMS-Ortoimatges/WMS-Sentinel-2-orthoimages
É pena ser só para a Catalunha…
Componentes e Processo geral
Os passos do processo e os componentes de software usado são:
- Pesquisa e descargas de imagens – Sentinelsat
- Processamento das imagens – GDAL/OGR
- Manutenção da BD de imagens – GDAL/OGR
- Serviço WMS Time – MapServer
- Visualizar o serviço em QGIS
- Visualizador Web – OpenLayers (Feito à pressa! No futuro, talvez um de jeito 😉
Pesquisa e descarga das imagens – Sentinelsat
A pesquisa das imagens e a filtragem das que interessam é feita usando a biblioteca python SentinelSat, que também inclui uma CLI. São apenas usadas imagens Sentinel-2, nível 2A (significa que têm já correção atmosférica).
O esquema de pesquisa é muito simples, e estou a procurar formas de o melhorar. É indicada uma data, e são pesquisadas 16 quadrículas ou tiles ou grânulos (não tenho paciência para ler o dicionário do Sentinel) que correspondem a Portugal, que tenham sido obtidas até 5 dias atrás, com nuvens até 70%. Também são recusadas imagens com menos de 75% do seu tamanho normal para evitar imagens com grandes áreas sem dados. Alguns mosaicos podem ter sido construídos manualmente com nuvens até 90%…
Um exemplo de pesquisa com Sentinelsat:
sentinelsat --start 20190818 --end 20190819 --sentinel 2 --instrument MSI --producttype S2MSI2A --cloud 40 --query "filename=/.+29S(N|P)C.+/" --user bla --password blabla
Este comando procura ficheiros com nomes que incluam “29SNC” ou “29SPC” entre os dias 18 e 19-08-2019, com cobertura de nuvens até 40%. A facilidade de usar expressões regulares é muito flexível. Os docs do Sentinelsat e da interface de pesquisa da ESA são muito bons.
Só são descarregadas 4 bandas das 13 disponíveis: B02 (Red), B03 (Green), B04 (Blue), e B08 (Infrared). Isto significa que para uma data são descarregados 64 ficheiros .jp2 (16 tiles x 4 bandas), cerca de 7GB no total. As 4 bandas são depois combinadas em combinações RGB e IRG.
O processo de pesquisa deveria ser melhorado… Gostava de evitar mosaicos com partes significativas sem dados (neste momento há tiles com 25% de área sem dados). Estou a pensar numa pesquisa sobre um maior período, ordenar as imagens por qualidade, e detectar que datas apresentam melhores coberturas. Ou seja, um processo quase inverso do actual… Outra opção é usar o footprint de cada imagem e analisar geometricamente qual a combinação com menores vazios, ou mesmo só fazer o mosaico se não houver vazios.
Bom tile 
Mau tile
Processamento das imagens – GDAL
Depois de descarregados, os ficheiros são processados com GDAL e OGR. As 4 bandas de cada quadrícula ou grânulo são convertidas para 8bit, sendo criadas uma imagem virtual RGB e outra IRG para cada Tile. Depois é aplicado um stretch “virtual” para melhorar o contraste. São criados tileindexes para reunir todas as imagens RGB e IRG do país para uma data. Em seguida tento explicar melhor…
Imagem IRG 
Imagem RGB
A compressão faz com que o tamanho das bandas passe de ~100MB (jp2000, 16bit) para ~6MB (tiff/jpeg, 8bit). E assim os 64 ficheiros tif só ocupam 600MB para todo Portugal Continental, numa data.
É importante referir que para passar de 16bit para 8bit usei a forma mais simples do parâmetro -scale, que efectivamente passa os valores min-max dos 16bits da imagem, para valores entre 0-255. Isto efectivamente já provoca uma alteração à cor e contraste da imagem. E, obviamente, alguma perda de informação. Estas imagens são apenas para visualização e não se recomendam para análise.
Foi necessário melhorar ainda mais o contraste para que as imagens não fiquem demasiado escuras. Como o MapServer não suporta grande coisa na simbologia de rasters, tive de resolver ao nível dos dados.
Encontrei muita informação sobre a opção scale, mas poucas soluções para ajuste do contraste. Algumas óptimas soluções permitem criar novas imagens melhoradas, mas obrigam a mais espaço em disco e mais tempo de processamento…
A solução foi criar um .vrt que faz um ajuste ao contraste dinamicamente através da opção -scale. É usado o método de ajuste pelo desvio padrão, aplicando-o a cada banda. Ou seja, em cada banda é obtida a média e o desvio padrão, e o ajuste é feito calculando novos mínimos e máximos (fazendo o mesmo que uma das opções de Contrast Enhancement do QGIS):
-scale 0 255 media-2.8*stdev media+2.8*stddev
Funciona muito bem, desde que não existam nuvens na imagem:
Original 8bit 
Contraste melhorado
Com nuvens, e respectivas sombras, tudo piora, como é de esperar… será necessário melhorar o processo para evitar estas áreas brancas e negras:
Criar base de dados das datas e imagens – OGR
Numa data temos assim 16 imagens RGB, e 16 imagens IRG. Todas são virtuais (.vrt), apens combinam 3 das 4 bandas, e não ocupam espaço. Agora queremos criar uma listagem que indique que datas já recolhemos, e quais as imagens que constituem cada data.
A base de de dados é apenas um shapefile com as quadrículas de todas as imagens vrt… o processo é simples e consiste em criar tileindexes… uma forma anciã de ver mosaicos de imagens e que ainda funciona em Mapserver.
Começamos por criar um tileindex das imagens numa data. Simples comando de gdaltindex. Um exemplo de índice rgb do mosaico para o dia 30/06/2019:
Índice com datas 
Índice com ficheiros
Na verdade, os índices são geograficamente todos iguais, porque uso sempre as mesmas 16 tiles. Mas o ideal seria pesquisar pela área de Portugal, e ver que tiles têm melhor cobertura na data escolhida.
A única coisa que varia entre datas são os ficheiros de imagem que são descarregados.
Bom, já temos um tileindex para cada dia descarregado, que é um shapefile com um atributo a indicar o caminho para cada imagem.
Como fazer uma base de dados de todos os mosaicos que já foram descarregados e existem no servidor? A resposta é simples: copiamos estes registos para um shapefile global usando ogr2ogr com a opção -update. E sempre que se constrói um mosaico para uma data nova, vamos inserir estes registos no shapefile global.
Aqui é usada a função SQL do OGR, que é absolutamente fantástica… permite executar SQL ao carregar dados para um shapefile, ou qualquer outro formato.
Assim, ao copiar as quads de um mosaico para o índice global aproveitamos para actualizar alguns campos extra:
- data do mosaico (campo time)
- data da imagem (campo dataimg)
- nome do ficheiro (campo location)
- nome do ficheiro com contraste melhorado (campo localviz)
Significa que sabemos as datas todas que recolhemos, e as imagens que as compõem. Tudo com shapefile!! (o shapefile é eterno!)
E é compatível com MapServer…
Só um exemplo do comando ogr para apagar imagens que já existam de uma data:
ogrinfo -dialect SQLITE tileindex_global_irg.shp -sql "DELETE FROM tileindex_global_irg where location like '%_20190830%_irg_%.vrt'"
Publicar as imagens num serviço WMS-Time
Usei o MapServer como servidor WMS com suporte do parâmetro Time.
O MapServer é fácil de configurar e manipular apenas com ficheiros de texto. A sua arquitectura é tão simples que apenas é necessário um nginx para o colocar na net. A exigência de memória é também muito reduzida. E claro, é um óptimo amigo do GDAL/OGR. Tudo o que era preciso…
Assim, foi criado um mapfile único com 4 layers:
- Índice das imagens, com label mostrando a data de cada imagem
- Índice das imagens, com label mostrando o nome do ficheiro (para vermos a tile respectiva se quisermos obter o original no site da ESA)
- Mosaico RGB
- Mosaico IRG
Todos apontam para o tileindex RGB ou para o tileindex IRG. Os índices são vectoriais, e os mosaicos são rasters. Simples.
O parâmetro TIME permite filtrar os dados para só mostrar aqueles que cumprem essa query. Ou seja, passamos uma data e o servidor devolve uma imagem onde todos os layers com TIME definido são filtrados por essa data.
No nosso caso, o campo usado para o filtro de data é o campo time, que indica a data de construção do mosaico. Por exemplo, este request mostra só imagens do mosaico RGB com data de 2019-06-30:
http://sentinelpt.viasig.com/wms/sentinelpt/?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetMap&FORMAT=image%2Fjpeg&TRANSPARENT=true&LAYERS=SentinelPT_RGB%2CIndice&TIME=2019-06-30&CRS=EPSG%3A3857&STYLES=&WIDTH=866&HEIGHT=538&BBOX=-868964.8014314906%2C4637845.355088675%2C-860690.4931196203%2C4642985.745240854
Há uma limitação ainda por resolver no serviço WMS: as datas disponíveis podem ser anunciadas pelo serviço. Mas neste momento estão fixadas:
Como ver no QGIS
Pois é… o QGIS não tem grande suporte para usar WMS-Time… mas funciona com alguns truques – basta indicar a data que queremos no url do serviço, e ligar as opções para ignorar os url’s devolvidos no capabilities doc:
Visualizador web
Bom, este visualizador é muito básico. Foi feito com base neste viewer baseado em OpenLayers: https://www.earder.com/tutorials/timeseries-with-geoserver-and-openlayers/.
Permite selecionar a data, tema RGB ou IRG, com ou sem labels. E ver o link WMS correspondente à data selecionada:
http://sentinelpt.viasig.com/
Melhorias??
Tantas, tantas…
- Criar transparência onde não há dados (zonas negras)
- Criar serviço de download (WCS)?
- Melhorar o ajuste de contraste automático actual, de acordo com cada imagem ou fazendo um match dos histogramas (neste momento, é feito um stretch de desvio padrão em todas as bandas)
- Excluir as nuvens do ajuste de contraste
- Criar serviço de tiles (WMTS)?
- Usar um visualizador web como deve ser (TerriaJS?)
- Obter dados para 2017
- Download dos footprints das imagens e usá-los na pesquisa e processamento
- Selecionar entre tiles disponíveis com menor área sem dados
- Reduzir a compressão dos ficheiros, aumentando a qualidade
Acho que nem eu próprio li isto tudo… mas fica para cábula futura.
Legalidades
Este é um projecto pessoal, sem qualquer garantia. Portanto, use por sua própria conta e risco. Contém dados de imagem de satélite Copernicus Sentinel-2 para vários anos, processados para efeitos de visualização e arquivo. Os dados originais são disponibilizados gratuitamente pela União Europeia para todos os fins. Mais informação pode ser consultada aqui: https://scihub.copernicus.eu.
Os dados deste projecto são disponibilizados sob a licença Atribuição 4.0 Internacional (CC BY 4.0)(link). Em resumo, esta licença permite qualquer uso dos dados, mesmo comercialmente, desde que seja indicada a sua fonte e não sejam impostas restrições adicionais aos dados.
