Sentinel 2 is a multispectral satellite with 13 bands of spatial resolutions from 10m to 60m launched by the European Spatial Agency (ESA). On the visible spectrum and near infrared (bands 2, 3, 4 and 8) the spatial resolution is 10m with a global coverage every five days. The mission is a constellation of two satellites: Sentinel-2A and Sentinel-2B with images from June 2015 and March 2017 respectively. The high resolution, temporal distribution and availability make Sentinel 2 images a valuable and free resource for agricultural monitoring, land cover classification or water quality. There are different ways to download Sentinel 2 images, our recommended option to explore, previsualize and download images is the Copernicus Open Access Hub: https://scihub.copernicus.eu/dhus/#/home

Rasterio is a Python library that allows to read, inspect, visualize and write geospatial raster data. The library uses GeoTIFF and other spatial raster formats and is capable of working with satellite imagery, digital elevation models, and drone imagery data products. Rasterio allows you to import a single band or multiband geospatial raster in a interactive Python enviroment as Jupyter notebook, the library can keep the “duality” of the geospatial raster, that means, it can handle the location and resolution parameters as well as the matrix values of the gridded elements.

This tutorial shows some basic procedures to explore a multiband Sentinel 2 granule with Python 3 and Rasterio on a Jupyter Notebook. The tutorial shows the commands to identify the raster array dimensions and the geospatial referencing parameters, make representation of each visible band and export band composites as true color and false color geoespatial rasters in Tiff format.

Geospatial rasters in Python, a long history with a happy end

From our experience, geospatial raster manipulation and analysis was a pretty hard task on a enviroment different from a GIS desktop program as QGIS. On a mathematical perspective a geospatial raster is a matrix array, were every cell is georeferenced to a extension on surface. As a GIS user, we treat geospatial rasters as images, that means, something that can be seen an interpreted by our eyes; however as a numerical model, a geospatial raster is a matrix array where values represent something from the field.

There was a missing link between the matrix array and the geoespatial attributes. It was not possible to have full control of the matrix values with the geospatial location. However, nowadays the georeferentiation and manipulation of raster arrays is pretty simple with Rasterio.

Multiple band (red, green, blue) representation with Rasterio and Matplotlib

Histogram from a Sentinel 2 image done with Rasterio plotting utilities

Representation of a processed false color image with Rasterio in QGIS 3

Useful links

In order to work with Python in Windows users we recommend to work with Anaconda that you can download free from: https://www.anaconda.com/download/

Rasterio and other geospatial libraries can be a little bit tricky to install, this is a tutorial to install these libraries without problems: https://youtu.be/4ybddFC80fU

More information about Rasterio: https://github.com/mapbox/rasterio

Tutorial

Python code

This is the python code used on this tutorial:

#import required libraries
import rasterio
from rasterio import plot
import matplotlib.pyplot as plt
%matplotlib inline

#import bands as separate 1 band raster
imagePath = '../Sentinel2/GRANULE/L1C_T11SKB_A007675_20180825T184430/IMG_DATA/'
band2 = rasterio.open(imagePath+'T11SKB_20180825T183909_B02.jp2', driver='JP2OpenJPEG') #blue
band3 = rasterio.open(imagePath+'T11SKB_20180825T183909_B03.jp2', driver='JP2OpenJPEG') #green
band4 = rasterio.open(imagePath+'T11SKB_20180825T183909_B04.jp2', driver='JP2OpenJPEG') #red
band8 = rasterio.open(imagePath+'T11SKB_20180825T183909_B08.jp2', driver='JP2OpenJPEG') #nir

#number of raster bands
band4.count

#number of raster columns
band4.width

#number of raster rows
band4.height

#plot band 
plot.show(band4)

#type of raster byte
band4.dtypes[0]

#raster sytem of reference
band4.crs

#raster transform parameters
band4.transform

#raster values as matrix array
band4.read(1)

#multiple band representation
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(12, 4))
plot.show(band2, ax=ax1, cmap='Blues')
plot.show(band3, ax=ax2, cmap='Greens')
plot.show(band4, ax=ax3, cmap='Reds')
fig.tight_layout()

#export true color image
trueColor = rasterio.open('../Output/SentinelTrueColor2.tiff','w',driver='Gtiff',
                         width=band4.width, height=band4.height,
                         count=3,
                         crs=band4.crs,
                         transform=band4.transform,
                         dtype=band4.dtypes[0]
                         )
trueColor.write(band2.read(1),3) #blue
trueColor.write(band3.read(1),2) #green
trueColor.write(band4.read(1),1) #red
trueColor.close()

src = rasterio.open(r"../Output/SentinelTrueColor2.tiff", count=3)
plot.show(src)

#export false color image
falseColor = rasterio.open('../Output/SentinelFalseColor.tiff', 'w', driver='Gtiff',
                          width=band2.width, height=band2.height,
                          count=3,
                          crs=band2.crs,
                          transform=band2.transform,
                          dtype='uint16'                   
                         )
falseColor.write(band3.read(1),3) #Blue
falseColor.write(band4.read(1),2) #Green
falseColor.write(band8.read(1),1) #Red
falseColor.close()

#generate histogram
trueColor = rasterio.open('../Output/SentinelTrueColor2.tiff')
plot.show_hist(trueColor, bins=50, lw=0.0, stacked=False, alpha=0.3, histtype='stepfilled', title="Histogram")