Merging multiple bands together through gdal...correctly

Question

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.

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.

And here is a image of the failure :p

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 :)

Answer 1

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.