使用gdal库进行两个栅格图像相除和相减以及（B1 - B2）/(B1 + B2)的操作并且分别存储到不同的栅格图像中的代码

以下是使用gdal库进行两个栅格图像相除和相减以及（B1 - B2）/(B1 + B2)的操作并且分别存储到不同的栅格图像中的Python代码示例：

import gdal
import numpy as np

# 读取栅格图像1
raster1 = gdal.Open("raster1.tif")
band1 = raster1.GetRasterBand(1)
array1 = band1.ReadAsArray()

# 读取栅格图像2
raster2 = gdal.Open("raster2.tif")
band2 = raster2.GetRasterBand(1)
array2 = band2.ReadAsArray()

# 计算相减结果
result_subtract = array1 - array2

# 计算相除结果
result_divide = np.divide(array1, array2, out=np.zeros_like(array1), where=array2!=0)

# 计算（B1 - B2）/(B1 + B2)结果
result_ratio = np.divide((array1 - array2), (array1 + array2), out=np.zeros_like(array1), where=array1!=0)

# 存储相减结果
driver = gdal.GetDriverByName("GTiff")
outRaster = driver.Create("result_subtract.tif", raster1.RasterXSize, raster1.RasterYSize, 1, band1.DataType)
outRaster.SetProjection(raster1.GetProjection())
outRaster.SetGeoTransform(raster1.GetGeoTransform())
outBand = outRaster.GetRasterBand(1)
outBand.WriteArray(result_subtract)
outBand.FlushCache()

# 存储相除结果
outRaster = driver.Create("result_divide.tif", raster1.RasterXSize, raster1.RasterYSize, 1, band1.DataType)
outRaster.SetProjection(raster1.GetProjection())
outRaster.SetGeoTransform(raster1.GetGeoTransform())
outBand = outRaster.GetRasterBand(1)
outBand.WriteArray(result_divide)
outBand.FlushCache()

# 存储（B1 - B2）/(B1 + B2)结果
outRaster = driver.Create("result_ratio.tif", raster1.RasterXSize, raster1.RasterYSize, 1, band1.DataType)
outRaster.SetProjection(raster1.GetProjection())
outRaster.SetGeoTransform(raster1.GetGeoTransform())
outBand = outRaster.GetRasterBand(1)
outBand.WriteArray(result_ratio)
outBand.FlushCache()

# 关闭栅格图像
raster1 = None
raster2 = None
outRaster = None