Raw Sentinel 2 jp2 в RGB геотиф

Я ищу способ объединить файлы полос Sentinel 2 jp2 ( B02, B03, B04 ) и исправить цвета RGB. Все должно быть сделано с помощью скриптов bash или python. Для моего примера я работаю над этими изображениями . В идеале решение будет близко к этому уроку.

Я могу объединить группы с помощью этой команды

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

Но по какой-то причине я не могу исправить цвета RGB с помощью команды imagemagic. На выходе получается ~ 700 МБ черного изображения.

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

В конце концов я хотел бы иметь файл geotiff, чтобы загрузить его на mapbox. Объяснение того, как следует выбирать convertпараметры, приветствуется.

Я разрабатываю приложение, которое должно угадать, какая часть спутникового изображения является сельскохозяйственной землей. Изображение сцены будет разрезано на более мелкие участки (возможно, 64x64) и будет классифицировано по CNN ( кадрирование или не кадрирование ). Я использую этот набор данных для обучения модели Inception-v3. Набор данных содержит изображения 64x64 RGB с пространственным разрешением 10 м.

Больше информации о 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

Это merged.tif до и после применения решения @ ben

Есть 2 части проблемы. Во-первых, вы хотите преобразовать из 16 бит в 8 бит, и опция -scale в gdal_translate это делает, как упоминалось в предыдущем ответе.

 -scale minOriginal maxOriginal minOutput maxOutput  

Вторая проблема - это проблема повышения контрастности: когда вы изменяете масштаб, вы хотите иметь высокую контрастность для интересующих вас пикселей. ПРЕДУПРЕЖДЕНИЕ. «Волшебный» контраст отсутствует, поскольку при изменении масштаба вы обычно теряете некоторую информацию : это делается для улучшения визуализации данных, и профессиональные программы делают это на лету, не записывая новый файл. Если вы хотите продолжить обработку ваших данных, ваш «черный» геотиф содержит ту же информацию, что и ваш jp2, и готов к обработке. Если вы вычисляете, например, индекс растительности, это следует делать с «исходными» значениями отражательной способности, а не с пересчитанными. При этом, вот несколько шагов для создания визуально улучшенного 8-битного изображения.

@ben дал вам общий метод для изменения масштаба отражения от 0-1 (умноженный на 10000 с этим продуктом) до 0-255. Это безопасно (не исключая), но только облака и некоторые голые почвы имеют очень высокую отражательную способность, поэтому вы не видите много на суше (кроме голых почв) и ничего в воде. Следовательно, контрастные улучшения, обычно применяемые к изображениям, заключаются в том, чтобы брать только подмножество полного диапазона. На всякий случай, вы можете использовать знания о том, что максимальная отражательная способность материала общей поверхности Земли обычно ниже 0,5 / 0,6 (см. Здесьдля некоторых примеров). Конечно, это предполагает, что ваше изображение было атмосферно скорректировано (изображения L2A). Тем не менее, диапазон отражательной способности отличается в каждой спектральной полосе, и у вас не всегда самые яркие поверхности Земли в вашей области интересов. Вот как выглядит «безопасный» метод (с максимальной отражательной способностью 0,4, как 4096, предложенный @RoVo)

С другой стороны, контраст можно оптимизировать для каждой полосы. Вы можете определить этот диапазон вручную (например, вас интересует цвет акварели и вы знаете максимальное ожидаемое значение отражения воды) или на основе статистики изображения. Обычно используемый метод заключается в сохранении приблизительно 95% значений и «отбрасывании» (слишком темный -> 0 или слишком яркий -> 255) остальных, что аналогично определению диапазона на основе среднего значения +/- 1,96 * среднеквадратичное отклонение. Конечно, это только приблизительное значение, поскольку оно предполагает нормальное распределение, но на практике работает довольно хорошо (за исключением случаев, когда у вас слишком много облаков или если статистика использует некоторые значения NoData).

Давайте возьмем вашу первую группу в качестве примера:

среднее = 320

стандартный = 536

95% доверительный интервал = [-731: 1372]

но, конечно, коэффициент отражения всегда больше нуля, поэтому вы должны установить минимум на 0.

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

А если у вас последняя версия gdal, вы можете использовать -scale_ {band #} (0 255 - выход по умолчанию, поэтому я не повторяю его), чтобы вам не приходилось разделять отдельные полосы. Также я использовал vrt вместо tif в качестве промежуточного файла (не нужно писать полное изображение: достаточно виртуального)

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

Обратите внимание, что на вашу статистику сильно влияют "артефакты", такие как облака и NoData. С одной стороны, дисперсия завышена, когда у вас есть экстремальные значения. С другой стороны, ваше среднее значение занижено, когда имеется большое количество «нулевых» значений (что делает автоматически контрастное изображение слишком ярким, как в примере), и оно будет завышено, если будет большинство облаков (что приведет к изображение слишком темное). Поэтому на данном этапе результаты не будут лучшими, которые вы могли бы получить.

Автоматизированным решением было бы установить фоновые и облачные значения для «nodata» и вычислить вашу статистику без NoData (см. Этот пост для получения подробной информации о вычислении статистики без NoData, и этот пример для примера, чтобы установить значения больше 4000 для NoData также ). Для одного изображения я обычно вычисляю статистику для максимально возможного безоблачного подмножества. Со статистикой из подмножества, где нет «NoData» (вверху слева от вашего изображения), это даст окончательный результат. Вы можете видеть, что диапазон составляет примерно половину «безопасного» диапазона, что означает, что у вас вдвое больше контраста:

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

Как последнее замечание, gdal_constrast_stretch выглядит хорошо, но я не проверял

Вы можете просто использовать TCI.jp2файл, который включен в SAFE.zipфайлы. Обратите внимание, что эти файлы недоступны в файлах S2 до октября 2016 года.

В качестве альтернативы вы можете конвертировать полосы, используя 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 4096является разумным значением для сцен Sentinel-2, а afaik также используется для изображений TCI.jp2. Опустите 4096, если вы хотите получить более легкий результат.

Если вы ищете решение, которое вы связали в вопросе, вам следует следовать и настроить сценарий обработки оболочки Landsat 8 , предоставленный для загрузки в руководстве.

В частности, как это делается, сначала вы можете захотеть изменить масштаб отдельных полос, например, следующим образом:

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

Обратите внимание, что гистограмма ваших изображений предполагает, что у вас на изображении только очень темные поверхности (это так?), Но обычно ваше изображение часового-2 будет отражать верхнюю часть атмосферы или отраженную поверхность, где значения обычно находятся в диапазоне от 0 и 10000 - если также не возможны более высокие значения, например, если у вас есть облака на изображении.

Затем вы можете объединить полосы и настроить внешний вид изображения:

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

Вот что происходит с моим изображением при этом: