python处理遥感数据（NVDI计算、辐射校正）

前言：作为遥感、GIS从业者，会用python将能提高你的核心竞争力，同时也证明你可以摆脱专业软件的束缚，迈入自主研发的新阶段。ps:其实ENVI、Arcgis等相关软件已经非常完善了，想要突破比较困难，但是这一小步的迈出也决定了你已跨入自主研发的道路。

1.辐射校正

1.1传感器端辐射校正

1.2大气校正

1.3地表辐射校正

2.代码

3.结束

1.辐射校正

辐射校正的目的是尽可能消除传感器自身条件、薄雾等大气条件、太阳位置和角度条件及某些不可避免的噪声引起的传感器测量值与目标的光谱反射率或者光谱辐亮度等物理量之间的差异。还原图像原来面目，进而保障后续图像处理工作的稳定。

所以辐射校正的内容也是由产生影响的因素决定的即：传感器端辐射校正、大气校正和地表辐射校正。

1.1传感器端辐射校正

（这里拿LANDSAT的TM图像校正方式作为举例，一般遥感软件在读数据时自动读取头文件中的参数信息，进行了这一步操作）

①辐亮度与DN值关系：

$L_{\lambda }=Gain*DN+Bias$

$L_{\lambda }$ ——图像的辐亮度(单位：W/m*sr);Gain——增益；DN——图像亮度值；Bias——偏置。

②辐亮度和大气上界反射率值关系：

$\rho =\frac{\pi \cdot L_{\lambda }\cdot d^{2}}{ESUN_{\lambda }\cdot cos\theta }$

$\rho$ ——反射率；d——日地距离；ESUN——太阳光谱辐射量； $\theta$ ——为太阳天顶角。

1.2大气校正

这里主要有：统计学方法、辐射传递方程计算法、波段对比法等。

寻找黑区域辐射，作为修正值。

$Lhazel=LI_{min}-LI_{1}$ %

$LI_{min}$ ——每一波段最小光谱辐射值； $L_{1}$ %——反射率为1%的黑体辐射值。

1.3地表辐射校正

$f\left ( x,y \right )=\frac{g\left ( x,y \right )}{cos\theta }$

$f(x,y)$ ——直射时图像； $g(x,y)$ ——高度角i时的图像； $\theta$ ——太阳天顶角。

综合以上：

地面反射率校正：

$\rho =\frac{\pi \cdot (L_{\lambda }-Lhazel)\cdot d^{2}}{ESUN_{\lambda }\cdot cos^{2}\theta }$

2.代码

利用GDAL库，结合公式和一些代码，自己总结了下，计算NDVI的简单代码：

import gdal
from osgeo import gdal_array
import  math
#路径
path=r'C:\Users\Desktop\test\path\\'Bias=1.18070871
Gains=-7.38070852
# 计算JD
Year = 2019
Month = 8
Day = 23
JD = Day - 32075 + 1461 * (Year + 4800 + (Month - 14) / 12) / 4 + 367 * (Month - 2 - (Month - 14) / 12 * 12) / 12 - 3 * ((Year + 4900 + (Month - 14) / 12) / 100) / 4
# 设置ESUNI值
ESUNI = 196.9
# 计算日地距离Dist
Dist = 1 - 0.01674 * math.cos((0.9856 * (JD - 4) * math.pi / 180))
# 计算太阳天顶角
cos = math.cos(math.radians(90 - 39.5))
inter = (math.pi * Dist * Dist) / (ESUNI * cos * cos)
# 大气校正参数设置
Lmini = -6.2
Lmax = 293.7
Qcal = 1
Qmax = 255
LIMIN = Lmini + (Qcal * (Lmax - Lmini) / Qmax)
LI = (0.01 * ESUNI * cos * cos) / (math.pi * Dist * Dist)
Lhazel = LIMIN - LIb4=path+'T50SPJ_20190823T025551_B02.jp2'
b8=path+'T50SPJ_20190823T025551_B08.jp2'Red=gdal_array.LoadFile(b4)
NIR=gdal_array.LoadFile(b8)Red=inter*(Bias+Gains*Red-Lhazel)
NIR=inter*(Bias+Gains*NIR-Lhazel)gdal_array.numpy.seterr(all="ignore")
ndvi=((NIR-Red)*1.0)/((NIR+Red)*1.0)
ndvi1=gdal_array.numpy.nan_to_num(ndvi)
target=path+'ndvi.tif'
out=gdal_array.SaveArray(ndvi1,target,format = "GTiff",prototype = b4)
out=None

3.结束

利用两个波段（近红、红），简单进行辐射校正后的ndvi计算，提取植被。

实际上，遥感图像处理对于细节的处理很多，针对不同传感器，参数和模型方法也不同，需要依据需求和场景灵活使用，本文只是简单依据流程进行的NDVI计算过程，之中的不足、错误还望指正。