GeoTIFF文件中地理空间数据混合数据类型的处理方案

问题场景 #

处理2字节双波段GeoTIFF影像时，需要生成Float32格式的比值派生波段（第一波段/第二波段），但直接堆叠保存会导致所有波段被统一转换为Float32类型。要求在单个GeoTIFF文件中保留原双波段的2字节数据类型，同时新增Float32格式的第三波段，使用Python的GDAL和Rasterio库实现。

尝试代码（Rasterio） #

with rasterio.open("image.tif") as src:
    data = src.read()
    profile = src.profile
    meta = src.meta.copy()
    band_count = src.count

third_band = (data[0].astype('float32'))/((data[1].astype('float32')) + 1e-6)

with rasterio.open("output_img.tif", "w", **meta) as dst:
  dst.write(data[0], 1)
  dst.write(data[1], 2)
  dst.write(third_band, 3)

解决方法 #

方法一：Rasterio逐波段指定数据类型 #

Rasterio默认元数据中的dtype是全局统一设置，需删除全局dtype配置，在写入每个波段时单独指定对应的数据类型：

import rasterio

with rasterio.open("image.tif") as src:
    # 单独读取原波段数据
    band1 = src.read(1)
    band2 = src.read(2)
    # 复制元数据并更新波段数
    meta = src.meta.copy()
    meta.update(count=3)
    # 删除全局dtype，允许逐波段设置类型
    del meta['dtype']

# 计算Float32格式的比值波段
third_band = (band1.astype('float32')) / (band2.astype('float32') + 1e-6)

with rasterio.open("output_img.tif", "w", **meta) as dst:
    # 写入原波段，指定原数据类型
    dst.write(band1, 1, dtype=src.dtypes[0])
    dst.write(band2, 2, dtype=src.dtypes[1])
    # 写入第三波段，指定Float32类型
    dst.write(third_band, 3, dtype='float32')

方法二：GDAL逐波段定义数据类型 #

GDAL支持在创建数据集时为每个波段单独指定数据类型，步骤如下：

from osgeo import gdal

# 打开原影像
src_ds = gdal.Open("image.tif")
if src_ds is None:
    raise Exception("无法打开目标影像文件")

# 获取原影像的地理参数
proj = src_ds.GetProjection()
geotrans = src_ds.GetGeoTransform()
width = src_ds.RasterXSize
height = src_ds.RasterYSize
# 获取原波段的原生数据类型
band1_dtype = src_ds.GetRasterBand(1).DataType
band2_dtype = src_ds.GetRasterBand(2).DataType

# 创建输出GeoTIFF数据集
driver = gdal.GetDriverByName("GTiff")
dst_ds = driver.Create("output_img.tif", width, height, 3, options=["COMPRESS=LZW"])
if dst_ds is None:
    raise Exception("无法创建输出影像文件")

# 设置投影和地理变换
dst_ds.SetProjection(proj)
dst_ds.SetGeoTransform(geotrans)

# 写入第一波段（保留原数据类型）
dst_band1 = dst_ds.GetRasterBand(1)
dst_band1.SetRasterDataType(band1_dtype)
dst_band1.WriteArray(src_ds.GetRasterBand(1).ReadAsArray())

# 写入第二波段（保留原数据类型）
dst_band2 = dst_ds.GetRasterBand(2)
dst_band2.SetRasterDataType(band2_dtype)
dst_band2.WriteArray(src_ds.GetRasterBand(2).ReadAsArray())

# 计算并写入第三波段（Float32类型）
band1_data = src_ds.GetRasterBand(1).ReadAsArray().astype('float32')
band2_data = src_ds.GetRasterBand(2).ReadAsArray().astype('float32')
third_band = band1_data / (band2_data + 1e-6)
dst_band3 = dst_ds.GetRasterBand(3)
dst_band3.SetRasterDataType(gdal.GDT_Float32)
dst_band3.WriteArray(third_band)

# 释放资源
src_ds = None
dst_ds = None

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