以下是修改后的代码,修复了缩进错误:
from osgeo import gdal # 导入 GDAL 库
import os # 导入操作系统相关库
def calculate_ndvi(input_tif_path, output_tif_path):
"""计算给定 TIF 文件的 NDVI 并输出为新的 TIF 文件
参数:
input_tif_path (str): 输入 TIF 文件的完整路径
output_tif_path (str): 输出 NDVI 文件的完整路径
"""
# 获取 GDAL 驱动
driver = gdal.GetDriverByName('GTiff')
# 打开输入 TIF 文件
ds_ = gdal.Open(input_tif_path)
if ds_ is None: # 如果无法打开文件,打印错误并返回
print(f"Error: Could not open file {input_tif_path}")
return
# 获取输入 TIF 文件的波段数量
num_bands = ds_.RasterCount
print(f"Number of bands: {num_bands}")
# 检查波段数量是否满足要求,如果不满足则打印错误并返回
if num_bands < 4 or num_bands < 5:
print("Error: File does not contain the required number of bands")
return
# 获取红色和近红外波段
band4_red = ds_.GetRasterBand(4)
band5_nir = ds_.GetRasterBand(5)
# 读取红色和近红外波段的数据
red4_array = band4_red.ReadAsArray()
nir5_array = band5_nir.ReadAsArray()
# 计算 NDVI
ndvi = (nir5_array - red4_array) / (nir5_array + red4_array)
# 创建输出 NDVI 文件
out_tif = driver.Create(output_tif_path,
band4_red.XSize,
band5_nir.YSize,
1,
gdal.GDT_Float32)
# 设置输出文件的投影和仿射变换信息
out_tif.SetProjection(ds_.GetProjection())
out_tif.SetGeoTransform(ds_.GetGeoTransform())
# 将 NDVI 数据写入输出文件的第一个波段
outband = out_tif.GetRasterBand(1)
outband.WriteArray(ndvi)
out_tif.FlushCache()
del out_tif # 关闭输出文件
print(f"Output file created: {output_tif_path}")
# 使用示例
input_tif_path=r'D:\yaolaoshi\我国分省土壤数据\吉林.tif'
output_tif_path=r'D:/yaolaoshi/ASTGTM2_N30E120.tif'
calculate_ndvi(input_tif_path, output_tif_path) # 调用 calculate_ndvi 函数
请注意,以上代码中的路径需要根据实际情况进行修改。
内容由零声教学AI助手提供,问题来源于学员提问