Python中使用OpenCV庫(kù)來(lái)進(jìn)行簡(jiǎn)單的氣象學(xué)遙感影像計(jì)算

下面我們就來(lái)看看OpenCV在Python編程下的應(yīng)用，我們來(lái)處理一下簡(jiǎn)單的氣象學(xué)計(jì)算，用python里面的opencv庫(kù)寫個(gè)腳本批處理圖像反射率的計(jì)算試試~

核心步驟就是 遙感影像光譜輻射定標(biāo) →大氣校正→計(jì)算反射率這三步了

1、遙感影像的光譜輻射定標(biāo)
由遙感器的靈敏度特征引起的輻射畸變主要由其光學(xué)系統(tǒng)或光電轉(zhuǎn)換系統(tǒng)的特征形成的,光電轉(zhuǎn)換系統(tǒng)的靈敏性特征通常很重復(fù),其校正一般是通過(guò)定期的地面測(cè)定值進(jìn)行的。
遙感器光譜輻射定標(biāo)時(shí)采用以下轉(zhuǎn)換算式:

遙感器各波段偏移與增益值從論文找了找后，找到這么一張表~

那么這么個(gè)函數(shù)就能定標(biāo)咯：

def computL(gain,Dn,bias): 
 return (gain*Dn+bias) 

2、遙感影像的大氣校正
任何一種依賴大氣物理模型的大氣校正方法都需要先進(jìn)行遙感器的輻射校準(zhǔn)。
公式是這個(gè)咯（Chavez P S，Jr. Image -Based Atmospheric Correction Revisited and Improved Photogrammetric Engineering and Remote Sensing, 1996,62,1025 -1036）

其中:Lhazel——大氣層光譜輻射值;LI，min——遙感器每一波段最小光譜輻射值;LI,1%——反射率為1%的黑體輻射值。

關(guān)于LI，min和LI,1%的計(jì)算公式就省略了啊，感興趣的同學(xué)可以自己去查查論文~

而計(jì)算Lhazel需要的參數(shù)可以從遙感圖像的頭文件中獲得一部分，還有一部分是固定的參數(shù)~這些都藏在ENVI的背后，不過(guò)自己寫腳本的時(shí)候找出他們還是廢了一番功夫的。

計(jì)算Lhazel的代碼如下：

#ESUN 
ESUNI71=196.9 
cos=math.cos(math.radians(90-41.3509605)) 
# 
Lmini=-6.2 
Lmax=293.7 
# 
Qcal=1 
Qmax=255 
LIMIN=Lmini+(Qcal*(Lmax-Lmini)/Qmax) 
LI=(0.01*ESUNI71*cos*cos)/(math.pi*D*D) 
Lhazel=LIMIN-LI 

3、計(jì)算遙感影像的反射率
根據(jù)太陽(yáng)輻射和大氣傳輸原理與過(guò)程,TM/ETM+數(shù)據(jù)地面反射率反演的數(shù)學(xué)模型可綜合表達(dá)為:

其中:ρ——地面相對(duì)反射率;D——日地天文單位距離;LsatI——傳感器光譜輻射值,即大氣頂層的輻射能量;LhazeI——大氣層輻射值;ESUNl——大氣頂層的太陽(yáng)平均光譜輻射,即大氣頂層太陽(yáng)輻照度;SZ——太陽(yáng)天頂角。

這里提一下其中兩個(gè)參數(shù)的計(jì)算公式：
日地天文單位距離 D=1 -0.01674 cos(0.9856×(JD-4)×π/180);
（JD為遙感成像的儒略日(Julian Day)，計(jì)算公式為：

JD=K-32075+1461*(I+4800+(J-14)/12)/4+367*(J-2-(J-14)/12*12)/12-3*((I+4900+(J-14)/12)/100)/4 

I、J、K分別為年、月、日

有了這些，最后就能直接算出來(lái)反射率啦，粗糙代碼如下，因?yàn)槭菍懼娴?，也沒(méi)怎么處理：
不過(guò)需要注意的是，遙感圖像進(jìn)行計(jì)算跟輸出的時(shí)候，需要使用uint16類型的數(shù)組來(lái)存儲(chǔ)的（uint8長(zhǎng)度不夠啊。）
一些參數(shù)涉及到浮點(diǎn)數(shù)計(jì)算，如果對(duì)處理結(jié)果有極高要求的話，最好使用專門的科學(xué)運(yùn)算庫(kù)（像我這種渣學(xué)校才不介意這些）

import cv2 
import numpy as np 
import math 
 
img1=cv2.imread('F:L71121040_04020030220_B10.TIF') 
#圖像格式轉(zhuǎn)換 
img10=cv2.cvtColor(img1,cv2.COLOR_BGR2GRAY) 
 
#計(jì)算JD 
I=2003 
J=2 
K=20 
JD=K-32075+1461*(I+4800+ (J-14)/12)/4+367*(J-2-(J-14)/12*12)/12-3*((I+4900+(J-14)/12)/100)/4 
#設(shè)置ESUNI值 
ESUNI71=196.9 
#計(jì)算日地距離D 
D=1-0.01674*math.cos((0.9856*(JD-4)*math.pi/180)) 
#計(jì)算太陽(yáng)天頂角 
cos=math.cos(math.radians(90-41.3509605)) 
 
inter=(math.pi*D*D)/(ESUNI71*cos*cos) 
 
#大氣校正參數(shù)設(shè)置 
Lmini=-6.2 
Lmax=293.7 
Qcal=1 
Qmax=255 
LIMIN=Lmini+(Qcal*(Lmax-Lmini)/Qmax) 
LI=(0.01*ESUNI71*cos*cos)/(math.pi*D*D) 
Lhazel=LIMIN-LI 
 
 
def copy(img,new1): 
 new1= np.zeros(img.shape,dtype='uint16') 
 new1[:,:] = img[:,:] 
 
def computL(gain,Dn,bias): 
 return (gain*Dn+bias) 
 
if __name__ == '__main__': 
 print 'D=',D 
 print 'cosZS=',cos 
 print 'Lhazel=',Lhazel 
 #計(jì)算圖像反射率 
 result=np.zeros(img.shape,dtype='uint16') 
 for i in range(0,img.shape(1)): 
 for j in range(0,img.shape(0)): 
 Lsat=computL(1.18070871,img10[i,j],-7.38070852) 
 result[i,j]=inter*(Lsat-Lhazel)*1000 
 
#保存圖像 
 cv2.imwrite("F:\result.tif", result) 
 cv2.namedWindow("Image") 
 cv2.imshow("Image", result) 
 cv2.waitKey(0) 

聲明：本網(wǎng)頁(yè)內(nèi)容旨在傳播知識(shí)，若有侵權(quán)等問(wèn)題請(qǐng)及時(shí)與本網(wǎng)聯(lián)系，我們將在第一時(shí)間刪除處理。TEL:177 7030 7066 E-MAIL:11247931@qq.com