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.
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
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
#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 ) 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")
You can download the input data and scripts for this tutorial here.