I'm using Landsat 8 imagery in Python and want to make true color images that would look good in a publication. I have many images, so ideally this would be an automated process where each image comes out looking not too dark, not overexposed, and with good contrast.

What I've been using is a linear stretch of 1.25%, but some of the images still look washed out.

from skimage import exposure
import matplotlib.pyplot as plt
import numpy as np

r = np.random.randint(0, 255, size=10000).reshape((100, 100))
g = np.random.randint(0, 255, size=10000).reshape((100, 100))
b = np.random.randint(0, 255, size=10000).reshape((100, 100))
rgb = np.dstack((r, g, b))

def linearStretch(input, percent):
    pLow, pHigh = np.percentile(input[~np.isnan(input)], (percent, 100 - percent))
    img_rescale = exposure.rescale_intensity(input, in_range=(pLow, pHigh))
    return img_rescale

img_rescaled = linearStretch(rgb, 1.25)

Is there a known image enhancement technique that can be used to reliably process images automatically without much/any individual tweaking?

All of my images have had clouds removed and replaced with NaN.


Here is a histogram equalization solution that is working better than linear stretch, but some images are still a little washed out, not as saturated as I'd like.

# Get RGB bands from Landsat image
with rasterio.open(landsat_path, 'r') as f:
    red, green, blue = f.read(4), f.read(3), f.read(2)
    red[red == -999999] = 0
    green[green == -999999] = 0
    blue[blue == -999999] = 0
    rgb = np.dstack((red, green, blue))

# Some of my images have NaNs for NoData so I have to remove them when creating the histogram
shape = rgb.shape
rgb_vector = rgb.reshape([rgb.shape[0] * rgb.shape[1], rgb.shape[2]])
rgb_vector = rgb_vector[~np.isnan(rgb_vector).any(axis=1)]

# View histogram of RGB values
fig = plt.figure(figsize=(10, 7))
for color, channel in zip('rgb_vector', np.rollaxis(rgb, axis=-1)):
    counts, centers = exposure.histogram(channel)
    plt.plot(centers[1::], counts[1::], color=color)

# Get cutoff values based on standard deviations. Ideally these would be on either side of each histogram peak and cutoff the tail. 
lims = []
for i in range(3):
    x = np.mean(rgb_vector[:, i])
    sd = np.std(rgb_vector[:, i])
    low = x-(1.75*sd)  # Adjust the coefficient here if the image doesn't look right
    high = x + (1.75 * sd)  # Adjust the coefficient here if the image doesn't look right
    if low < 0:
        low = 0
    if high > 1:
        high = 1
    lims.append((low, high))

r = exposure.rescale_intensity(rgb[:, :, 0], in_range=lims[0])
g = exposure.rescale_intensity(rgb[:, :, 1], in_range=lims[1])
b = exposure.rescale_intensity(rgb[:, :, 2], in_range=lims[2])
rgb_enhanced = np.dstack((r, g, b))

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.