关于#Python#计算NDVI的问题

用Python计算Landsat8 OLI影像文件的NDVI值，所得结果为负数

代码如下

# 导入包
from osgeo import gdal
import numpy as np
from osgeo import gdal_array
import os

# 读取
FilePath = 'D:\\Co\\4\\output'
tif_name = os.listdir(FilePath)  # 读取FilePath下的目录名,包含扩展名
SavePath = 'D:\\Co\\4\\output\\outputNDVI'
print(tif_name)

ReadName = 'D:\\Co\\4\\output\\dqjz_daiyue'
SaveName = 'D:\\Co\\4\\output\\outputNDVI\\d9.tif'
dataset = gdal.Open(ReadName)
cols = dataset.RasterXSize  # 列数
rows = dataset.RasterYSize  # 行数  空间坐标 可以理解为行列号

# 计算 NDVI=（NIR-RED）/（NIR+RED）
# 提取红波段
band_red = dataset.GetRasterBand(3)
data_red = band_red.ReadAsArray(0, 0, cols, rows).astype(np.float64)   # 根据位置获取像素值
# 提取近红外波段
band_nir = dataset.GetRasterBand(4)
data_nir = band_nir.ReadAsArray(0, 0, cols, rows).astype(np.float64)
ndvi = ((data_nir - data_red) * 0.1) / (((data_nir + data_red)+ 1E-6) * 0.1)
# 警示错误的运算
gdal_array.numpy.seterr(all="warn")

# 创建栅格数据集
driver = gdal.GetDriverByName('GTiff')  # 将GTiff图纸实体化为driver，括号内为文件类型
dataset_out = driver.Create(SaveName, cols, rows, 1, gdal.GDT_Float64)  # 这里的1指的是创建的栅格数据文件只能容纳一个波段
# 定义仿射矩阵和投影信息
dataset_out.SetGeoTransform(dataset.GetGeoTransform())
dataset_out.SetProjection(dataset.GetProjection())
# 写入像素值
band_out = dataset_out.GetRasterBand(1)
band_out.WriteArray(ndvi)  # 把数组写入到栅格数据相应的波段中
# 清空缓存，关闭数据集
dataset.FlushCache()
dataset_out.FlushCache()  # 把内存中的数据写入到文件中
del dataset
del dataset_out

运行结果及报错内容

在ENVI中打开后作对比，如下

img

请问是哪里出现了问题？百思不得其解。