6

I'm using some Sentinel-2 satellite images in python. Now I have no issues using the newer ones (past 2016). But I need to use some from 2016. These are not preprocessed in the same way by the European Space Agency!

Normally when you download a tile, you usually get a .jp2 file for each of the satellites bands. But in the newer version, they preprocess a RGB version for you alongside the normal bands. This version works awesome for me in python. However, in order to create an RGB version of the older images, I need to combine the three bands (4,3,2 - R,G,B) into 1 file. Gdal_merge handles this quite well at first. When I open the image it looks great! But upon reading it into python, I immediately notice something is off. The image shows up as a purely white picture with some blue streaks across it. Now I went to gdalinfo both the newer working example and the older version I stitched together myself, and this is the output.

enter image description here

As you may notice, the dimensions seems fine at first. But the bands are not of the right type, nor of the right color. So I'm clearly doing something wrong when I merge the files.

This here is the command I use to merge 3 bands into 1 .jp2 file.

gdal_merge.py -o outfile.jp2 -separate B04.jp2 BO3.jp2 BO2.jp2

Now as mentioned. This creates a file, and the file looks beautiful when i open it in QGIS. But its useless to me in python.

Here is a screendump of the python import.

img is the ESA preprocessed image.

img1 is my bastardized gdal_merge import. enter image description here

And here is a image of the failure :p enter image description here

Now it seems to me that I'm lacking some basic understanding of this sort of image manipulation. So with the combined wisdom of Stackoverflow - what can I do to correctly stitch my bands together into a sexy RGB, that can be properly read by the rasterio module.

Thanks in advance :)

Mars
  • 341
  • 1
  • 3
  • 12

1 Answers1

12

What seems to be happening is that you're creating an uint8 stack out of the original uint16 data, so all your values are basically turning into 255, the maximum of uint8.

To fix this, just add -ot uint16 to your call, and everything should be working.

Regarding combining bands, the way I usually go about this is using gdalbuildvrt which creates a virtual dataset from your input files. This .vrt file is only a couple of kb in size, and can be subsequently used for any further GDAL processing (and could potentially read by rasterio):

gdalbuildvrt -separate stack.vrt B04.jp2 BO3.jp2 BO2.jp2

There's many other options you can specify, such as a common resolution (say you want to stack the 10 and 20 meter bands), nodata, target extent, etc.

If you want an actual GeoTIFF, just run it through gdal_translate:

gdal_translate stack.vrt stack.tif

Again, gdal_translate has a ton of cool options, just have a look at the documentation.

Val
  • 6,585
  • 5
  • 22
  • 52
  • Thank you for explaining it to me! I had suspected the -ot had something to do it. But when I tried it earlier, it gave me an error. But shouldnt it be -ot uint8? I want the image to resemble the newer pictures which are done in uint8 – Mars Sep 16 '18 at 10:08
  • I used byte and it did it "right" now. However the color conversion is completely wrong :( – Mars Sep 16 '18 at 11:13
  • 1
    You used byte? The data is uint16, so you would need to scale the data before using byte – Val Sep 16 '18 at 11:38
  • Yeah the byte format is the format of the S2 composite. I want to mimick that. How would I go about scaling the values down? – Mars Sep 17 '18 at 10:28
  • The basic idea is that you get your values in a 0-1 range and then multiply it with the scale factor (255) in this case ... there's ample info on google, try looking for "scale to byte", "scale uint16 to uint8", ... etc. I'm sure there's questions here on SO on that, if not you can always ask yourself ;) – Val Sep 17 '18 at 11:15
  • 1
    I found a way to finally get it right. https://gis.stackexchange.com/questions/257836/raw-sentinel-2-jp2-to-rgb-geotiff/296204#296204 # 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 – Mars Sep 17 '18 at 16:36
  • Thanks. In my case `gdal_merge.py` was crashing with MemoryError and your solution worked without any special memory requirements. – Martin Ždila Mar 09 '19 at 21:47