Geotiff Raw Sentinel 2 jp2 to RGB

Szukam sposobu na połączenie plików zespołu Sentinel 2 jp2 ( B02, B03, B04 ) i naprawienie kolorów RGB. Wszystko powinno być zrobione za pomocą skryptu bash lub python. Na mój przykład pracuję nad tymi obrazami . Idealnie rozwiązanie będzie zbliżone do tego samouczka.

Jestem w stanie połączyć pasma za pomocą tego polecenia

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

Ale z jakiegoś powodu nie mogę naprawić kolorów RGB za pomocą polecenia imagemagic. Dane wyjściowe to ~ 700 MB czarnego obrazu.

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

W końcu chciałbym mieć plik geotiff, aby załadować go na mapbox. Wyjaśnienie, jak należy wybrać convertparametry, jest mile widziane.

Tworzę aplikację, która powinna odgadnąć, która część zdjęcia satelitarnego to grunty rolne. Obraz sceny zostanie pocięty na mniejsze łatki (może 64 x 64) i będzie klasyfikowany według CNN ( kadrowanie lub brak kadrowania ). Używam tego zestawu danych do trenowania modelu Inception-v3. Zestaw danych zawiera obrazy RGB 64 x 64 o rozdzielczości przestrzennej 10 m.

Więcej informacji o scalone.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

Jest to scalone.tif przed i po zastosowaniu rozwiązania @ ben

Problem obejmuje 2 części. Po pierwsze, chcesz przekonwertować z 16 bitów na 8 bitów, a robi to opcja -scale programu gdal_translate, jak wspomniano w poprzedniej odpowiedzi.

 -scale minOriginal maxOriginal minOutput maxOutput  

Drugi problem to problem poprawy kontrastu: podczas przeskalowywania chcesz mieć wysoki kontrast dla interesujących Cię pikseli. OSTRZEŻENIE: Nie ma „magicznego” kontrastu, ponieważ podczas przeskalowywania zwykle tracisz pewne informacje : robi się to w celu poprawy wizualizacji danych, a profesjonalne oprogramowanie robi to w locie bez pisania nowego pliku. Jeśli chcesz dalej przetwarzać swoje dane, twoja „czarna” geotiff zawiera te same informacje co twój jp2 i jest gotowa do przetworzenia. W przypadku obliczania np. Wskaźnika wegetacji należy to zrobić z „oryginalnymi” wartościami współczynnika odbicia, a nie z przeskalowanymi wartościami. Biorąc to pod uwagę, oto kilka kroków, aby stworzyć wizualnie ulepszony obraz 8-bitowy.

@ben dał ci ogólną metodę przeskalowania współczynnika odbicia od 0-1 (pomnożonego przez 10000 przez ten produkt) do 0-255. Jest to bezpieczne (bez wyjątku), ale tylko chmury i niektóre gołe gleby mają naprawdę wysokie współczynniki odbicia, więc nie widać dużo na lądzie (z wyjątkiem gołych gleb) i niczego w wodzie. Dlatego ulepszenia kontrastu powszechnie stosowane na zdjęciach polegają na przyjęciu tylko podzbioru pełnego zakresu. Po bezpiecznej stronie możesz wykorzystać wiedzę, że maksymalny współczynnik odbicia wspólnego materiału powierzchni Ziemi wynosi zwykle poniżej 0,5 / 0,6 (patrz tutajdla niektórych przykładów). Oczywiście zakłada to, że twój obraz został poprawiony pod względem atmosferycznym (obrazy L2A). Jednak zakres odbicia różni się w każdym paśmie widmowym i nie zawsze masz najjaśniejsze powierzchnie Ziemi w swoim obszarze zainteresowania. Oto jak wygląda „bezpieczna” metoda (z maksymalnym współczynnikiem odbicia wynoszącym 0,4, jak 4096 sugerowany przez @RoVo)

Z drugiej strony kontrast można zoptymalizować dla każdego pasma. Możesz zdefiniować ten zakres ręcznie (np. Interesuje Cię kolor wody i znasz maksymalną oczekiwaną wartość współczynnika odbicia wody) lub na podstawie statystyk obrazu. Powszechnie stosowana metoda polega na utrzymywaniu około 95% wartości i „odrzucaniu” (zbyt ciemnych -> 0 lub zbyt jasnych -> 255) pozostałych, co jest podobne do definiowania zakresu na podstawie wartości średniej +/- 1,96 * odchylenie standardowe. Oczywiście jest to tylko przybliżenie, ponieważ zakłada rozkład normalny, ale działa całkiem dobrze w praktyce (z wyjątkiem sytuacji, gdy masz za dużo chmur lub jeśli statystyki wykorzystują pewne wartości NoData).

Weźmy twój pierwszy zespół jako przykład:

średnia = 320

standardowe = 536

95% przedział ufności = [-731: 1372]

ale oczywiście współczynnik odbicia jest zawsze większy od zera, dlatego należy ustawić minimum na 0.

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

A jeśli masz najnowszą wersję gdal, możesz użyć opcji -scale_ {band #} (0 255 to wyjście domyślne, więc nie powtarzam tego), abyś nie musiał dzielić pojedynczych pasm. Użyłem również vrt zamiast tif jako pliku pośredniego (nie trzeba pisać pełnego obrazu: wystarczy wirtualny)

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

Zwróć uwagę, że na twoje statystyki silnie wpływają „artefakty”, takie jak chmury i NoData. Z jednej strony wariancja jest przeszacowana, gdy masz ekstremalne wartości. Z drugiej strony, Twoja średnia jest zaniżona, gdy występuje duża liczba wartości „zerowych” (co sprawia, że ​​automatycznie kontrastowany obraz jest zbyt jasny, jak na przykładzie) i byłby zawyżony, gdyby była większość chmur (co spowodowałoby, że obraz za ciemny). Na tym etapie wyniki nie byłyby zatem najlepsze, jakie można uzyskać.

Zautomatyzowanym rozwiązaniem byłoby ustawienie wartości tła i chmury na „nodata” i obliczenie statystyk bez NoData (zobacz ten post, aby uzyskać szczegółowe informacje na temat obliczania statystyk bez NoData, a ten przykład, aby ustawić wartości większe niż 4000 na NoData, a także ). W przypadku pojedynczego obrazu zwykle obliczam statystyki dotyczące największego możliwego podzbioru wolnego od chmury. W przypadku statystyk z podzbioru, w którym nie ma „NoData” (w lewym górnym rogu obrazu), daje to końcowy wynik. Widać, że zasięg stanowi mniej więcej połowę „bezpiecznego” zakresu, co oznacza, że ​​masz dwa razy większy kontrast:

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

Na koniec uwaga: gdal_constrast_stretch wygląda dobrze, ale nie testowałem

Możesz po prostu użyć TCI.jp2pliku zawartego w SAFE.zipplikach. Pamiętaj, że te pliki nie są dostępne w plikach S2 przed październikiem 2016 r

Alternatywnie możesz przekonwertować pasma za pomocą 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 4096jest rozsądną wartością dla scen Sentinel-2 i afaik jest również używany dla obrazów TCI.jp2. Obniż 4096, jeśli chcesz otrzymać jaśniejszy wynik.

Jeśli szukasz rozwiązania takiego jak to, które podłączyłeś w pytaniu, powinieneś postępować i dostosować skrypt powłoki Landsat 8, który jest dostępny do pobrania w samouczku.

W szczególności, jak to się tam dzieje, najpierw możesz przeskalować pojedyncze pasma, np. W następujący sposób:

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

Zauważ, że histogram twoich zdjęć sugeruje, że masz na sobie tylko bardzo ciemne powierzchnie (czy tak jest w przypadku?), Ale zwykle twój obraz wartownika-2 będzie odbijany od góry atmosfery lub na powierzchni, gdzie wartości zwykle mieszczą się w zakresie od 0 i 10000 - chyba że możliwe są również wyższe wartości, np. jeśli masz chmury na obrazie.

Następnie możesz połączyć pasma i dostroić wygląd obrazu:

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

Oto, co dzieje się z moim obrazem: