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()