Tag Archives: Sentinel

SentinelPT WMS Time Machine

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).

http://sentinelpt.viasig.com/

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!

Rick and Morty a melhor série de todos os tempos!
Melhor invenção de todos os tempos

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:

  1. Pesquisa e descargas de imagens – Sentinelsat
  2. Processamento das imagens – GDAL/OGR
  3. Manutenção da BD de imagens – GDAL/OGR
  4. Serviço WMS Time – MapServer
  5. Visualizar o serviço em QGIS
  6. 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.

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…

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:

Com nuvens, e respectivas sombras, tudo piora, como é de esperar… será necessário melhorar o processo para evitar estas áreas brancas e negras:

Falta de contraste na presença de nuvens

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:

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)
Exemplo de identify no Índice de imagens RGB

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:

  1. Índice das imagens, com label mostrando a data de cada imagem
  2. Índice das imagens, com label mostrando o nome do ficheiro (para vermos a tile respectiva se quisermos obter o original no site da ESA)
  3. Mosaico RGB
  4. 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:

Período anunciado no WMS-Time é fictício…

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:

Usar o serviço WMS-Time no QGIS

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.