Pan-sharpening of Landsat by RGB <-> HSV

• 2017/10/14 Kenlo Nasahara
• ipython3 on Ubuntu Linux

import numpy as np
from osgeo import gdal, gdal_array
from PIL import Image
import matplotlib.pyplot as plt     # for display images
import cv2                          # openCV library, for RGB <-> HSV transform

# read Landsat 8 RGB bands 
ds=gdal.Open("LC08_L1TP_107035_20170320_20170328_01_T1_B2.TIF", gdal.GA_ReadOnly)
Lb=ds.GetRasterBand(1).ReadAsArray()  # Landsat blue band
ds=gdal.Open("LC08_L1TP_107035_20170320_20170328_01_T1_B3.TIF", gdal.GA_ReadOnly)
Lg=ds.GetRasterBand(1).ReadAsArray()  # Landsat green band
ds=gdal.Open("LC08_L1TP_107035_20170320_20170328_01_T1_B4.TIF", gdal.GA_ReadOnly)
Lr=ds.GetRasterBand(1).ReadAsArray()  # Landsat red band

# combine RGB bands into one array
Lrgb = np.zeros((Lb.shape[0], Lb.shape[1], 3), dtype=np.uint16)
Lrgb[:,:,0]=Lr
Lrgb[:,:,1]=Lg
Lrgb[:,:,2]=Lb

# read Landsat 8 panchromatic band
ds=gdal.Open("LC08_L1TP_107035_20170320_20170328_01_T1_B8.TIF", gdal.GA_ReadOnly)
Lp=ds.GetRasterBand(1).ReadAsArray()

# Get image size for fine lattice (panchromatic band)
fwidth = ds.RasterXSize
fheight = ds.RasterYSize

# interpolate RGB image at fine lattice
Lrgbi=np.empty((fheight, fwidth, 3))
## fine pixel at the center of the coarse pixel:
Lrgbi[0::2,0::2,:]= Lrgb
## fine pixel at the corners of four coarse pixels:
Lrgbi[1::2,1::2,:]=(Lrgb[0:-1,0:-1,:]+Lrgb[1:,0:-1,:]+Lrgb[0:-1,1:,:]+Lrgb[1:,1:,:])/4
## fine pixel at the boundary of two adjucent pixels
Lrgbi[1::2,0::2,:]=(Lrgb[0:-1,:,   :]+Lrgb[1:,   :,:])/2
## fine pixel at the boundary of two adjucent pixels
Lrgbi[0::2,1::2,:]=(Lrgb[:   ,0:-1,:]+Lrgb[:,   1:,:])/2

# RGB -> HSV transform
sourceRGB=Lrgbi.astype(np.float32)
hsv=cv2.cvtColor(sourceRGB, cv2.COLOR_RGB2HSV)

# replace V with panchromatic band
panV=Lp.astype(np.float32)
coef=np.max(hsv[:,:,2])/np.max(panV)
hsv[:,:,2]=panV*coef

# HSV -> RGB transform
destinRGB=cv2.cvtColor(hsv, cv2.COLOR_HSV2RGB)

display image

# original image (non-sharpening)
img=sourceRGB[5000:5250,5000:5250,:]
plt.imshow(img/np.max(img))
plt.show()

# result image (pan-sharpening)
img=destinRGB[5000:5250,5000:5250,:]
plt.imshow(img/np.max(img))
plt.show()