python+gdal+numpy计算ndvi

from osgeo import gdal
import numpy as np


ds: gdal.Dataset = gdal.Open('gf.tif')
# 红波段
red = ds.GetRasterBand(3).ReadAsArray() * 0.0001
# 近红外波段
nir = ds.GetRasterBand(4).ReadAsArray() * 0.0001

np.seterr(all="ignore")

ndvi = (nir - red) / (nir + red)
# NAN——>0
nan_index = np.isnan(ndvi)
ndvi[nan_index] = 0
ndvi = ndvi.astype(np.float32)
# 创建tif
tif_driver = gdal.GetDriverByName('GTiff')
out_ds = tif_driver.Create('ndvi.tif', ds.RasterXSize, ds.RasterYSize, 1, gdal.GDT_Float32)
# 设置投影坐标
out_ds.SetProjection(ds.GetProjection())
out_ds.SetGeoTransform(ds.GetGeoTransform())
# 写入数据
out_band = out_ds.GetRasterBand(1)
out_band.WriteArray(ndvi)