# -*- coding: utf-8 -*-

"""
@File    : 16_波段计算-NDVI为例.py
@Author  : fungis@163.com
@notice  :
"""

import os

from osgeo import gdal_array as ga

# Landset8 NDVI=(float(b5)-float(b4))/(float(b5)+float(b4))
# 红光波段
RED_FILE = r'E:\gisData\LC08_L1TP_123037_20200209_20200211_01_T1\LC08_L1TP_123037_20200209_20200211_01_T1_B4.TIF'
# 近红外波段
NIR_FILE = r'E:\gisData\LC08_L1TP_123037_20200209_20200211_01_T1\LC08_L1TP_123037_20200209_20200211_01_T1_B5.TIF'
# 输出路径
result_target = r'./results/images/ndvi2.tif'

RED_data = ga.LoadFile(RED_FILE).astype("float32")
NIR_data = ga.LoadFile(NIR_FILE).astype("float32")

ga.numpy.seterr(all="ignore")
ndvi = (NIR_data - RED_data) / (RED_data + NIR_data)
result_array = ga.numpy.nan_to_num(ndvi)
if not os.path.exists(result_target[:result_target.rindex('/')]):
    # os.mkdir创建路径
    os.makedirs(result_target[:result_target.rindex('/')])

out = ga.SaveArray(result_array, result_target, format="GTiff", prototype=RED_FILE)  # 投影信息prototype