Geotiff Raw Sentinel 2 jp2 para RGB

Estou procurando uma maneira de mesclar arquivos de banda jp2 do Sentinel 2 ( B02, B03, B04 ) e corrigir cores RGB. Tudo deve ser feito com scripts bash ou python. Para o meu exemplo, trabalho nessas imagens . Idealmente, a solução estará próxima deste tutorial.

Eu sou capaz de mesclar as bandas com este comando

gdal_merge.py -separate -co PHOTOMETRIC=RGB -o merged.tif B04.jp2 B03.jp2 B02.jp2

Mas, por algum motivo, não consigo corrigir as cores RGB com o comando imagemagic. A saída é ~ 700MB de imagem em preto.

convert -channel B -gamma 1.05 -channel RGB -sigmoidal-contrast 20,20% -modulate 100,150 merged.tif merged-cc.tif

Eventualmente, eu gostaria de ter um arquivo geotiff para enviá-lo no mapbox. Explicação de como se deve escolher convertparâmetros é bem-vinda.

Estou desenvolvendo um aplicativo que deve adivinhar qual parte da imagem de satélite é de terras agrícolas. Uma imagem de cena será cortado em pequenos pedaços (talvez 64x64) e vai ser classificadas pela CNN ( colheita ou não-colheita ). Eu uso esse conjunto de dados para treinar o modelo Inception-v3. O conjunto de dados contém imagens RGB de 64x64 com resolução espacial de 10m.

Mais informações sobre merged.tif

Band 1 Block=10980x1 Type=UInt16, ColorInterp=Red
  Metadata:
    STATISTICS_MAXIMUM=4818
    STATISTICS_MEAN=320.61101402206
    STATISTICS_MINIMUM=0
    STATISTICS_STDDEV=536.76609312554
Band 2 Block=10980x1 Type=UInt16, ColorInterp=Green
  Metadata:
    STATISTICS_MAXIMUM=4206
    STATISTICS_MEAN=350.98505344194
    STATISTICS_MINIMUM=0
    STATISTICS_STDDEV=534.43264268631
Band 3 Block=10980x1 Type=UInt16, ColorInterp=Blue
  Metadata:
    STATISTICS_MAXIMUM=3801
    STATISTICS_MEAN=364.44611471973
    STATISTICS_MINIMUM=0
    STATISTICS_STDDEV=544.55509661709

Isso é mesclado.tif antes e depois de aplicar a solução da @ ben

Existem 2 partes do problema. A primeira é que você deseja converter de 16 bits para 8 bits, e a opção -scale de gdal_translate faz isso, conforme mencionado na resposta anterior.

 -scale minOriginal maxOriginal minOutput maxOutput  

O segundo problema é um problema de aprimoramento de contraste: quando você redimensiona, deseja ter um alto contraste para os pixels em que está interessado. AVISO: Não há contraste "mágico" porque, quando você redimensiona, geralmente perde algumas informações : isso é feito para melhorar a visualização dos dados, e os softwares profissionais fazem isso rapidamente, sem precisar escrever um novo arquivo. Se você deseja continuar processando seus dados, seu geotiff "preto" contém as mesmas informações que seu jp2 e está pronto para ser processado. Se você calcular, por exemplo, o índice de vegetação, isso deve ser feito com os valores de refletância "originais", não com os redimensionados. Dito isto, aqui estão algumas etapas para criar uma imagem de 8 bits visualmente aprimorada.

O @ben forneceu um método genérico para redimensionar a refletância de 0-1 (multiplicado por 10000 com este produto) para 0-255. Isso é seguro (sem exclusão), mas apenas nuvens e alguns solos nus têm refletâncias realmente altas, para que você não veja muito em terra (exceto solos nus) e nada na água. Portanto, os aprimoramentos de contraste comumente aplicados nas imagens consistem em capturar apenas um subconjunto de toda a faixa. No lado seguro, você pode usar o conhecimento de que a refletância máxima do material comum da superfície da Terra geralmente é inferior a 0,5 / 0,6 (veja aquipara alguns exemplos). Obviamente, isso pressupõe que sua imagem foi corrigida atmosférica (imagens L2A). No entanto, o alcance da refletância difere em cada banda espectral e você nem sempre tem as superfícies mais brilhantes da Terra em sua área de interesse. Aqui está a aparência do método "seguro" (com uma refletância máxima de 0,4, como o 4096 sugerido por @RoVo)

Por outro lado, o contraste pode ser otimizado para cada banda. Você pode definir esse intervalo manualmente (por exemplo, se interessa pela cor da água e conhece o valor máximo esperado de refletância da água) ou com base nas estatísticas da imagem. Um método comumente usado consiste em manter aproximadamente 95% dos valores e em "descartar" (muito escuro -> 0 ou muito claro -> 255) o restante, o que é semelhante à definição do intervalo com base no valor médio +/- 1,96 * desvio padrão. Obviamente, isso é apenas uma aproximação, pois assume uma distribuição normal, mas funciona muito bem na prática (exceto quando você tem muitas nuvens ou se as estatísticas usam alguns valores NoData).

Vamos pegar sua primeira banda como exemplo:

média = 320

std = 536

Intervalo de confiança de 95% = [-731: 1372]

mas é claro que a refletância é sempre maior que zero, portanto, você deve definir o mínimo em 0.

gdal_translate -scale 0 1372 0 255 -ot Byte  B01.jp2 B01-scaled.tif  

E se você possui uma versão recente do gdal, pode usar -scale_ {band #} (0 255 é a saída padrão, portanto não a repito) para que você não precise dividir as bandas únicas. Também usei vrt em vez de tif como um arquivo intermediário (não há necessidade de escrever uma imagem completa: basta uma virtual)

gdalbuildvrt -separate stack.vrt B04.jp2 B03.jp2 B02.jp2
gdal_translate -scale_1 0 1372 -scale_2 0 1397 -scale_3 0 1430 -ot Byte  stack.vrt im_rescaled.tif

Observe que suas estatísticas são fortemente afetadas por "artefatos", como nuvens e NoData. Por um lado, a variação é superestimada quando você tem valores extremos. Por outro lado, sua média é subestimada quando há uma grande quantidade de valores "zero" (tornando a imagem contrastada automaticamente muito brilhante como no exemplo) e seria superestimada se houvesse uma maioria de nuvens (o que tornaria a imagem muito escura). Nesse estágio, os resultados não seriam os melhores que você poderia obter.

Uma solução automatizada seria definir valores de segundo plano e nuvem como "nodata" e calcular suas estatísticas sem o NoData (consulte esta postagem para obter detalhes sobre estatísticas de computação sem o NoData, e este por exemplo, para definir valores maiores que 4000 no NoData também ) Para uma única imagem, costumo calcular as estatísticas do maior subconjunto possível de nuvens. Com estatísticas de um subconjunto em que não há "NoData" (canto superior esquerdo da sua imagem), isso fornece o resultado final. Você pode ver que o intervalo é cerca da metade do intervalo "seguro", o que significa que você tem o dobro de contraste:

gdal_translate -scale_1 38 2225 -scale_2 553 1858 -scale_3 714 1745 -ot Byte  stack.vrt im_rescaled.tif

Como última observação, gdal_constrast_stretch parece bom, mas eu não testei

Você pode simplesmente usar o TCI.jp2arquivo que está incluído nos SAFE.ziparquivos. Observe que esses arquivos não estão disponíveis nos arquivos S2 antes de outubro de 2016

Como alternativa, você pode converter as bandas usando GDAL:

# Merge bands
gdalbuildvrt -separate TCI.vrt B04.jp2 B03.jp2 B02.jp2

# Convert to uncompressed GeoTiff
gdal_translate -ot Byte -co TILED=YES -scale 0 4096 0 255 TCI.vrt TCI.tif

# _OR_ Convert to JPEG - compressed GeoTiff
gdal_translate -ot Byte -co TILED=YES -co COMPRESS=JPEG -co PHOTOMETRIC=YCBCR -scale 0 4096 0 255 TCI.vrt TCI.tif

-scale 0 4096é um valor razoável para cenas do Sentinel-2 e também é usado para as imagens TCI.jp2. Abaixe o 4096 se desejar receber um resultado mais claro.

Se você está procurando uma solução como a que você vinculou na pergunta, deve seguir e ajustar o skript do shell de processamento do Landsat 8, fornecido para download no tutorial.

Particularmente, como é feito lá, primeiro você pode querer redimensionar as bandas únicas, por exemplo, da seguinte maneira:

gdal_translate -ot Byte -scale 0 10000 0 255 B04.jp2 B04-scaled.tif 
gdal_translate -ot Byte -scale 0 10000 0 255 B03.jp2 B03-scaled.tif
gdal_translate -ot Byte -scale 0 10000 0 255 B02.jp2 B02-scaled.tif

Observe que o histograma de suas imagens sugere que você só possui superfícies muito escuras na imagem (é esse o caso?), Mas geralmente sua imagem sentinela-2 será refletida no topo da atmosfera ou na superfície, onde os valores geralmente variam entre 0 e 10000 - a menos que valores mais altos também sejam possíveis, por exemplo, se houver nuvens na imagem.

Em seguida, você pode mesclar as bandas e ajustar a aparência da imagem:

gdal_merge.py -v -ot Byte -separate -of GTiff -co PHOTOMETRIC=RGB -o RGB-scaled.tif B04-scaled.tif B03-scaled.tif B02-scaled.tif
convert -channel B -gamma 1.05 -channel RGB -sigmoidal-contrast 20,40% -modulate 100,150 RGB-scaled.tif RGB-scaled-cc.tif

Isto é o que acontece com a minha imagem ao fazer isso: