如何统计多幅GeoTIFF影像中各像素的覆盖次数？

你的核心问题是原代码计算的是全局有效像素总数，而非单个像素的覆盖次数。要实现每个像素的覆盖计数，需要逐幅读取影像，对每幅的有效区域进行累加。以下是两种可行的实现方案：

方案一：使用Rasterio实现 #

1. 定义所有待处理的GeoTIFF文件路径列表
2. 以第一幅影像为基准，创建与影像尺寸一致的计数数组，初始化为0
3. 遍历每幅影像，读取数据并识别有效像素（非NoData值），对计数数组对应位置累加1
4. 将无覆盖区域设为NaN，最后生成热力图

示例代码：

import rasterio
import numpy as np
import matplotlib.pyplot as plt

# 所有待处理的GeoTIFF文件路径
tiff_files = ["image1.tif", "image2.tif", "image3.tif"]

# 读取第一幅影像获取元数据和尺寸
with rasterio.open(tiff_files[0]) as src:
    meta = src.meta.copy()
    count_array = np.zeros((src.height, src.width), dtype=np.int16)
    nodata = src.nodata

# 遍历每幅影像进行计数累加
for tiff in tiff_files:
    with rasterio.open(tiff) as src:
        data = src.read(1)
        # 标记有效像素（非NoData）并累加
        valid_mask = ~np.isnan(data) if np.isnan(nodata) else (data != nodata)
        count_array[valid_mask] += 1

# 将无覆盖的区域设为NaN
count_array[count_array == 0] = np.nan

# 生成热力图
plt.figure(figsize=(10, 8))
im = plt.imshow(count_array, cmap='viridis', interpolation='nearest')
plt.colorbar(im, label='像素覆盖次数')
plt.title('GeoTIFF像素覆盖次数热力图')
plt.axis('off')
plt.show()

# 可选：将计数结果保存为新的GeoTIFF
meta.update(dtype='int16', nodata=np.nan)
with rasterio.open('coverage_count.tif', 'w', **meta) as dst:
    dst.write(count_array, 1)

方案二：使用GDAL实现 #

如果习惯使用GDAL，可采用以下方式实现：

from osgeo import gdal
import numpy as np
import matplotlib.pyplot as plt

tiff_files = ["image1.tif", "image2.tif", "image3.tif"]

# 打开第一幅影像获取基础信息
ds = gdal.Open(tiff_files[0])
band = ds.GetRasterBand(1)
count_array = np.zeros((band.YSize, band.XSize), dtype=np.int16)
nodata = band.GetNoDataValue()

# 遍历所有影像累加计数
for tiff in tiff_files:
    ds = gdal.Open(tiff)
    band = ds.GetRasterBand(1)
    data = band.ReadAsArray()
    valid_mask = (data != nodata)
    count_array[valid_mask] += 1
    ds = None  # 释放资源

# 处理无覆盖区域
count_array[count_array == 0] = np.nan

# 绘制热力图
plt.figure(figsize=(10, 8))
im = plt.imshow(count_array, cmap='viridis', interpolation='nearest')
plt.colorbar(im, label='像素覆盖次数')
plt.title('GeoTIFF像素覆盖次数热力图')
plt.axis('off')
plt.show()

# 可选：保存计数结果为新GeoTIFF
driver = gdal.GetDriverByName('GTiff')
out_ds = driver.Create('coverage_count_gdal.tif', count_array.shape[1], count_array.shape[0], 1, gdal.GDT_Int16)
out_ds.SetProjection(ds.GetProjection())
out_ds.SetGeoTransform(ds.GetGeoTransform())
out_band = out_ds.GetRasterBand(1)
out_band.SetNoDataValue(np.nan)
out_band.WriteArray(count_array)
out_band.FlushCache()
out_ds = None

关键注意事项 #

• 确保所有GeoTIFF影像具有相同的投影、分辨率和地理范围，否则需先通过rasterio.reproject或GDAL的gdalwarp工具进行对齐处理
• 处理NoData值时，需注意不同影像的NoData定义（如-9999、NaN等），确保有效像素识别准确
• 计数数组的dtype可根据影像数量调整，若覆盖次数超过32767，需改用np.int32

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