Raw Sentinel 2 jp2 a geotiff RGB

Estoy buscando una manera de fusionar los archivos de banda Sentinel 2 jp2 ( B02, B03, B04 ) y corregir los colores RGB. Todo debe hacerse con script bash o python. Para mi ejemplo, trabajo en estas imágenes . Idealmente, la solución estará cerca de este tutorial.

Puedo 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. La salida es ~ 700 MB de imagen en negro.

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

Finalmente, me gustaría tener un archivo geotiff para cargarlo en mapbox. La explicación de cómo se deben elegir los convertparámetros es bienvenida.

Estoy desarrollando una aplicación que debería adivinar qué parte de la imagen de satélite es tierra agrícola. Una imagen de escena se cortará en parches más pequeños (quizás 64x64) y se clasificará por CNN ( recortada o no recortada ). Utilizo este conjunto de datos 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 acerca de 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

Esto se fusiona.tif antes y después de aplicar la solución de @ben

Hay 2 partes del problema. La primera es que desea convertir de 16 bits a 8 bits, y la opción de escala de gdal_translate lo hace, como se mencionó en la respuesta anterior.

 -scale minOriginal maxOriginal minOutput maxOutput  

El segundo problema es un problema de mejora de contraste: cuando cambia la escala, desea tener un alto contraste para los píxeles que le interesan. ADVERTENCIA: No hay contraste "mágico" porque, cuando cambia la escala, generalmente pierde alguna información : se hace para mejorar la visualización de los datos, y los softwares profesionales lo hacen sobre la marcha sin escribir un nuevo archivo. Si desea procesar aún más sus datos, su geotiff "negro" contiene la misma información que su jp2 y está listo para ser procesado. Si calcula, por ejemplo, el índice de vegetación, esto debe hacerse con los valores de reflectancia "originales", no con los valores reescalados. Dicho esto, aquí hay algunos pasos para crear una imagen de 8 bits visualmente mejorada.

@ben le dio un método genérico para reescalar la reflectancia de 0-1 (multiplicado por 10000 con este producto) a 0-255. Esto es seguro (sin exclusión), pero solo las nubes y algunos suelos desnudos tienen reflejos muy altos, 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 comúnmente en las imágenes consisten en tomar solo un subconjunto del rango completo. En el lado seguro, puede usar el conocimiento de que la reflectancia máxima del material de la superficie de la Tierra común generalmente está por debajo de 0.5 / 0.6 (consulte aquípara algunos ejemplos) Por supuesto, esto supone que su imagen ha sido corregida atmosféricamente (imágenes L2A). Sin embargo, el rango de reflectancia difiere en cada banda espectral y no siempre tiene las superficies terrestres más brillantes en su área de interés. Así es como se ve el método "seguro" (con una reflectancia máxima de 0.4, como el 4096 sugerido por @RoVo)

Por otro lado, el contraste podría optimizarse para cada banda. Puede definir este rango manualmente (por ejemplo, si está interesado en el color del agua y conoce el valor máximo de reflectancia esperado del agua) o en función de las estadísticas de la imagen. Un método comúnmente usado consiste en mantener aproximadamente el 95% de los valores y "descartar" (demasiado oscuro -> 0 o demasiado brillante -> 255) el resto, que es similar a definir el rango basado en el valor medio +/- 1.96 * Desviación Estándar. Por supuesto, esto es solo una aproximación porque supone una distribución normal, pero funciona bastante bien en la práctica (excepto cuando tiene demasiadas nubes o si las estadísticas hacen uso de algunos valores NoData).

Tomemos su primera banda como ejemplo:

media = 320

std = 536

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

pero, por supuesto, la reflectancia siempre es mayor que cero, por lo tanto, debe establecer el mínimo en 0.

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

Y si tiene una versión reciente de gdal, puede usar -scale_ {band #} (0 255 es la salida predeterminada, por lo que no lo repito) para que no necesite dividir bandas individuales. También usé 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

Tenga en cuenta que sus estadísticas se ven fuertemente afectadas por "artefactos" como nubes y NoData. Por un lado, la varianza se sobreestima cuando tiene valores extremos. Por otro lado, su promedio se subestima cuando hay una gran cantidad de valores "cero" (haciendo que la imagen contrastada automáticamente sea demasiado brillante como en el ejemplo) y se sobreestimaría si hubiera una mayoría de nubes (lo que haría que imagen demasiado oscura). Por lo tanto, en esta etapa, los resultados no serían los mejores que podría obtener.

Una solución automatizada sería establecer valores de fondo y de nube en "nodata" y calcular sus estadísticas sin NoData (consulte esta publicación para obtener detalles sobre cómo calcular estadísticas sin NoData, y esta también como ejemplo para establecer valores superiores a 4000 en NoData ) Para una sola imagen, generalmente calculo las estadísticas en el subconjunto más grande posible sin nubes. Con estadísticas de un subconjunto donde no hay "NoData" (arriba a la izquierda de su imagen), esto da el resultado final. Puede ver que el rango es aproximadamente la mitad del rango "seguro", lo que significa que tiene el doble de contraste:

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

Como último comentario, gdal_constrast_stretch se ve bien pero no lo he probado

Simplemente puede usar el TCI.jp2archivo que está incluido en los SAFE.ziparchivos. Tenga en cuenta que estos archivos no están disponibles en archivos S2 antes de octubre de 2016

Alternativamente, puede convertir las 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 4096es un valor razonable para las escenas Sentinel-2 y afaik también se usa para las imágenes TCI.jp2. Baje el 4096 si desea recibir un resultado más ligero.

Si está buscando una solución como la que ha vinculado en la pregunta, debe seguir y ajustar el skript de shell de procesamiento Landsat 8 que se proporciona para descargar 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 solo tiene superficies muy oscuras en su imagen (¿es este el caso?), Pero generalmente su imagen centinela-2 será la reflectancia de la parte superior de la atmósfera o en la superficie, donde los valores comúnmente oscilan entre 0 y 10000, a menos que también sean posibles valores más altos, por ejemplo, si tiene nubes en la imagen.

Luego puede fusionar las bandas y ajustar la apariencia 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 sucede a mi imagen cuando hago esto: