Python中使用OpenCV库来进行简单的气象学遥感影像计算
OpenCV的全称是Open Source Computer Vision Library,是一个跨平台的计算机视觉库。OpenCV是由英特尔公司发起并参与开发,以BSD许可证授权发行,可以在商业和研究领域中免费使用。OpenCV可用于开发实时的图像处理、计算机视觉以及模式识别程序。该程序库也可以使用英特尔公司的IPP进行加速处理。 下面我们就来看看OpenCV在Python编程下的应用,我们来处理一下简单的气象学计算,用python里面的opencv库写个脚本批处理图像反射率的计算试试~ 核心步骤就是 遥感影像光谱辐射定标 →大气校正→计算反射率这三步了 1、遥感影像的光谱辐射定标 遥感器各波段偏移与增益值从论文找了找后,找到这么一张表~ 那么这么个函数就能定标咯: def computL(gain,Dn,bias): return (gain*Dn+bias) 2、遥感影像的大气校正 其中: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
其中:ρ――地面相对反射率;D――日地天文单位距离;LsatI――传感器光谱辐射值,即大气顶层的辐射能量;LhazeI――大气层辐射值;ESUNl――大气顶层的太阳平均光谱辐射,即大气顶层太阳辐照度;SZ――太阳天顶角。 这里提一下其中两个参数的计算公式: 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分别为年、月、日 有了这些,最后就能直接算出来反射率啦,粗糙代码如下,因为是写着玩的,也没怎么处理: 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,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) (编辑:李大同) 【声明】本站内容均来自网络,其相关言论仅代表作者个人观点,不代表本站立场。若无意侵犯到您的权利,请及时与联系站长删除相关内容! |