Raw Sentinel 2 jp2 ke geotiff RGB

Saya sedang mencari cara untuk menggabungkan file band Sentinel 2 jp2 ( B02, B03, B04 ) dan memperbaiki warna RGB. Semua harus dilakukan dengan skrip bash atau python. Sebagai contoh saya, saya bekerja pada gambar-gambar ini . Idealnya solusinya akan dekat dengan tutorial ini .

Saya dapat menggabungkan band dengan perintah ini

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

Tapi untuk beberapa alasan saya tidak bisa memperbaiki warna RGB dengan perintah imagemagic. Outputnya adalah ~ 700MB gambar hitam.

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

Akhirnya saya ingin memiliki file geotiff untuk mengunggahnya di mapbox. Penjelasan tentang bagaimana seseorang harus memilih convertparameter dipersilahkan.

Saya sedang mengembangkan aplikasi yang harus menebak bagian mana dari citra satelit adalah tanah pertanian. Gambar pemandangan akan dipotong menjadi tambalan yang lebih kecil (mungkin 64x64) dan akan diklasifikasikan oleh CNN ( tanaman atau bukan tanaman ). Saya menggunakan dataset ini untuk melatih model Inception-v3. Dataset berisi gambar 64x64 RGB dengan resolusi spasial 10m.

Info lebih lanjut tentang 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

Ini digabung.tif sebelum dan sesudah menerapkan solusi @ ben

Ada 2 bagian dari masalah. Yang pertama adalah Anda ingin mengonversi dari 16 bit ke 8bit, dan opsi -scale dari gdal_translate melakukannya, seperti yang disebutkan dalam jawaban sebelumnya.

 -scale minOriginal maxOriginal minOutput maxOutput  

Masalah kedua adalah masalah peningkatan kontras: ketika Anda mengubah skala, Anda ingin memiliki kontras yang tinggi untuk piksel yang Anda minati. PERINGATAN: Tidak ada kontras "ajaib" karena, ketika Anda mengubah skala, Anda biasanya kehilangan beberapa informasi : itu dilakukan untuk meningkatkan visualisasi data, dan perangkat lunak profesional melakukan ini dengan cepat tanpa menulis file baru. Jika Anda ingin memproses lebih lanjut data Anda, geotiff "hitam" Anda berisi informasi yang sama dengan jp2 Anda dan siap untuk diproses. Jika Anda menghitung, mis. Indice vegetasi, ini harus dilakukan dengan nilai reflektansi "asli", bukan yang direccaled. Yang sedang berkata, berikut adalah beberapa langkah untuk membuat gambar 8bit yang ditingkatkan secara visual.

@ben memberi Anda metode umum untuk mengubah skala pemantulan dari 0-1 (dikalikan dengan 10.000 dengan produk ini) ke 0-255. Ini aman (tanpa pengecualian), tetapi hanya awan dan beberapa tanah gundul yang memiliki pantulan sangat tinggi, sehingga Anda tidak melihat banyak hal di darat (kecuali tanah gundul) dan tidak ada apa pun di air. Oleh karena itu, peningkatan kontras yang biasa diterapkan pada gambar terdiri dari pengambilan hanya sebagian dari rentang penuh. Di sisi yang aman, Anda dapat menggunakan pengetahuan bahwa maksimum pemantulan material permukaan bumi pada umumnya di bawah 0,5 / 0,6 (lihat di siniuntuk beberapa contoh). Tentu saja, ini mengasumsikan bahwa gambar Anda telah diperbaiki secara atmosfer (gambar L2A). Namun, rentang pantulan berbeda di setiap pita spektral dan Anda tidak selalu memiliki permukaan Bumi paling terang di bidang yang Anda minati. Berikut adalah bagaimana metode "aman" terlihat (dengan pemantulan maksimal 0,4, seperti 4096 yang disarankan oleh @RoVo)

Di sisi lain, kontras dapat dioptimalkan untuk setiap band. Anda dapat menentukan kisaran ini secara manual (mis. Anda tertarik dengan warna air dan Anda tahu nilai reflektansi maksimum yang diharapkan dari air) atau berdasarkan statistik gambar. Metode yang umum digunakan terdiri dari menjaga sekitar 95% dari nilai dan "membuang" (terlalu gelap -> 0 atau terlalu terang -> 255) sisanya, yang mirip dengan mendefinisikan rentang berdasarkan nilai rata-rata +/- 1.96 * standar deviasi. Tentu saja, ini hanya perkiraan karena mengasumsikan distribusi normal, tetapi berfungsi cukup baik dalam praktiknya (kecuali ketika Anda memiliki terlalu banyak cloud atau jika statistik menggunakan beberapa nilai NoData).

Mari ambil band pertama Anda sebagai contoh:

berarti = 320

std = 536

Interval kepercayaan 95% = [-731: 1372]

tetapi tentu saja pantulan selalu lebih besar dari nol, oleh karena itu Anda harus menetapkan minimum pada 0.

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

Dan jika Anda memiliki versi gdal terbaru, maka Anda dapat menggunakan -scale_ {band #} (0 255 adalah output default, jadi saya tidak mengulanginya) sehingga Anda tidak perlu membagi band tunggal. Saya juga menggunakan vrt bukan tif sebagai file perantara (tidak perlu menulis gambar penuh: yang virtual sudah cukup)

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

Perhatikan bahwa statistik Anda sangat dipengaruhi oleh "artefak" seperti cloud dan NoData. Di satu sisi, varians terlalu tinggi ketika Anda memiliki nilai ekstrim. Di sisi lain, rata-rata Anda diremehkan ketika ada sejumlah besar nilai "nol" (membuat gambar yang dikontraskan secara otomatis terlalu terang seperti pada contoh) dan itu akan terlalu berlebih jika ada mayoritas awan (yang akan membuat gambar terlalu gelap). Pada tahap ini, hasilnya bukan yang terbaik yang bisa Anda dapatkan.

Solusi otomatis adalah mengatur nilai latar belakang dan cloud menjadi "nodata" dan menghitung statistik Anda tanpa NoData (lihat posting ini untuk detail tentang statistik penghitungan tanpa NoData, dan ini sebagai contoh untuk menetapkan nilai yang lebih besar dari 4000 ke NoData juga) ). Untuk satu gambar, saya biasanya menghitung statistik pada subset bebas awan terbesar yang mungkin. Dengan statistik dari subset di mana tidak ada "NoData" (kiri atas gambar Anda), ini memberikan hasil akhir. Anda dapat melihat bahwa kisarannya sekitar setengah dari kisaran "aman", yang berarti Anda memiliki kontras dua kali lebih banyak:

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

Sebagai komentar terakhir, gdal_constrast_stretch terlihat bagus tetapi saya belum menguji

Anda cukup menggunakan TCI.jp2file yang disertakan dalam SAFE.zipfile. Perhatikan bahwa file-file ini tidak tersedia dalam file S2 sebelum Oktober 2016

Atau Anda dapat mengonversi band menggunakan 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 4096adalah nilai yang wajar untuk adegan Sentinel-2 dan afaik juga digunakan untuk gambar TCI.jp2. Turunkan 4096 jika Anda ingin menerima hasil yang lebih ringan.

Jika Anda mencari solusi seperti yang Anda tautkan dalam pertanyaan, Anda harus mengikuti dan menyesuaikan skrip shell pemrosesan Landsat 8 yang disediakan untuk diunduh dalam tutorial.

Khususnya, seperti yang dilakukan di sana, Anda mungkin ingin mengubah skala band tunggal, misalnya sebagai berikut:

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

Perhatikan bahwa histogram gambar Anda menunjukkan bahwa Anda hanya memiliki permukaan yang sangat gelap pada gambar Anda (apakah ini masalahnya?) Tetapi biasanya gambar sentinel-2 Anda akan menjadi atmosfer-atas atau reflektansi permukaan yang nilai umumnya berkisar antara 0 dan 10000 - kecuali nilai yang lebih tinggi juga dimungkinkan misalnya jika Anda memiliki awan di gambar.

Kemudian Anda dapat menggabungkan band dan menyempurnakan tampilan gambar:

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

Inilah yang terjadi pada gambar saya ketika melakukan ini: