使用xr.concat堆叠TIF生成NC文件时图像变形的解决方法

使用以下Python代码将TIF文件合并为NC文件时，处理某数据集正常，但处理第二个数据集时出现图像变形：原TIF尺寸为1529×2777，生成的NC文件尺寸变为3058×5554，且两个数据集投影一致。

原代码：

import glob
import xarray as xr
import rioxarray as rxr
import pandas as pd

filenames = glob.glob('*.tif')
print(filenames)
transformed_filenames = [filename.replace("-", "") for filename in filenames]
print(transformed_filenames)
def time_index_from_filenames(filenames):
    '''helper function to create a pandas DatetimeIndex
       Filename example: 20150520_0164.tif'''
    return pd.DatetimeIndex([pd.Timestamp(f[:8]) for f in filenames])


time = xr.Variable('time', time_index_from_filenames(transformed_filenames))
chunks = {'x': 1529, 'y': 2777, 'band': 1}
da = xr.concat([rxr.open_rasterio(f) for f in filenames], dim=time)

da = da.drop_vars('band') 
da.to_netcdf('OC3_AR.nc')

问题根源 #

即使投影一致，部分TIF文件的空间坐标仍可能存在差异：比如y轴像素排列方向相反（从上到下/从下到上）、x/y轴分辨率细微偏差，或者坐标范围不完全匹配。xr.concat在遇到这种情况时会自动扩展维度以兼容所有数据，最终导致尺寸翻倍。

修复步骤 #

1. 对齐所有TIF的空间维度 #

以第一个文件为空间模板，将所有后续文件的坐标、分辨率和尺寸对齐到该模板，确保合并时维度不会被扩展。修改后的代码如下：

import glob
import xarray as xr
import rioxarray as rxr
import pandas as pd

filenames = glob.glob('*.tif')
print(filenames)
transformed_filenames = [filename.replace("-", "") for filename in filenames]
print(transformed_filenames)

def time_index_from_filenames(filenames):
    '''helper function to create a pandas DatetimeIndex
       Filename example: 20150520_0164.tif'''
    return pd.DatetimeIndex([pd.Timestamp(f[:8]) for f in filenames])

# 读取第一个文件作为空间模板，确保后续文件对齐到该模板的空间参数
template_ds = rxr.open_rasterio(filenames[0])
time = xr.Variable('time', time_index_from_filenames(transformed_filenames))

# 逐个读取并对齐文件
aligned_datasets = []
for f in filenames:
    ds = rxr.open_rasterio(f)
    # 对齐到模板的空间维度（坐标、分辨率、尺寸）
    ds_aligned = ds.rio.reproject_match(template_ds)
    aligned_datasets.append(ds_aligned)

# 合并对齐后的数据集
da = xr.concat(aligned_datasets, dim=time)
da = da.drop_vars('band') 
da.to_netcdf('OC3_AR.nc')

2. 提前检查空间参数（可选） #

如果想定位具体是哪个文件导致的问题，可以在合并前验证所有文件的空间参数：

for f in filenames:
    ds = rxr.open_rasterio(f)
    print(f"文件: {f}")
    print(f"x轴范围: {ds.x.min().values:.4f} ~ {ds.x.max().values:.4f}, 分辨率: {ds.x[1]-ds.x[0]:.4f}")
    print(f"y轴范围: {ds.y.min().values:.4f} ~ {ds.y.max().values:.4f}, 分辨率: {ds.y[1]-ds.y[0]:.4f}")
    print(f"图像尺寸: {ds.shape}\n")

如果发现某文件的y轴分辨率为负数（说明像素方向与其他文件相反），可以单独反转y轴：ds = ds.reindex(y=ds.y[::-1])，再进行对齐合并。

内容的提问来源于stack exchange，提问作者Aurélien Lengrand