Visualizing multi-band satellite images in Python

In this tutorial, I explain how to visualize multi-band Landsat 8 satellite imagery as true- and false-color composites using Python and matplotlib. The basic procedure is: 1) read individual raster bands as array 2) rescale to 0 and 1 using minimum and maximum values (or percentiles or mean + - standard deviation for contrast enhancement) 3) stack arrays using numpy.dstack 4) plot RGB composite with plt.imshow()

Code:

from osgeo import gdal
import numpy as np
import matplotlib.pyplot as plt

def scaleMinMax(x):
return((x - np.nanmin(x))/(np.nanmax(x) - np.nanmin(x)))
def scaleCCC(x):
return((x - np.nanpercentile(x, 2))/(np.nanpercentile(x, 98) - np.nanpercentile(x,2)))
def scaleStd(x):
return((x - (np.nanmean(x)-np.nanstd(x)*2))/((np.nanmean(x)+np.nanstd(x)*2) - (np.nanmean(x)-np.nanstd(x)*2)))

ds = gdal.Open("Landsat8_20200423_B1-7.tif")

r = ds.GetRasterBand(4).ReadAsArray()
g = ds.GetRasterBand(3).ReadAsArray()
b = ds.GetRasterBand(2).ReadAsArray()

ds = None

rMinMax = scaleMinMax(r)
gMinMax = scaleMinMax(g)
bMinMax = scaleMinMax(b)

rgbMinMax = np.dstack((rMinMax,gMinMax,bMinMax))
plt.figure()
plt.imshow(rgbMinMax)
plt.show()

rCCC = scaleCCC(r)
gCCC = scaleCCC(g)
bCCC = scaleCCC(b)

rgbCCC = np.dstack((rCCC,gCCC,bCCC))
plt.figure()
plt.imshow(rgbCCC)
plt.show()

rStd = scaleStd(r)
gStd = scaleStd(g)
bStd = scaleStd(b)

rgbStd = np.dstack((rStd,gStd,bStd))
plt.figure()
plt.imshow(rgbStd)
plt.show()