下面我們就來看看OpenCV在Python編程下的應用,我們來處理一下簡單的氣象學計算,用python里面的opencv庫寫個腳本批處理圖像反射率的計算試試~
核心步驟就是 遙感影像光譜輻射定標 →大氣校正→計算反射率這三步了
1、遙感影像的光譜輻射定標
由遙感器的靈敏度特征引起的輻射畸變主要由其光學系統或光電轉換系統的特征形成的,光電轉換系統的靈敏性特征通常很重復,其校正一般是通過定期的地面測定值進行的。
遙感器光譜輻射定標時采用以下轉換算式:
遙感器各波段偏移與增益值從論文找了找后,找到這么一張表~
那么這么個函數就能定標咯:
def computL(gain,Dn,bias): return (gain*Dn+bias)
2、遙感影像的大氣校正
任何一種依賴大氣物理模型的大氣校正方法都需要先進行遙感器的輻射校準。
公式是這個咯(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%的黑體輻射值。
關于LI,min和LI,1%的計算公式就省略了啊,感興趣的同學可以自己去查查論文~
而計算Lhazel需要的參數可以從遙感圖像的頭文件中獲得一部分,還有一部分是固定的參數~這些都藏在ENVI的背后,不過自己寫腳本的時候找出他們還是廢了一番功夫的。
計算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、計算遙感影像的反射率
根據太陽輻射和大氣傳輸原理與過程,TM/ETM+數據地面反射率反演的數學模型可綜合表達為:
其中:ρ——地面相對反射率;D——日地天文單位距離;LsatI——傳感器光譜輻射值,即大氣頂層的輻射能量;LhazeI——大氣層輻射值;ESUNl——大氣頂層的太陽平均光譜輻射,即大氣頂層太陽輻照度;SZ——太陽天頂角。
這里提一下其中兩個參數的計算公式:
日地天文單位距離 D=1 -0.01674 cos(0.9856×(JD-4)×π/180);
(JD為遙感成像的儒略日(Julian Day),計算公式為:
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分別為年、月、日
有了這些,最后就能直接算出來反射率啦,粗糙代碼如下,因為是寫著玩的,也沒怎么處理:
不過需要注意的是,遙感圖像進行計算跟輸出的時候,需要使用uint16類型的數組來存儲的(uint8長度不夠啊。)
一些參數涉及到浮點數計算,如果對處理結果有極高要求的話,最好使用專門的科學運算庫(像我這種渣學校才不介意這些)
import cv2 import numpy as np import math img1=cv2.imread('F:L71121040_04020030220_B10.TIF') #圖像格式轉換 img10=cv2.cvtColor(img1,cv2.COLOR_BGR2GRAY) #計算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 #設置ESUNI值 ESUNI71=196.9 #計算日地距離D D=1-0.01674*math.cos((0.9856*(JD-4)*math.pi/180)) #計算太陽天頂角 cos=math.cos(math.radians(90-41.3509605)) inter=(math.pi*D*D)/(ESUNI71*cos*cos) #大氣校正參數設置 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 #計算圖像反射率 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)
聲明:本網頁內容旨在傳播知識,若有侵權等問題請及時與本網聯系,我們將在第一時間刪除處理。TEL:177 7030 7066 E-MAIL:11247931@qq.com