Pan-sharpening of Landsat by Gram Schmidt

  • 2017/10/15 Kenlo Nasahara
  • ipython3 on Ubuntu Linux
In [1]:
import numpy as np
from osgeo import gdal, gdal_array
import matplotlib.pyplot as plt     # for display images
In [2]:
def inprod(x, y):                   # inner product. "np.dot" is useless for large vector.
   xv=np.reshape(x,x.size)          # x vector
   yv=np.reshape(y,y.size)          # y vector
   return(np.dot(xv,yv))
In [3]:
def showBGRimage(img):              # make and display RGB image from BGR array
   rgb=np.empty([img.shape[0],img.shape[1],3])
   rgb[:,:,0]=img[:,:,2]
   rgb[:,:,1]=img[:,:,1]
   rgb[:,:,2]=img[:,:,0]
   plt.imshow(rgb/np.max(rgb))
   plt.show()
In [4]:
# read Landsat 8 30m bands 
ds=gdal.Open("LC08_L1TP_107035_20170320_20170328_01_T1_B1.TIF", gdal.GA_ReadOnly)
cwidth = ds.RasterXSize               # "coarse lattice width"
cheight = ds.RasterYSize              # "coarse lattice height"
Lorig=np.empty((cheight, cwidth, 6))  # "Landsat original"

ds=gdal.Open("LC08_L1TP_107035_20170320_20170328_01_T1_B1.TIF", gdal.GA_ReadOnly)
Lorig[:,:,1]=ds.GetRasterBand(1).ReadAsArray()
ds=gdal.Open("LC08_L1TP_107035_20170320_20170328_01_T1_B2.TIF", gdal.GA_ReadOnly)
Lorig[:,:,2]=ds.GetRasterBand(1).ReadAsArray()
ds=gdal.Open("LC08_L1TP_107035_20170320_20170328_01_T1_B3.TIF", gdal.GA_ReadOnly)
Lorig[:,:,3]=ds.GetRasterBand(1).ReadAsArray()
ds=gdal.Open("LC08_L1TP_107035_20170320_20170328_01_T1_B4.TIF", gdal.GA_ReadOnly)
Lorig[:,:,4]=ds.GetRasterBand(1).ReadAsArray()
ds=gdal.Open("LC08_L1TP_107035_20170320_20170328_01_T1_B5.TIF", gdal.GA_ReadOnly)
Lorig[:,:,5]=ds.GetRasterBand(1).ReadAsArray()
ds=gdal.Open("LC08_L1TP_107035_20170320_20170328_01_T1_B6.TIF", gdal.GA_ReadOnly)
In [5]:
# simulate pan band. This must be changed with a response function. Future work!
Lorig[:,:,0]=(Lorig[:,:,1]+Lorig[:,:,2]+Lorig[:,:,3]+Lorig[:,:,4]+Lorig[:,:,5])/5
In [6]:
# cutout
cx0=2500     # "coarse lattice x0"
cy0=2500     # "coarse lattice y0"
cx1=2625
cy1=2625
Lc=Lorig[cy0:cy1, cx0:cx1, :]     # Landsat coarse image
L=Lc.copy()                       # Landsat input image
L[L==65535]=0                     # suppress influence of null value
showBGRimage(L[:,:,2:5])          # watch and check the input image! B2=blue, B3=green, B4=red
In [7]:
# remove DC (bias) component
bias=[np.mean(L[:,:,i]) for i in range(0,L.shape[2])]
for i in range(0,L.shape[2]):
    L[:,:,i]=L[:,:,i]-bias[i]
In [8]:
# Gram-Schmidt transformation: L |-> Lg
Lg=np.empty(L.shape)                    # Landsar Gram-Schmidt image
u=np.zeros([L.shape[2],L.shape[2]])     # coefficients of Gram-Schmidt transformation
for k in range (0,L.shape[2]):
  Lg[:,:,k]=L[:,:,k].copy()
  u[k,k]=1
  for i in range(0,k):
     u[i,k]=inprod(Lg[:,:,k],L[:,:,i])/inprod(Lg[:,:,i],Lg[:,:,i])
     Lg[:,:,k]+=-Lg[:,:,i]*u[i,k]
In [9]:
# upscaling (coarse lattice -> fine lattice)
fwidth=(cx1-cx0)*2-1           # "fine lattice width". "-1" is important!!!
fheight=(cy1-cy0)*2-1          # "fine lattice height". "-1" is important!!!
# interpolate 30m image at fine lattice. Be careful of Landsat pixel alignments!
Lgintp=np.empty((fheight, fwidth, Lorig.shape[2]))
## fine pixel at the center of the coarse pixel:
Lgintp[0::2,0::2,:]= Lg
## fine pixel at the corners of four coarse pixels:
Lgintp[1::2,1::2,:]=(Lg[0:-1,0:-1,:]+Lg[1:,0:-1,:]+Lg[0:-1,1:,:]+Lg[1:,1:,:])/4
## fine pixel at the boundary of two adjucent pixels
Lgintp[1::2,0::2,:]=(Lg[0:-1,:,   :]+Lg[1:,   :,:])/2
## fine pixel at the boundary of two adjucent pixels
Lgintp[0::2,1::2,:]=(Lg[:   ,0:-1,:]+Lg[:,   1:,:])/2
In [10]:
# read Landsat 8 panchromatic band
ds=gdal.Open("LC08_L1TP_107035_20170320_20170328_01_T1_B8.TIF", gdal.GA_ReadOnly)
Lorigpan=ds.GetRasterBand(1).ReadAsArray()      # "Landsat original panchromatic"
fx0=cx0*2     # "fine lattice x0"
fy0=cy0*2     # "fine lattice y0"
fx1=fx0+fwidth
fy1=fy0+fheight
Lpan=Lorigpan[fy0:fy1, fx0:fx1].copy()
Lpan[Lpan==65535]=0                    # suppress influence of null value
bias_Lpan=np.mean(Lpan)
Lpan=Lpan-bias_Lpan
In [11]:
# inverse Gram-Schmidt transformation
Lgintp[:,:,0]=Lpan.copy()
Lf=np.empty(Lgintp.shape)   # Landsat fine image (as a result of Gram-Schmidt inverse)
for k in range (1,Lf.shape[2]):
  Lf[:,:,k]=Lgintp[:,:,k]+bias[k]
  for i in range(0,k):
     Lf[:,:,k]+=Lgintp[:,:,i]*u[i,k]
In [12]:
# display image
showBGRimage(Lc[:,:,2:5])    # original image (non-sharpening)
showBGRimage(Lf[:,:,2:5])    # result image (pan-sharpening)