パンシャープン
pan-sharpening
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()
グラムシュミットを用いたパンシャープン
#!/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()
Keyword(s):
References:[リモセンの基礎]