如何在Python中将numpy数组转为GeoTIFF？求解决Qcal转TIFF问题

我明白你在把处理好的Qcal（NumPy数组）转成TIFF格式时遇到了碰壁的情况，结合你已经导入gdal的代码上下文，我给你一套可行的落地方案，包含完整代码和关键注意事项：

核心思路 #

利用GDAL库创建TIFF文件驱动，将NumPy数组直接写入其中，同时可以保留原影像的地理参考信息（如果需要的话），确保生成的TIFF能在GIS工具中正常使用。

完整代码实现 #

1. 导入依赖库 #

import numpy as np
from osgeo import gdal, osr

2. 定义转换函数 #

这个函数可以处理单波段和多波段的NumPy数组，还支持传入地理变换和投影信息：

def numpy_to_tiff(array, output_path, geotransform=None, projection=None):
    # 获取数组的行列数
    rows, cols = array.shape[:2]
    # 判断波段数，单波段数组自动转为(1, rows, cols)格式
    bands = array.shape[0] if len(array.shape) == 3 else 1
    if bands == 1:
        array = array.reshape((1, rows, cols))
    
    # 匹配NumPy和GDAL的数据类型，避免数据失真
    dtype_mapping = {
        np.uint8: gdal.GDT_Byte,
        np.uint16: gdal.GDT_UInt16,
        np.int16: gdal.GDT_Int16,
        np.float32: gdal.GDT_Float32,
        np.float64: gdal.GDT_Float64
    }
    data_type = dtype_mapping.get(array.dtype, gdal.GDT_Float32)
    
    # 创建TIFF文件
    driver = gdal.GetDriverByName('GTiff')
    dataset = driver.Create(output_path, cols, rows, bands, data_type)
    
    # 设置地理变换和投影（如果有原影像的参考信息）
    if geotransform is not None:
        dataset.SetGeoTransform(geotransform)
    if projection is not None:
        srs = osr.SpatialReference()
        srs.ImportFromWkt(projection)
        dataset.SetProjection(srs.ExportToWkt())
    
    # 将数组写入每个波段
    for band_idx in range(bands):
        dataset.GetRasterBand(band_idx + 1).WriteArray(array[band_idx])
    
    # 确保数据写入磁盘并释放资源
    dataset.FlushCache()
    dataset = None

3. 结合你的Landsat处理流程使用 #

先补全你读取MTL文件和获取Qcal的代码，再调用转换函数：

# 读取MTL文件获取辐射定标参数
mtl_path = 'LC08_L1TP_180028_20170623_20170630_01_T1_MTL.txt'
mtl_params = {}
with open(mtl_path, 'r') as f:
    for line in f:
        line = line.strip()
        if '=' in line:
            key, val = line.split('=', 1)
            mtl_params[key.strip()] = val.strip().strip('"')

# 提取所需参数
RADIANCE_MULT_BAND_10 = float(mtl_params['RADIANCE_MULT_BAND_10'])
RADIANCE_ADD_BAND_10 = float(mtl_params['RADIANCE_ADD_BAND_10'])
K1_CONSTANT_BAND_10 = float(mtl_params['K1_CONSTANT_BAND_10'])
K2_CONSTANT_BAND_10 = float(mtl_params['K2_CONSTANT_BAND_10'])

# 读取原B10影像得到Qcal数组（或者用你自己计算好的Qcal）
b10_dataset = gdal.Open('LC08_L1TP_180028_20170623_20170630_01_T1_B10.TIF')
Qcal = b10_dataset.ReadAsArray()
# 获取原影像的地理信息（可选，但建议保留）
geo_transform = b10_dataset.GetGeoTransform()
projection = b10_dataset.GetProjection()
b10_dataset = None

# 调用函数生成TIFF
output_tiff_path = 'Qcal_result.tif'
numpy_to_tiff(Qcal, output_tiff_path, geo_transform, projection)

关键注意事项 #

• 数据类型匹配：一定要确保NumPy数组的类型和GDAL的类型对应，上面的函数里已经做了映射处理，避免出现数据截断或异常。
• 地理信息保留：如果你需要生成的TIFF能在GIS软件中正确定位，一定要传入原影像的geo_transform和projection，这些信息可以直接从原Landsat波段文件中读取。
• 资源释放：处理GDAL数据集后，记得将其设为None，避免内存泄漏。

内容的提问来源于stack exchange，提问作者Boris Altynnik