Handle .JP2 Sentinel data

Question

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?

Answer 1

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

Answer 2

The easiest way is probably to download the sentinel toolbox software, with which you can clip (if necessary), stack and convert the bands to Gtiff. However in my experience this software is quite slow, and can hog system resources while running.

A better way might be to use gdal utilities, if you have it installed on your system? The gdal_translate command line tool can convert the jp2 to a GTiff:

gdal_translate -of "Gtiff" in_file_name.jp2 out_file_name.tif

You can also access the gdal utilities from python:

from osgeo import gdal

in_image = gdal.Open("in_file_name.jp2")

driver = gdal.GetDriverByName("GTiff")

out_image = driver.CreateCopy("out_file_name.jp2", in_image, 0)

in_image = None
out_image = None

If you want to stack and convert all the images in a python script, you can use the Remote Sensing and GIS Library (RSGISLib) for python. This is only available on Unix systems though, so if you have windows you will need to use a virtual machine.

import rsgislib

images = ['sentinel_band_1.jp2', 'sentinel_band_2.jp2' ...]

band_names = ['band_1_name', 'band_2_name'...]

out_image = 'sentinel_image_stack.tif'

data_type = rsgislib.TYPE_16UINT

imageutils.stackImageBands(images, band_names, out_image, None, 0, 'Gtiff', data_type) 

Answer 3

ArcGIS sometimes gives an error message when you attempt to add .jp2 files using certain commands. However, when you just drag and drop the file from a folder in your computer to the table of contents, ArcGIS can accept the .jp2 files.

Locate the .jp2 images in your systems folder then "drag and drop" to the Arcgis Table of Contents (the list of files in your map).