基于Xarray计算雾频率地图：排除无效值的实现难题

需求说明 #

使用rioxarray.open_rasterio读取多幅geoTIFF格式的雾图，计算每个像素的雾频率：雾频率 = 该像素出现雾的次数 / 该像素的有效观测次数。其中像素值定义为：1=有雾，0=无雾，-9999=无效值（不计入统计）。

原尝试用where过滤无效值后，累加结果出现异常值1.79769e+308（float32最大值），需要修正处理逻辑。

原实现代码 #

# open all fog maps and create a list:
folder = "E:/Jasper/Studium/BA_Thesis/MODIS_data/MODIS_2021_data/2021_06/fog_frequency"
list_of_maps = glob.glob(folder + '/fog_map*.tif', recursive=True)  # all files that start with "fog_map"

# make list with all different filenames (dates) in this folder:
maps = []  # initialize empty list for all file names
for i in range(0, np.size(list_of_maps)):
    # files naming convention "fog_map_YYYYMMDD_HHMMSS.tif":
    maps.append(list_of_maps[i].split('fog_map_')[1][0:8])

# find out how many dates are in the folder:
maps = np.unique(maps)  # remove duplicates from array
print(maps)
print('\ndata from {} different dates in this folder\n'.format(np.size(maps)))

# create fog_sum xarray dataArray to have something to start out with and later subtract it again:
fog_sum = rioxarray.open_rasterio("E:/Jasper/Studium/BA_Thesis/MODIS_data/MODIS_2021_data/2021_06/fog_frequency/fog_map_20210601.tif")
fog_sum_subtract = rioxarray.open_rasterio("E:/Jasper/Studium/BA_Thesis/MODIS_data/MODIS_2021_data/2021_06/fog_frequency/fog_map_20210601.tif")

# add all fog maps:
for i in range(0, np.size(list_of_maps)):
    # open data sets:
    fog_map = rioxarray.open_rasterio(list_of_maps[i], engine='rasterio')
    # fog_map = fog_map.where(fog_map >= 0)
    fog_sum = fog_sum + fog_map

# subtract original fog map and export as geoTIFF:
fog_sum = fog_sum - fog_sum_subtract
fog_sum.rio.to_raster("E:/Jasper/Studium/BA_Thesis/MODIS_data/MODIS_2021_data/2021_06/fog_frequency/fog_sum.tif",
                      driver="GTiff")

问题原因 #

直接用where(fog_map >=0)将无效值转为nan后，与初始包含-9999的fog_sum累加时，会触发数值溢出（-9999与nan的运算逻辑导致异常值）。此外，原代码未单独统计有效像素数量，无法直接计算频率。

修正后的实现代码 #

import glob
import numpy as np
import rioxarray

folder = "E:/Jasper/Studium/BA_Thesis/MODIS_data/MODIS_2021_data/2021_06/fog_frequency"
list_of_maps = glob.glob(folder + '/fog_map*.tif', recursive=True)

# 初始化统计数组：先读取第一幅图作为模板
first_map = rioxarray.open_rasterio(list_of_maps[0])
fog_count = first_map.where(first_map == 1, 0)  # 初始雾计数：雾的位置为1，其他0
valid_count = first_map.where(first_map >= 0, 0).notnull().astype(int)  # 初始有效像素计数

# 循环处理剩余所有图
for map_path in list_of_maps[1:]:
    fog_map = rioxarray.open_rasterio(map_path)
    # 过滤无效值，转换为可统计的格式
    filtered_map = fog_map.where(fog_map >= 0)
    # 累加雾出现的次数：仅当像素为1时加1，否则加0
    fog_count += filtered_map.where(filtered_map == 1, 0)
    # 累加有效像素次数：有效像素（0/1）加1，无效加0
    valid_count += filtered_map.notnull().astype(int)

# 计算雾频率：雾次数/有效次数，避免除以0的情况
fog_frequency = fog_count / valid_count
fog_frequency = fog_frequency.where(valid_count > 0, np.nan)  # 无有效观测的像素设为nan

# 保留原空间参考信息
fog_frequency.rio.set_crs(first_map.rio.crs, inplace=True)
fog_frequency.rio.set_transform(first_map.rio.transform(), inplace=True)

# 导出为geoTIFF
fog_frequency.rio.to_raster(
    folder + "/fog_frequency.tif",
    driver="GTiff",
    dtype="float32"
)

print("雾频率地图已生成")

代码说明 #

• 分别维护fog_count（雾出现次数）和valid_count（有效观测次数）两个数组，避免无效值干扰累加
• 每次循环仅对有效像素进行统计，无效值通过where转为nan后，用notnull()判断是否计入有效次数
• 计算频率时处理除以0的情况，无有效观测的像素设为nan
• 保留原数据的空间参考信息，确保导出的TIFF地理坐标正确

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