I have downloaded several Sentinel 2 granules, with single bands raster data in .jp2 format, and I want to create true and false color composite images out of them. Using ArcGIS I could not open any single band raster (error says "invalid dataset"). Using QGIS I can open them but I cannot composite (I have tried using both "merge" and "create virtual raster" tool). I also tried to convert into geotiff using python scripts, but I am not that good in python, thus it did not work and probably I would need step by step guidance to do that. Does anyone have some ideas to solve the problem?
[GIS] Handle .JP2 Sentinel data
jpeg 2000sentinel-2
Related Solutions
It looks like this is a problem with the conda packaging of gdal. I see that the current gdal version on conda is 2.1.3, not 2.1.0, so hopefully this packaged version works. Consider removing the current gdal and reinstalling.
I looked at the contents of the gdal conda package I mentioned and for 64-bit linux python 2.7 numpy 1.12 the files gdaljp2abstractdataset.h, gdaljp2metadata.h, and gdaljp2metadatagenerator.h are included. Haven't installed it to confirm whether the driver exists.
There are 2 parts of the problem. The first is that you want to convert from 16 bits to 8bit, and the -scale option of gdal_translate does it, as mentioned in the previous answer.
-scale minOriginal maxOriginal minOutput maxOutput
The second problem is a contrast enhancement problem: when you rescale, you want to have a high contrast for the pixels that you are interested in. WARNING : There is no "magic" contrast because, when you rescale, you usually loose some information: it is done to improve the visualization of the data, and professionnal softwares do this on the fly without writing a new file. If you want to further process your data, your "black" geotiff contain the same information as your jp2 and is ready to be processed. If you compute, e.g., vegetation indice, this should be done with the "original" reflectance values, not the rescaled ones. That being said, here are some steps to create a visually enhanced 8bit image.
@ben gave you a generic method to rescale the reflectance from 0-1 (multiplied by 10000 with this product) to 0-255. This is safe (no exclusion), but only clouds and some bare soils have really high reflectances, so you don't see much on land (except bare soils) and nothing in water. Therefore, contrast enhancements commonly applied on images consist in taking only a subset of the full range. On the safe side, you can use the knowledge that the max reflectance of common Earth surface material is usually below 0.5/0.6 (see here for some examples). Of course, this assumes that your image has been atmospherically corrected (L2A images). However, the range of reflectance differs in each spectral band and you don't always have the brightest Earth surfaces in your area of interest. Here is how the "safe" method looks like (with a max reflectance of 0.4, like the 4096 suggested by @RoVo)
On the other hand, contrast could be optimized for each band. You can define this range manually (e.g. you are interested in water colour and you know the maximum expected reflectance value of water) or based on image statistics. A commonly used method consists in keeping approximately 95% of the values and "discarding" (too dark -> 0 or too bright -> 255) the rest, which is similar to defining the range based on the mean value +/- 1.96*standard deviation. Of course, this is only an approximation because it assumes a normal distribution, but it works quite good in practice (except when you have too many clouds or if the stats make use of some NoData values).
Lets take your first band as example :
mean = 320
std = 536
95% confidence interval = [-731:1372]
but of course the reflectance is always larger than zero, therefore you must set the minimum at 0.
gdal_translate -scale 0 1372 0 255 -ot Byte B01.jp2 B01-scaled.tif
And if you have a recent version of gdal, then you can use -scale_{band#} (0 255 is the default output, so I don't repeat it) so that you don't need to split single bands. Also I used vrt instead of tif as an intermediate file (no need to write a full image: a virtual one is enough)
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
Note that your statistics are strongly affected by "artifacts" such as clouds and NoData. On one side, the variance is overestimated when you have extreme values. On the other side, your average is understimated when there is a large amount of "zero" values (making the automatically contrasted image too bright like on the example) and it would be overstimated if there was a majority of clouds (which would make the image too dark). At this stage, the results would therefore not be the best you could get.
An automated solution would be to set background and cloud values to "nodata" and compute your stats without the NoData (see this post for details on computing stats without NoData, and this one for an example to set values larger than 4000 to NoData as well). For a single image, I usually compute the statistics on the largest possible cloud-free subset. With stats from a subset where there are no "NoData" (top left of your image), this give the final result. You can see that the range is about half of the "safe" range, which means that you have twice as much contrast:
gdal_translate -scale_1 38 2225 -scale_2 553 1858 -scale_3 714 1745 -ot Byte stack.vrt im_rescaled.tif
As a last remark, gdal_constrast_stretch looks good but I haven't tested
Best Answer
The latest version of GDAL (2.1) has the ability to read Sentinel 2 data (see http://gdal.org/frmt_sentinel2.html) in the SAFE file structure. It reads as sub-datasets so you need to run
gdalinfofirst on the main XML file to show all subdatasets (one for each resolution and UTM zone) then to select the one you want. You can then create a GeoTiff from the selected dataset using the following command:(Example from GDAL website). The GeoTiff will be a mosaic of all .jp2 tiles within the given UTM zone comprising all bands for the selected resolution and can be loaded into QGIS/ArcMap.
Rather than using
gdalinfoyou can get a list of all the subdatasets using the GDAL Python bindings:Then pass each to
gdal_translate. I used a similar approach in a script I wrote (see https://spectraldifferences.wordpress.com/2016/08/16/convert-sentinel-2-data-using-gdal/ for more details.)You can install the latest version of GDAL under Windows (or Linux and OS X) using conda (http://conda.pydata.org/miniconda.html) through the conda-forge channel. Once miniconda has been installed run the following steps.