How do I plot a satellite image using rasterio? I'm able to plot the individual bands but not the entire image. Is there a way to combine the separate R, G, B numpy arrays to create a standard RGB image?


The Rasterio Plotting documentation describes how to visualize multiband imagery. For example, using 4-band NAIP imagery:

import rasterio
from rasterio.plot import show
src = rasterio.open("path/to/your/image/m_3511642_sw_11_1_20140704.tif")


enter image description here

To visualize specific band combination use the following approach (source). In this case, I am creating a false color composite image using the NIR band:

import rasterio
import numpy as np
import matplotlib.pyplot as plt

# Open the file:
raster = rasterio.open('path/to/your/image/m_3511642_sw_11_1_20140704.tif')

# Normalize bands into 0.0 - 1.0 scale
def normalize(array):
    array_min, array_max = array.min(), array.max()
    return (array - array_min) / (array_max - array_min)

# Convert to numpy arrays
nir = raster.read(4)
red = raster.read(3)
green = raster.read(2)

# Normalize band DN
nir_norm = normalize(nir)
red_norm = normalize(red)
green_norm = normalize(green)

# Stack bands
nrg = np.dstack((nir_norm, red_norm, green_norm))

# View the color composite

enter image description here

  • Thanks a lot @Aaron. But following your code snippet resulted in the outputting of the blue channel of the image. The documentation of the rasterio plot function states that if the dataset is of raster format, display the rgb image as defined in the colorinterp metadata, or default to first band. Perhaps this is what happened as the first band is blue and my image might not be of the correct format. So how would I go about plotting the rgb image? Do I have to convert my file to a particular format or is there a function that takes it into account? – Ashwin Nair Dec 17 '18 at 6:42
  • @AshwinNair I’ll update my answer, but in the meantime check this out: automating-gis-processes.github.io/CSC/notebooks/L5/… – Aaron Dec 17 '18 at 19:57

i used earthpy to carry out the layerstacking

import os from glob import glob import matplotlib.pyplot as plt import rasterio as rio from rasterio.plot import plotting_extent import geopandas as gpd import earthpy as et import earthpy.spatial as es import earthpy.plot as ep from pathlib import Path

#fetch the data

landsat=Path(r'C:\Users\AFRO TEOP\programming\georaster\data')


landsat_bands.sort() landsat_bands=landsat_bands[3:6] #place the right array to source rgb bands


array_stack, extent = es.stack(landsat_bands, nodata=-9999) extent = plotting_extent(array_stack[0], extent["transform"])

ep.plot_rgb( array_stack, rgb=(4,3,2), # place the combination of rgb ax=ax, stretch=True, extent=extent, ) plt.show()

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.