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 gdalinfo
first 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:
gdal_translate SENTINEL2_L1C:S2A_OPER_MTD_SAFL1C_PDMC_20150818T101440_R022_V20150813T102406_20150813T102406.xml:10m:EPSG_32632 \
output_10m.tif
(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 gdalinfo
you can get a list of all the subdatasets using the GDAL Python bindings:
from osgeo import gdal
dataset = gdal.Open('S2/S2.xml', gdal.GA_ReadOnly)
subdatasets = dataset.GetSubDatasets()
dataset = None
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.
conda create -n gdal2 -c conda-forge gdal
source activate gdal2
I can't reproduce your issue with GDAL 2.2-dev (yesterday's trunk version from gisinternals.com).
I downloaded a rando Sentinel2 zip file from https://scihub.copernicus.eu/dhus/#/home. It happened to be "S2A_MSIL1C_20170129T110311_N0204_R094_T31UCT_20170129T110306.zip".
I unzipped the file and changed my working directory into
C:\data\S2A_MSIL1C_20170129T110311_N0204_R094_T31UCT_20170129T110306.SA
FE\GRANULE\L1C_T31UCT_A008385_20170129T110306\IMG_DATA>
I used this command:
C:\data\sentinel\S2A_MSIL1C_20170129T110311_N0204_R094_T31UCT_20170129T110306.SA
FE\GRANULE\L1C_T31UCT_A008385_20170129T110306\IMG_DATA>gdal_translate T31UCT_201
70129T110311_B08.jp2 c:\data\sentinel.tif
Input file size is 10980, 10980
0...10...20...30...40...50...60...70...80...90...100 - done.
Gdalinfo report from the resulting tiff file:
C:\data>gdalinfo sentinel.tif
Driver: GTiff/GeoTIFF
Files: sentinel.tif
Size is 10980, 10980
Coordinate System is:
PROJCS["WGS 84 / UTM zone 31N",
GEOGCS["WGS 84",
DATUM["WGS_1984",
SPHEROID["WGS 84",6378137,298.257223563,
AUTHORITY["EPSG","7030"]],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0,
AUTHORITY["EPSG","8901"]],
UNIT["degree",0.0174532925199433,
AUTHORITY["EPSG","9122"]],
AUTHORITY["EPSG","4326"]],
PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",3],
PARAMETER["scale_factor",0.9996],
PARAMETER["false_easting",500000],
PARAMETER["false_northing",0],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]],
AXIS["Easting",EAST],
AXIS["Northing",NORTH],
AUTHORITY["EPSG","32631"]]
Origin = (300000.000000000000000,5800020.000000000000000)
Pixel Size = (10.000000000000000,-10.000000000000000)
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
INTERLEAVE=BAND
Corner Coordinates:
Upper Left ( 300000.000, 5800020.000) ( 0d 3'56.87"E, 52d18'50.48"N)
Lower Left ( 300000.000, 5690220.000) ( 0d 7'45.16"E, 51d19'40.99"N)
Upper Right ( 409800.000, 5800020.000) ( 1d40'33.33"E, 52d20'35.00"N)
Lower Right ( 409800.000, 5690220.000) ( 1d42'16.48"E, 51d21'21.89"N)
Center ( 354900.000, 5745120.000) ( 0d53'37.96"E, 51d50'16.88"N)
Band 1 Block=10980x1 Type=UInt16, ColorInterp=Gray
Image Structure Metadata:
NBITS=15
Best Answer
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.
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 :
but of course the reflectance is always larger than zero, therefore you must set the minimum at 0.
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)
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:
As a last remark, gdal_constrast_stretch looks good but I haven't tested