Python中使用OpenCV库来进行简单的气象学遥感影像计算

作者:GISLU 时间:2021-02-02 09:45:49 

OpenCV的全称是Open Source Computer Vision Library,是一个跨平台的计算机视觉库。OpenCV是由英特尔公司发起并参与开发,以BSD许可证授权发行,可以在商业和研究领域中免费使用。OpenCV可用于开发实时的图像处理、计算机视觉以及模式识别程序。该程序库也可以使用英特尔公司的IPP进行加速处理。
OpenCV用C++语言编写,它的主要接口也是C++语言,但是依然保留了大量的C语言接口。该库也有大量的Python, Java and MATLAB/OCTAVE (版本2.5)的接口。这些语言的API接口函数可以通过在线文档获得。现在也提供对于C#, Ch,Ruby的支持。
在Windows上编译OpenCV中与摄像输入有关部分时,需要DirectShow SDK中的一些基类。该SDK可以从预先编译的Microsoft Platform SDK (or DirectX SDK 8.0 to 9.0c / DirectX Media SDK prior to 6.0)的子目录Samples\Multimedia\DirectShow\BaseClasses获得。

下面我们就来看看OpenCV在Python编程下的应用,我们来处理一下简单的气象学计算,用python里面的opencv库写个脚本批处理图像反射率的计算试试~

核心步骤就是 遥感影像光谱辐射定标 →大气校正→计算反射率这三步了

1、遥感影像的光谱辐射定标
由遥感器的灵敏度特征引起的辐射畸变主要由其光学系统或光电转换系统的特征形成的,光电转换系统的灵敏性特征通常很重复,其校正一般是通过定期的地面测定值进行的。
遥感器光谱辐射定标时采用以下转换算式:

Python中使用OpenCV库来进行简单的气象学遥感影像计算

遥感器各波段偏移与增益值从论文找了找后,找到这么一张表~

Python中使用OpenCV库来进行简单的气象学遥感影像计算

那么这么个函数就能定标咯:


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)

Python中使用OpenCV库来进行简单的气象学遥感影像计算

其中: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+数据地面反射率反演的数学模型可综合表达为:

Python中使用OpenCV库来进行简单的气象学遥感影像计算

其中:ρ——地面相对反射率;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)
标签:Python,OpenCV
0
投稿

猜你喜欢

  • python开发之thread实现布朗运动的方法

    2021-04-10 00:17:36
  • asp.net下利用js实现返回上一页的实现方法小集

    2024-05-28 15:37:22
  • python3实现ftp服务功能(服务端 For Linux)

    2022-05-18 06:05:03
  • Python海龟绘图详解

    2022-07-26 02:02:02
  • 详细解读tornado协程(coroutine)原理

    2022-08-21 00:20:18
  • 解决SQL Server日志文件损坏严重的问题

    2009-02-05 15:55:00
  • Python+OpenCV 实现简单的高斯滤波(推荐)

    2021-07-18 20:31:26
  • python实现层次聚类的方法

    2023-05-03 22:26:13
  • mysql 正则表达式查询含有非数字和字符的记录

    2024-01-15 11:26:21
  • 巧用Dreamweaver4文件库更新网站

    2007-02-03 11:31:00
  • python使用requests实现发送带文件请求功能

    2023-11-03 14:23:13
  • Vision Transformer图像分类模型导论

    2022-09-09 08:37:39
  • 通过自学python能找到工作吗

    2021-07-24 04:26:33
  • 基于python生成器封装的协程类

    2022-09-25 04:57:53
  • Python Socket编程入门教程

    2022-03-08 01:08:49
  • 介绍Python的Urllib库的一些高级用法

    2023-03-29 08:28:06
  • python单向循环链表实例详解

    2023-05-25 01:29:40
  • 详解Django配置优化方法

    2023-11-01 05:46:58
  • node.js中的fs.writeFileSync方法使用说明

    2024-05-13 09:29:45
  • vue-cli4.5.x快速搭建项目

    2024-04-27 15:52:18
  • asp之家 网络编程 m.aspxhome.com