用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中打开后作对比，如下
请问是哪里出现了问题？百思不得其解。
