Raw Sentinel 2 jp2 zu RGB Geotiff

11

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 Vor nach

gkiko
quelle
1
Wie hoch ist die Bittiefe von Merged.tif und Min, Average und Max im Histogramm? Überprüfen Sie mitgdalinfo -hist merged.tif
user30184
@ user30184 Ich habe die angeforderten Informationen zu meiner Frage hinzugefügt
gkiko
Ich habe versucht, jp2 in geotiff zu konvertieren und dann die Farbkorrektur anzuwenden, aber ich bekomme immer noch ein schwarzes Bild
gkiko
Warum verwenden Sie nicht einfach das Bild TCI.jp2, das im Grunde dasselbe ist wie -scale 0 4096 0 255?
pLumo
1
Für Ernte / Nicht-Ernte können Sie vielleicht diese esa-sen2agri.org/resources/software verwenden, anstatt Ihre eigene Anwendung von Grund auf neu zu
erstellen

Antworten:

8

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)

Geben Sie hier die Bildbeschreibung ein

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.

Geben Sie hier die Bildbeschreibung ein

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

Geben Sie hier die Bildbeschreibung ein

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

Radouxju
quelle
Das Problem dabei ist, dass jedes Granulat eine andere Helligkeit hat. Je nachdem, was er erreichen möchte, ist es besser, eine feste Skala zu verwenden. -scale 0 4096 0 255erzeugt eine ziemlich gute Ausgabe, wenn wir keine Wolkentexturen benötigen ...
pLumo
@RoVo Ich bin damit einverstanden, dass dies zu Helligkeitswerten führt und dass Sie auf hellen Oberflächen wie Sand möglicherweise den Kontrast verlieren. Dies basiert jedoch auf den Statistiken des vom OP zusammengeführten Bildes. Sie werden keinen unterschiedlichen Kontrast auf dem Granulat haben. Normalerweise ist der Bereich in Rot, Grün und Blau viel kleiner als der Bereich in NIR. Aus diesem Grund ist es sinnvoll, für jedes Band einen anderen Kontrast zu verwenden.
Radouxju
7

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.

pLumo
quelle
5

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:

Geben Sie hier die Bildbeschreibung ein

ben
quelle
1
Ich habe meine Frage aktualisiert. Wie soll ich entscheiden, welche Parameter bei der Farbkorrektur von geoTIFF verwendet werden sollen?
Gkiko
Achten Sie beim Skalieren von Werten vom Eingabe- zum Ausgabebild immer auf den Max- und Min-Wert im Eingabebild. Zum Beispiel sollte der Skalierungsparameter für das erste Band wie folgt lauten: -scale 0 4818 0 255.
Milos Miletic