7

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?

1

3 Answers 3

9

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

show(src)

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
plt.imshow(nrg)

enter image description here

2
  • 1
    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?
    – ashnair1
    Dec 17, 2018 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, 2018 at 19:57
3

Adding to the answers already here. This will allow you to plot using show (so the CRS is maintained) and normalize using a percentile clip (e.g. 2nd & 98th) and account for NAN values. The example takes a 23 band Sentinel image, extracts bands 2 (blue), 3 (green), and 2 (red), applies a percentile clip, writes a new rgb image and plots it with show.

def pct_clip(array,pct=[2,98]):
    array_min, array_max = np.nanpercentile(array,pct[0]), np.nanpercentile(array,pct[1])
    clip = (array - array_min) / (array_max - array_min)
    clip[clip>1]=1
    clip[clip<0]=0
    return clip

with rio.open("C:\\Users\\User\\Google Drive\\\earthengine/Van_Sentinel.tif") as src:
    with rio.open(
            'RGB_Temp.tif', 'w+',
            driver='GTiff',
            dtype= rio.float32,
            count=3,
            crs = src.crs,
            width=src.width,
            height=src.height,
            transform=src.transform,
        ) as dst:
        V = pct_clip(src.read(4))
        dst.write(V,1)
        V = pct_clip(src.read(3))
        dst.write(V,2)
        V = pct_clip(src.read(2))
        dst.write(V,3)
    
fig,ax=plt.subplots()
with rio.open("RGB_Temp.tif") as src2:
    show(src2.read(),transform=src2.transform,ax=ax)

enter image description here

0

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=glob('**/*')

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


#layerstack

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 and acknowledge you have read our privacy policy.

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