Raw Sentinel 2 jp2 a RGB geotiff

Question

Estoy buscando una forma de fusionar Centinela 2 archivos de banda jp2( B02, B03, B04 ) y fijar los colores RGB. Todo debe hacerse con bash o python script. Para mi ejemplo trabajo en estas imágenes . Lo ideal es que la solución se acerque a este tutorial.

Soy capaz de fusionar las bandas con este comando

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

Pero por alguna razón no puedo arreglar los colores RGB con el comando imagemagic. El resultado es una imagen negra de ~700MB.

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

Eventualmente me gustaría tener un archivo geotiff para subirlo a mapbox. Explicación de cómo se debe elegir convert los parámetros son bienvenidos.

Estoy desarrollando una aplicación que debe adivinar qué parte de la imagen de satélite es tierra agrícola. Una imagen de la escena será cortada en parches más pequeños (tal vez 64x64) y será clasificada por CNN( cultivo o no cultivo ). Yo uso este para entrenar el modelo Inception-v3. El conjunto de datos contiene imágenes RGB de 64x64 con una resolución espacial de 10 m.

---

Más información 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

Este es merged.tif antes y después de aplicar la solución de @ben

Answer 1

Hay dos partes del problema. La primera es que quieres convertir de 16 bits a 8 bits, y la opción -scale de gdal_translate lo hace, como se menciona en la respuesta anterior.

 -scale minOriginal maxOriginal minOutput maxOutput  

El segundo problema es un problema de mejora del contraste: cuando se reescala, se quiere tener un alto contraste para los píxeles que interesan. No hay un contraste "mágico" porque, al reescalar, se pierde algo de información. @ben te dio un método genérico para reescalar la reflectancia de 0-1 (multiplicada por 10000 con este producto) a 0-255. Esto es seguro (sin exclusión), pero sólo las nubes y algunos suelos desnudos tienen reflectancias realmente altas, por lo que no se ve mucho en la tierra (excepto los suelos desnudos) y nada en el agua. Por lo tanto, las mejoras de contraste que se aplican habitualmente en las imágenes consisten en tomar sólo un subconjunto del rango completo. Se puede definir este rango manualmente (por ejemplo, si se está interesado en el color del agua y se conoce el valor máximo esperado de reflectancia del agua) o basándose en las estadísticas de la imagen. un uso común es entonces conservar aproximadamente el 95% de los valores y descartar el resto, lo que es similar a definir el rango basándose en el valor medio +/- 1,9*desviación estándar. De esto es sólo una aproximación porque asume una distribución normal, pero funciona bastante bien en la práctica (ecept cuando usted tiene demasiadas nubes de si las estadísticas utilizan algunos valores NoData.

Tomemos como ejemplo su primera banda:

media = 320

std = 536

Intervalo de confianza del 95% = [-731:1372]

pero, por supuesto, la reflectancia es siempre mayor que cero, por lo que se podría detectar el mínimo en 0.

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

Ten en cuenta que la varianza está sobrestimada debido a la presencia de nubes, por lo que tu imagen seguirá pareciendo un poco demasiado oscura. Yo suelo calcular las estadísticas sobre el mayor subconjunto posible sin nubes.

Y si tienes una versión reciente de gdal, entonces puedes usar -scale_{band#} (0 255 es la salida por defecto, así que no lo repito) para no tener que dividir bandas individuales. También he utilizado vrt en lugar de tif como archivo intermedio (no es necesario escribir una imagen completa: una virtual es suficiente)

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

Como última observación, gdal_constrast_stretch tiene buena pinta pero no lo he probado

Answer 2

Puede utilizar simplemente el TCI.jp2 que su incluyó en el SAFE.zip archivos. Tenga en cuenta que estos archivos no están disponibles en los archivos S2 antes de octubre de 2016

También puedes convertir las bandas con 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 es un valor razonable para las escenas de Sentinel-2 y, al parecer, también se utiliza para las imágenes TCI.jp2. Baje el valor de 4096 si quiere recibir un resultado más claro.

Answer 3

Si buscas una solución como la que enlazas en la pregunta debes seguir y ajustar el Esquema del procesamiento de Landsat 8 que se proporciona para su descarga en el tutorial.

En particular, como se hace allí, es posible que primero desee reescalar las bandas individuales, por ejemplo, de la siguiente manera:

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

Tenga en cuenta que el histograma de sus imágenes sugiere que sólo tiene superficies muy oscuras en su imagen (¿es este el caso?), pero normalmente su imagen sentinel-2 será la reflectancia de la parte superior de la atmósfera o de la superficie, donde los valores suelen oscilar entre 0 y 10000 - a menos que también sean posibles valores más altos, por ejemplo, si tiene nubes en la imagen.

A continuación, puedes fusionar las bandas y afinar el aspecto de la imagen:

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

Esto es lo que le ocurre a mi imagen al hacer esto: