Raw Sentinel 2 jp2 zu RGB Geotiff

Ich suche nach einer Möglichkeit, Sentinel 2- JP2-Banddateien ( B02, B03, B04 ) zusammenzuführen und RGB-Farben zu korrigieren . Alles sollte mit Bash- oder Python-Skript durchgeführt werden. In meinem Beispiel arbeite ich an diesen Bildern . Im Idealfall liegt die Lösung in der Nähe dieses Tutorials.

Ich kann die Bands mit diesem Befehl zusammenführen

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

Aber aus irgendeinem Grund kann ich RGB-Farben nicht mit dem Befehl imagemagic reparieren. Die Ausgabe ist ~ 700 MB Schwarzbild.

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

Schließlich möchte ich eine Geotiff-Datei haben, um sie auf Mapbox hochzuladen. Eine Erklärung, wie man convertParameter auswählen sollte, ist willkommen.

Ich entwickle eine Anwendung, die erraten sollte, welcher Teil des Satellitenbildes landwirtschaftliche Flächen sind. Ein Szenenbild wird in kleinere Bereiche (möglicherweise 64 x 64) geschnitten und nach CNN ( Zuschneiden oder Nicht-Zuschneiden ) klassifiziert . Ich verwende diesen Datensatz, um das Inception-v3-Modell zu trainieren. Der Datensatz enthält 64 x 64 RGB-Bilder mit einer räumlichen Auflösung von 10 m.

Weitere Infos zu fusioned.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

Dies wird vor und nach dem Anwenden der @ ben-Lösung zusammengeführt

Es gibt 2 Teile des Problems. Das erste ist, dass Sie von 16 Bit in 8 Bit konvertieren möchten, und die Option -scale von gdal_translate tut dies, wie in der vorherigen Antwort erwähnt.

 -scale minOriginal maxOriginal minOutput maxOutput  

Das zweite Problem ist ein Problem der Kontrastverbesserung: Wenn Sie neu skalieren, möchten Sie einen hohen Kontrast für die Pixel haben, an denen Sie interessiert sind. WARNUNG: Es gibt keinen "magischen" Kontrast, da Sie beim erneuten Skalieren normalerweise einige Informationen verlieren : Dies geschieht, um die Visualisierung der Daten zu verbessern, und professionelle Software erledigt dies im laufenden Betrieb, ohne eine neue Datei zu schreiben. Wenn Sie Ihre Daten weiterverarbeiten möchten, enthält Ihr "schwarzer" Geotiff dieselben Informationen wie Ihr jp2 und kann verarbeitet werden. Wenn Sie z. B. den Vegetationsindex berechnen, sollte dies mit den "ursprünglichen" Reflexionswerten erfolgen, nicht mit den neu skalierten. Abgesehen davon sind hier einige Schritte, um ein visuell verbessertes 8-Bit-Bild zu erstellen.

@ben gab Ihnen eine generische Methode, um das Reflexionsvermögen von 0-1 (multipliziert mit 10000 mit diesem Produkt) auf 0-255 neu zu skalieren. Dies ist sicher (kein Ausschluss), aber nur Wolken und einige kahle Böden haben ein wirklich hohes Reflexionsvermögen, sodass Sie an Land nicht viel sehen (außer kahlen Böden) und nichts im Wasser. Daher bestehen Kontrastverbesserungen, die üblicherweise auf Bilder angewendet werden, darin, nur eine Teilmenge des gesamten Bereichs zu erfassen. Auf der sicheren Seite können Sie das Wissen nutzen, dass das maximale Reflexionsvermögen von gewöhnlichem Erdoberflächenmaterial normalerweise unter 0,5 / 0,6 liegt (siehe hier)für einige Beispiele). Dies setzt natürlich voraus, dass Ihr Bild atmosphärisch korrigiert wurde (L2A-Bilder). Der Reflexionsgradbereich ist jedoch in jedem Spektralband unterschiedlich und Sie haben nicht immer die hellsten Erdoberflächen in Ihrem interessierenden Bereich. So sieht die "sichere" Methode aus (mit einem maximalen Reflexionsgrad von 0,4, wie der von @RoVo vorgeschlagene 4096)

Andererseits könnte der Kontrast für jedes Band optimiert werden. Sie können diesen Bereich manuell definieren (z. B. wenn Sie an Wasserfarben interessiert sind und den maximal erwarteten Reflexionsgrad von Wasser kennen) oder basierend auf Bildstatistiken. Eine häufig verwendete Methode besteht darin, ungefähr 95% der Werte beizubehalten und den Rest zu "verwerfen" (zu dunkel -> 0 oder zu hell -> 255), was der Definition des Bereichs basierend auf dem Mittelwert +/- 1,96 * ähnlich ist. Standardabweichung. Dies ist natürlich nur eine Annäherung, da es eine Normalverteilung voraussetzt, aber in der Praxis funktioniert es recht gut (außer wenn Sie zu viele Wolken haben oder wenn die Statistiken einige NoData-Werte verwenden).

Nehmen wir als Beispiel Ihre erste Band:

Mittelwert = 320

std = 536

95% -Konfidenzintervall = [-731: 1372]

Aber natürlich ist das Reflexionsvermögen immer größer als Null, daher müssen Sie das Minimum auf 0 setzen.

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

Und wenn Sie eine aktuelle Version von gdal haben, können Sie -scale_ {band #} verwenden (0 255 ist die Standardausgabe, ich wiederhole sie also nicht), damit Sie keine einzelnen Bänder teilen müssen. Außerdem habe ich vrt anstelle von tif als Zwischendatei verwendet (es ist nicht erforderlich, ein vollständiges Bild zu schreiben: ein virtuelles reicht aus)

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

Beachten Sie, dass Ihre Statistiken stark von "Artefakten" wie Clouds und NoData beeinflusst werden. Auf der einen Seite wird die Varianz überschätzt, wenn Sie Extremwerte haben. Auf der anderen Seite wird Ihr Durchschnitt unterschätzt, wenn es eine große Anzahl von "Null" -Werten gibt (wodurch das automatisch kontrastierte Bild wie im Beispiel zu hell wird), und es würde überschätzt, wenn es eine Mehrheit der Wolken gäbe (was die Bild zu dunkel). Zu diesem Zeitpunkt wären die Ergebnisse daher nicht die besten, die Sie erhalten könnten.

Eine automatisierte Lösung Set Hintergrund und Cloud - Werte auf „nodata“ sein würde und Ihre Statistiken berechnen , ohne den NoData (siehe diesen Beitrag für Details über die Berechnung der Statistiken ohne NoData, und dies für ein Beispiel , um Werte größer als 4000 bis NoData als auch ). Für ein einzelnes Bild berechne ich normalerweise die Statistiken für die größtmögliche wolkenfreie Teilmenge. Mit Statistiken aus einer Teilmenge, in der keine "NoData" (oben links in Ihrem Bild) vorhanden sind, erhalten Sie das Endergebnis. Sie können sehen, dass die Reichweite etwa die Hälfte der "sicheren" Reichweite beträgt, was bedeutet, dass Sie doppelt so viel Kontrast haben:

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

Als letzte Bemerkung sieht gdal_constrast_stretch gut aus, aber ich habe es nicht getestet

Sie können einfach die TCI.jp2Datei verwenden, die in den SAFE.zipDateien enthalten ist. Beachten Sie, dass diese Dateien vor Oktober 2016 nicht in S2-Dateien verfügbar sind

Alternativ können Sie die Bänder mit GDAL konvertieren:

# 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 4096ist ein angemessener Wert für Sentinel-2-Szenen und afaik wird auch für die TCI.jp2-Bilder verwendet. Senken Sie den 4096, wenn Sie ein leichteres Ergebnis erhalten möchten.

Wenn Sie nach einer Lösung suchen, wie Sie sie in der Frage verlinkt haben, sollten Sie das im Tutorial zum Herunterladen bereitgestellte Landsat 8-Verarbeitungs-Shell-Skript befolgen und anpassen .

Insbesondere möchten Sie, wie dort ausgeführt, zunächst die einzelnen Bänder neu skalieren, z. B. wie folgt:

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

Beachten Sie, dass das Histogramm Ihrer Bilder darauf hindeutet, dass Sie nur sehr dunkle Oberflächen in Ihrem Bild haben (ist dies der Fall?), Aber normalerweise ist Ihr Sentinel-2-Bild ein Reflexionsgrad der Atmosphäre oder der Oberfläche, bei dem die Werte üblicherweise zwischen 0 liegen und 10000 - es sei denn, höhere Werte sind ebenfalls möglich, z. B. wenn Sie Wolken im Bild haben.

Anschließend können Sie die Bänder zusammenführen und das Erscheinungsbild des Bildes optimieren:

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

Folgendes passiert mit meinem Bild: