とらりもんHOME  Index  Search  Changes  Login

とらりもん - パンシャープン Diff

  • Added parts are displayed like this.
  • Deleted parts are displayed like this.

pan-sharpening

{{toc}}

! パンシャープンとは

パンシャープンを体感するのはこちら

*[[Landsat画像を用いたパンシャープン画像の作成]]

!GDALコマンドでパンシャープン
$ gdal_pansharpen.py LC08_L1TP_107035_20170320_20170328_01_T1_B8.TIF LC08_L1TP_107035_20170320_20170328_01_T1_B2.TIF LC08_L1TP_107035_20170320_20170328_01_T1_B3.TIF LC08_L1TP_107035_20170320_20170328_01_T1_B4.TIF outi.tif


! いろんなパンシャープンを理解する
ここでは古典的なパンシャープン手法を、基本的な考え方(数式)とプログラムを交えて理解することを目標とする。

!!PCAを用いたパンシャープン

!!HSVを用いたパンシャープン
#!/usr/bin/env python3
# pan-sharpen by HSV transformation for Landsat data
# 2017/10/13 Kenlo Nishida Nasahara
# 2019/03/18 Kenlo Nishida Nasahara
import numpy as np
from osgeo import gdal, gdal_array
from PIL import Image
import matplotlib.pyplot as plt
import cv2

#--------
# parameter settings
fn_head="LC08_L1TP_107035_20170320_20170328_01_T1_"
cx0, cx1 = 2500, 2600     # "coarse lattice x range"
cy0, cy1 = 2500, 2600     # "coarse lattice y range"

#--------
def showRGBimage(img):              # make and display RGB image from BGR array
    plt.imshow(img/np.max(img))

#---------
print("Reading Landsat RGB images ...")
ds=gdal.Open(fn_head+"B2.TIF", gdal.GA_ReadOnly)
Lb=ds.GetRasterBand(1).ReadAsArray()
ds=gdal.Open(fn_head+"B3.TIF", gdal.GA_ReadOnly)
Lg=ds.GetRasterBand(1).ReadAsArray()
ds=gdal.Open(fn_head+"B4.TIF", gdal.GA_ReadOnly)
Lr=ds.GetRasterBand(1).ReadAsArray()

# coarse coordinate
gt = ds.GetGeoTransform()
cwidth = ds.RasterXSize
cheight = ds.RasterYSize
cminx = gt[0]
cminy = gt[3] + cwidth*gt[4] + cheight*gt[5]
cmaxx = gt[0] + cwidth*gt[1] + cheight*gt[2]
cmaxy = gt[3]
cdx=(cmaxx-cminx)/cwidth
cdy=(cmaxy-cminy)/cheight

Lrgb = np.zeros((Lb.shape[0], Lb.shape[1], 3), dtype=np.uint16)
Lrgb[:,:,0]=Lr
Lrgb[:,:,1]=Lg
Lrgb[:,:,2]=Lb

#---------
print("Reading Landsat panchromatic image")
ds=gdal.Open(fn_head+"B8.TIF", gdal.GA_ReadOnly)
Lp=ds.GetRasterBand(1).ReadAsArray()
# fine coordinate
gt = ds.GetGeoTransform()
fwidth = ds.RasterXSize
fheight = ds.RasterYSize
fminx = gt[0]
fminy = gt[3] + fwidth*gt[4] + fheight*gt[5]
fmaxx = gt[0] + fwidth*gt[1] + fheight*gt[2]
fmaxy = gt[3]
fdx=(fmaxx-fminx)/fwidth
fdy=(fmaxy-fminy)/fheight

#---------
print("Landsat RGB image interpolation")
# this may be wring ... interpolation may not be necessary...
Lrgbi=np.empty((fheight, fwidth, 3))
Lrgbi[0::2,0::2,:]= Lrgb
Lrgbi[1::2,1::2,:]=(Lrgb[0:-1,0:-1,:]+Lrgb[1:,0:-1,:]+Lrgb[0:-1,1:,:]+Lrgb[1:,1:,:])/4
Lrgbi[1::2,0::2,:]=(Lrgb[0:-1,:,   :]+Lrgb[1:,   :,:])/2
Lrgbi[0::2,1::2,:]=(Lrgb[:   ,0:-1,:]+Lrgb[:,   1:,:])/2

sourceRGB=Lrgbi[cy0*2:cy1*2,cx0*2:cx1*2,:].astype(np.float32)
sourceV=Lp[cy0*2:cy1*2,cx0*2:cx1*2].astype(np.float32)

#---------
print("RGB -> HSV transformation")
hsv=cv2.cvtColor(sourceRGB, cv2.COLOR_RGB2HSV)
coef=np.max(hsv[:,:,2])/np.max(sourceV)

#---------
print("HSV -> RGB transformation")
hsv[:,:,2]=sourceV*coef
destinRGB=cv2.cvtColor(hsv, cv2.COLOR_HSV2RGB)

#---------
print("Display image...")
plt.subplot(121)
plt.title("original")
showRGBimage(sourceRGB)
plt.subplot(122)
plt.title("pan-sharpen")
showRGBimage(destinRGB)
plt.show()


[[data4download/L8_pansharpen_HSV.jpg]]

!!グラムシュミットを用いたパンシャープン

#!/usr/bin/env python3
# pan-sharpen by Gram-Schmidt orthogonalization for Landsat data
# 2017/10/17 Kenlo Nishida Nasahara
# 2019/03/18 Kenlo Nishida Nasahara
# just run this script with no parameters!
# then you will get a figure with two images: before and after the pan-sharpening

import numpy as np
from osgeo import gdal, gdal_array
import matplotlib.pyplot as plt     # for display images

#--------
# parameter settings
NULV=65535
fn_head="LC08_L1TP_107035_20170320_20170328_01_T1_"
cx0, cx1 = 2500, 2600     # "coarse lattice x range"
cy0, cy1 = 2500, 2600     # "coarse lattice y range"

#--------
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], rgb[:,:,1], rgb[:,:,2] = img[:,:,2], img[:,:,1], img[:,:,0]
    plt.imshow(rgb/np.max(rgb))
#   plt.show()

#---------
print("Reading images ...")

ds=gdal.Open(fn_head+"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(fn_head+"B1.TIF", gdal.GA_ReadOnly)
Lorig[:,:,1]=ds.GetRasterBand(1).ReadAsArray()
ds=gdal.Open(fn_head+"B2.TIF", gdal.GA_ReadOnly)
Lorig[:,:,2]=ds.GetRasterBand(1).ReadAsArray()
ds=gdal.Open(fn_head+"B3.TIF", gdal.GA_ReadOnly)
Lorig[:,:,3]=ds.GetRasterBand(1).ReadAsArray()
ds=gdal.Open(fn_head+"B4.TIF", gdal.GA_ReadOnly)
Lorig[:,:,4]=ds.GetRasterBand(1).ReadAsArray()
ds=gdal.Open(fn_head+"B5.TIF", gdal.GA_ReadOnly)
Lorig[:,:,5]=ds.GetRasterBand(1).ReadAsArray()
ds=gdal.Open(fn_head+"B6.TIF", gdal.GA_ReadOnly)

#---------
print("Simulating 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

#---------
print("Cutting out images...")
Lc=Lorig[cy0:cy1, cx0:cx1, :]     # Landsat coarse image
L=Lc.copy()                       # Landsat input image
L[L==NULV]=0                     # suppress influence of null value
#showBGRimage(L[:,:,2:5])          # watch and check the input image! B2=blue, B3=green, B4=red

#---------
print("Removing 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]

#---------
print("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

#---------
print("Reading panchromatic band...")
ds=gdal.Open(fn_head+"B8.TIF", gdal.GA_ReadOnly)
Lorigpan=ds.GetRasterBand(1).ReadAsArray()      # "Landsat original panchromatic"
fx0, fy0 = cx0*2, cy0*2     # "fine lattice x0 and y0"
fx1, fy1 = fx0+fwidth, fy0+fheight
Lpan=Lorigpan[fy0:fy1, fx0:fx1].copy()
Lpan[Lpan==NULV]=0                    # suppress influence of null value
bias_Lpan=np.mean(Lpan)
Lpan=Lpan-bias_Lpan

#---------
print("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
plt.subplot(121)
plt.title("original")
showBGRimage(Lc[:,:,2:5])    # original image (non-sharpening)
plt.subplot(122)
plt.title("pan-sharpen")
showBGRimage(Lf[:,:,2:5])    # result image (pan-sharpening)
plt.show()

[[data4download/L8_pansharpen_GramSchmidt.jpg]]