import numpy as np
from osgeo import gdal, gdal_array
import matplotlib.pyplot as plt # for display images
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))
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()
# 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)
# 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
# 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
# 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]
# 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]
# 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
# 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
# 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]
# display image
showBGRimage(Lc[:,:,2:5]) # original image (non-sharpening)
showBGRimage(Lf[:,:,2:5]) # result image (pan-sharpening)