使用Rasterio保存numpy数组至GeoTIFF全为零的问题求助

处理Landsat 8数据计算植被指数（NIRv）时，已完成波段读取与指数计算，绘图显示结果正常，但保存为GeoTIFF后重新打开读取到的全是零值。数据的空间特征（尺寸、CRS、transform等）均未改变，可复现代码如下：

示例云优化文件路径 #

cogs = ['s3://usgs-landsat/collection02/level-2/standard/oli-tirs/2019/026/030/LC08_L2SP_026030_20190505_20200828_02_T1/LC08_L2SP_026030_20190505_20200828_02_T1_SR_B4.TIF',
 's3://usgs-landsat/collection02/level-2/standard/oli-tirs/2019/026/030/LC08_L2SP_026030_20190505_20200828_02_T1/LC08_L2SP_026030_20190505_20200828_02_T1_SR_B5.TIF',
 's3://usgs-landsat/collection02/level-2/standard/oli-tirs/2019/026/030/LC08_L2SP_026030_20190505_20200828_02_T1/LC08_L2SP_026030_20190505_20200828_02_T1_QA_PIXEL.TIF']

导入依赖包 #

import matplotlib.pyplot as plt
import numpy as np
import rasterio as rio
from rasterio.session import AWSSession
import boto3

读取数据 #

aws_session = AWSSession(boto3.Session(region_name='us-west-2'), requester_pays=True)
with rio.Env(aws_session,
            AWS_NO_SIGN_REQUEST='NO',
            GDAL_DISABLE_READDIR_ON_OPEN='TRUE'):
    with rio.open(cogs[0]) as src:
        profile = src.profile
        red_raw = src.read(1)
    with rio.open(cogs[1]) as src_nir:
        nir_raw = src_nir.read(1)
    with rio.open(cogs[2]) as src_qa:
        qa = src_qa.read(1)

数据处理与云掩膜 #

cloud_free_vals = [21824, 21952, 22080, 22208, 23888, 24144, 30048, 30304,
       54596, 54724, 54852, 54980, 56660, 56916, 62820, 63076]
cloud_mask = np.isin(qa, cloud_free_vals)

red = np.ma.masked_equal(red_raw,0.0)
red = np.ma.array(red, mask=cloud_mask) 
red = red * 0.0000275 + -0.2
red[red>1] = 1
red[red<0] = 0

nir = nir_raw * 0.0000275 + -0.2
nir = np.ma.masked_less(nir, 0.0)
nir[nir>1] = 1

nirv = (nir-red)/(nir+red)*nir

原保存代码（导致零值问题） #

save_example_file = rio.open(
    'test.tif',
    'w',
    **profile)
save_example_file.write(nirv.data, 1)
save_example_file.close()

验证代码 #

test2 = rio.open('test.tif', 'r')
plot_example = test2.read(1)
plt.figure()
plt.imshow(plot_example)
plt.colorbar()
plt.show()

问题原因分析 #

原始Landsat SR波段的profile中数据类型为uint16，但计算得到的nirv是浮点型数组（值范围通常在0~1之间）。直接使用原profile保存时，Rasterio会将浮点值强制转换为uint16类型，所有小于1的浮点值都会被截断为0，最终导致保存的文件全为零值。

此外，原profile中的nodata值是针对uint16设置的，与浮点型数据不兼容，掩码区域的处理也未正确映射到输出文件的nodata。

解决方案 #

修改输出文件的profile，匹配浮点型数据的类型与nodata设置，同时正确处理掩码区域：

修改后的保存代码 #

# 复制原始profile并修改关键参数
output_profile = profile.copy()
# 设置浮点型数据类型（float32平衡精度与存储空间）
output_profile['dtype'] = 'float32'
# 设置浮点型对应的nodata值
output_profile['nodata'] = np.nan
# 确保单波段输出
output_profile['count'] = 1

# 将掩码区域的值设置为nodata
nirv_data = np.where(nirv.mask, np.nan, nirv.data).astype(np.float32)

# 保存文件
with rio.open('test_fixed.tif', 'w', **output_profile) as dst:
    dst.write(nirv_data, 1)

验证修改后的结果 #

# 重新读取验证
with rio.open('test_fixed.tif') as src:
    fixed_data = src.read(1)

plt.figure()
plt.imshow(fixed_data, cmap='viridis')
plt.colorbar()
plt.title('Fixed NIRv Result')
plt.show()

# 检查非零值数量
print(f"非零值数量：{np.count_nonzero(~np.isnan(fixed_data))}")

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