Raw Sentinel 2 jp2 vers géotiff RGB

Je cherche un moyen de fusionner les fichiers de bande Sentinel 2 jp2 ( B02, B03, B04 ) et de corriger les couleurs RVB. Tout doit être fait avec un script bash ou python. Pour mon exemple, je travaille sur ces images . Idéalement, la solution sera proche de ce tutoriel.

Je peux fusionner les groupes avec cette commande

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

Mais pour une raison quelconque, je ne peux pas corriger les couleurs RVB avec la commande imagemagic. La sortie est une image noire de ~ 700 Mo.

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

Finalement, j'aimerais avoir un fichier géotiff pour le télécharger sur mapbox. Une explication sur la façon de choisir les convertparamètres est la bienvenue.

Je développe une application qui devrait deviner quelle partie de l'image satellite est une terre agricole. Une image de scène sera découpée en plus petits patchs (peut-être 64x64) et sera classée par CNN ( recadrée ou non recadrée ). J'utilise cet ensemble de données pour former le modèle Inception-v3. L'ensemble de données contient des images RVB 64x64 avec une résolution spatiale de 10 m.

Plus d'infos sur 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

Ceci est merged.tif avant et après l'application de la solution @ ben

Il y a 2 parties du problème. Le premier est que vous voulez convertir de 16 bits en 8 bits, et l'option -scale de gdal_translate le fait, comme mentionné dans la réponse précédente.

 -scale minOriginal maxOriginal minOutput maxOutput  

Le deuxième problème est un problème d'amélioration du contraste: lorsque vous redimensionnez, vous souhaitez avoir un contraste élevé pour les pixels qui vous intéressent. AVERTISSEMENT: Il n'y a pas de contraste "magique" car, lorsque vous redimensionnez, vous perdez généralement des informations : il est fait pour améliorer la visualisation des données, et les logiciels professionnels le font à la volée sans écrire de nouveau fichier. Si vous souhaitez poursuivre le traitement de vos données, votre géotiff "noir" contient les mêmes informations que votre jp2 et est prêt à être traité. Si vous calculez, par exemple, l'indice de végétation, cela devrait être fait avec les valeurs de réflectance "d'origine", pas celles remises à l'échelle. Cela étant dit, voici quelques étapes pour créer une image 8 bits visuellement améliorée.

@ben vous a donné une méthode générique pour redimensionner la réflectance de 0-1 (multiplié par 10000 avec ce produit) à 0-255. C'est sûr (pas d'exclusion), mais seuls les nuages ​​et certains sols nus ont des réflectances vraiment élevées, donc vous ne voyez pas beaucoup sur terre (sauf les sols nus) et rien dans l'eau. Par conséquent, les améliorations de contraste couramment appliquées aux images consistent à ne prendre qu'un sous-ensemble de la gamme complète. Du côté sûr, vous pouvez utiliser la connaissance que la réflectance maximale des matériaux de surface de la Terre communs est généralement inférieure à 0,5 / 0,6 (voir icipour quelques exemples). Bien sûr, cela suppose que votre image a été corrigée atmosphérique (images L2A). Cependant, la plage de réflectance diffère dans chaque bande spectrale et vous n'avez pas toujours les surfaces terrestres les plus brillantes dans votre zone d'intérêt. Voici à quoi ressemble la méthode "sûre" (avec une réflectance maximale de 0,4, comme la 4096 suggérée par @RoVo)

En revanche, le contraste pourrait être optimisé pour chaque bande. Vous pouvez définir cette plage manuellement (par exemple, vous êtes intéressé par la couleur de l'eau et vous connaissez la valeur de réflectance maximale attendue de l'eau) ou en fonction des statistiques de l'image. Une méthode couramment utilisée consiste à conserver environ 95% des valeurs et à "rejeter" (trop sombre -> 0 ou trop clair -> 255) le reste, ce qui revient à définir la plage en fonction de la valeur moyenne +/- 1,96 * écart-type. Bien sûr, ce n'est qu'une approximation car il suppose une distribution normale, mais cela fonctionne assez bien dans la pratique (sauf lorsque vous avez trop de nuages ​​ou si les statistiques utilisent certaines valeurs NoData).

Prenons l'exemple de votre premier groupe:

moyenne = 320

std = 536

Intervalle de confiance à 95% = [-731: 1372]

mais bien sûr, la réflectance est toujours supérieure à zéro, vous devez donc définir le minimum à 0.

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

Et si vous avez une version récente de gdal, vous pouvez utiliser -scale_ {band #} (0 255 est la sortie par défaut, donc je ne la répète pas) afin que vous n'ayez pas besoin de diviser des bandes uniques. J'ai également utilisé vrt au lieu de tif comme fichier intermédiaire (pas besoin d'écrire une image complète: une image virtuelle suffit)

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

Notez que vos statistiques sont fortement affectées par des "artefacts" tels que les nuages ​​et NoData. D'un côté, la variance est surestimée lorsque vous avez des valeurs extrêmes. D'un autre côté, votre moyenne est sous-estimée quand il y a une grande quantité de valeurs "zéro" (rendant l'image contrastée automatiquement trop lumineuse comme dans l'exemple) et elle serait surestimée s'il y avait une majorité de nuages ​​(ce qui rendrait le image trop sombre). À ce stade, les résultats ne seraient donc pas les meilleurs que vous puissiez obtenir.

Une solution automatisée consisterait à définir les valeurs d'arrière-plan et de cloud sur "nodata" et à calculer vos statistiques sans NoData (voir cet article pour plus de détails sur le calcul des statistiques sans NoData, et celui-ci pour un exemple pour définir des valeurs supérieures à 4000 à NoData également. ). Pour une seule image, je calcule généralement les statistiques sur le plus grand sous-ensemble possible sans nuage. Avec les statistiques d'un sous-ensemble où il n'y a pas de "NoData" (en haut à gauche de votre image), cela donne le résultat final. Vous pouvez voir que la plage est environ la moitié de la plage "sûre", ce qui signifie que vous avez deux fois plus de contraste:

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

Comme dernière remarque, gdal_constrast_stretch a l' air bien mais je n'ai pas testé

Vous pouvez simplement utiliser le TCI.jp2fichier qui est inclus dans les SAFE.zipfichiers. Notez que ces fichiers ne sont pas disponibles dans les fichiers S2 avant octobre 2016

Vous pouvez également convertir les bandes à l'aide de 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 4096est une valeur raisonnable pour les scènes Sentinel-2 et afaik également utilisé pour les images TCI.jp2. Abaissez le 4096 si vous souhaitez obtenir un résultat plus clair.

Si vous cherchez une solution comme celle que vous avez liée à la question, vous devez suivre et ajuster le skript shell de traitement Landsat 8 qui est fourni pour téléchargement dans le tutoriel.

En particulier, comme cela se fait là-bas, vous voudrez peut-être d'abord redimensionner les bandes individuelles, par exemple comme suit:

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

Notez que l'histogramme de vos images suggère que vous n'avez que des surfaces très sombres dans votre image (est-ce le cas?) Mais généralement votre image sentinelle-2 sera la réflectance en haut de l'atmosphère ou en surface où les valeurs varient généralement entre 0 et 10000 - sauf si des valeurs plus élevées sont également possibles, par exemple si vous avez des nuages ​​dans l'image.

Ensuite, vous pouvez fusionner les bandes et affiner l'apparence de l'image:

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

Voici ce qui arrive à mon image en faisant cela: