Python转nc/hdf文件为栅格TIFF格式-CMIP6常见错误

背景

在做气候变化相关研究的时候，从官网下载了CMIP6数据，发现大多为NetCDF (nc)或hdf格式。为了把它的每个波段数据批量转换为我所熟悉的栅格TIFF格式，笔者首先通过网络和ChatGPT搜索了基于Python的处理方法。检查时发现，这些代码考虑得都不够周到，甚至会出现生成的栅格上下颠倒的情况。

中间的空白竖线并非数据缺失，该数据的经度范围为0°~360°，0°左侧小部分西经地区在数值上经度为负数，未能显示。

查阅资料后，发现网络上的代码在使用gdal库的SetGeoTransform方法设置地图参数时，没有考虑数据的实际情况。

# 有缺陷的代码
# 代码前后部分略去，lons, lats, lons, lats为从原nc/hdf数据中提取的经纬度序列，data为二维气象数据
driver = gdal.GetDriverByName('GTiff')
ds_out = driver.Create(output_path, data.shape[1], data.shape[0], 1, gdal.GDT_Float32)
ds_out.SetProjection(srs.ExportToWkt())
# 经纬度范围
[lonMin, latMax, lonMax, latMin] = [min(lons), max(lats), max(lons), min(lats)]
N_lat = len(lats)
N_lon = len(lons)
# 分辨率
lon_resolution = (lonMax - lonMin) / (N_lon - 1)
lat_resolution = (latMax - latMin) / (N_lat - 1)
# Gdal设置影像的显示范围
geotransform = [lonMin, lon_resolution, 0, latMax, 0, -lat_resolution]
ds_out.SetGeoTransform(geotransform)
# 将数据写入 tif 文件
out_band = ds_out.GetRasterBand(1)
out_band.WriteArray(data)

这种方法设置geotransform参数时，有两个没能考虑到的问题：

1. gdal的官网确实提到，第1和第4个参数应当是数据范围的左上角。但这只限于气象数据是从左上角开始排列的情况下。展示lons和lats数据可以发现，CMIP6气象数据应该是从左下角作为起点开始排列，如果把起始点设置为左上角，会出现上下颠倒的情况
2. 第1和第4个参数指的是角落栅格的顶点坐标，查阅SSP的元数据可知，数据中的经纬度应该是栅格的中心点（也符合像元居中配准的习惯），因此第1和第四个参数要额外减去半个分辨率长度

发现这些问题后，我对代码进行了调整，不管气象数据是从地图哪个角落作为起点排列，都能正常输出结果。
首先通过file_info函数查看数据信息，从而修改nc2tiff函数里的字段名

nc文件转tiff栅格数据

以CMIP6的历史日降雨数据为例，将nc文件里的每个波段分别转换为tiff格式数据，并保存于指定文件夹内。

from netCDF4 import Dataset
import xarray as xr
from osgeo import gdal, osr
# filepath: nc数据在电脑里的保存路径
# output_folder：转换后的tiff格式文件
def file_info(filepath):
    # 查看nc数据的情况
    # 打开nc文件
    nc_obj = Dataset(filepath)
    # 查看nc文件有哪些变量
    print(nc_obj.variables.keys())


def nc2tiff(filepath, output_folder):
    # 打开 nc 文件
    ds = xr.open_dataset(filepath)

    # 把pr变量中的数据抽取出来
    pr = ds['pr']

    # 获取nc文件中的lon和lat信息
    lons = ds['lon'].values
    lats = ds['lat'].values

    # 设置 GeoTIFF 的 spatial reference
    srs = osr.SpatialReference()
    srs.ImportFromEPSG(4326)  # WGS84

    # 遍历所有的时间点
    for i in range(len(ds['time'])):
        # 把当前时间点的变量数据抽取出来
        data = pr[i].values

        # 创建 GeoTIFF 文件
        driver = gdal.GetDriverByName('GTiff')

        # 生成文件名
        t_str = str(ds['time'][i].values)  # e.g. 1950-01-01 12:00:00
        t_str = t_str.replace(':', '-')  # 之后需要用时间命名生成的tiff文件，但是':'不能作为文件名
        print(t_str)
        output_path = output_folder + '\\' + '{}.tiff'.format(t_str)

        ds_out = driver.Create(output_path, data.shape[1], data.shape[0], 1, gdal.GDT_Float32)

        # 设置文件的 spatial reference 和 geotransform
        ds_out.SetProjection(srs.ExportToWkt())

        lon_resolution = (lons[-1] - lons[0]) / (len(lons) - 1)
        lat_resolution = (lats[-1] - lats[0]) / (len(lats) - 1)

        # geotransform参数设置
        # 0：起点X坐标; 1：X向分辨率; 2：X向仿射变换角度; 3：起点Y坐标; 4：Y向仿射变换角度; 5：Y向分辨率
        # 注意1：官方文档称geotransform参数应该以左上角为基准，但查阅lons和lats数据后发现pr数据的顺序是从地图左下角开始排列的，
        # 如果盲目输入左上角坐标，写入pr数据后得到的tif文件会上下颠倒
        # 注意2：我们的数据应该是栅格中心的坐标，因此地图顶点的坐标应该要偏移半个像元尺寸（像元中心配准）
        geotransform = [lons[0] - lon_resolution / 2, lon_resolution, 0, lats[0] - lat_resolution / 2, 0,
                        lat_resolution]  # 不管降雨数据的起点在地图哪个角，都可以通用
        print(geotransform)
        ds_out.SetGeoTransform(geotransform)

        # 将数据写入 tif 文件
        out_band = ds_out.GetRasterBand(1)
        out_band.WriteArray(data)

        # 保存文件
        ds_out.FlushCache()

hdf数据转tiff栅格

以CMIP6的SSP126日降雨数据为例，将hdf文件里的每个波段分别转换为tiff格式数据，并保存于指定文件夹内。首先通过file_info查看数据字段，再以此调整hdf2tiff函数

import h5py
from osgeo import gdal, osr
import os

def file_info(filepath):
    # 查看hdf文件的关键词信息
    # 打开HDF5文件
    with h5py.File(filepath, 'r') as f:

        # 列出顶层组
        print("Keys: %s" % f.keys())
        top_item = list(f.keys())[0]  # 如果有多个项，可能需要调整索引

        # 获取顶层项
        item = f[top_item]

        # 检查这个项是否是数据集
        if isinstance(item, h5py.Dataset):
            # 如果是数据集，直接打印形状、类型和部分数据
            print("Data shape : ", item.shape)
            print("Data type: ", item.dtype)
            # 通过读取组和其子项目
            print("Group keys: ", item.keys())


def hdf2tiff(dhf_filepath, output_folder):
    # 打开你的 HDF5 文件
    file = h5py.File(dhf_filepath, 'r')

    # 输出 TIFF 文件的目录
    print(output_folder)

    # 获取时间数据集
    time_dataset = file['time']
    time_data = time_dataset[:]

    # 获取降雨数据集
    pr_dataset = file['pr']
    lon_dataset = file['lon']
    lat_dataset = file['lat']

    # 遍历所有时间点，将每个时间点的降雨数据转换为一个独立的 TIFF 文件
    for i, time in enumerate(time_data):
        # 获取该时间点的降雨数据
        pr_data = pr_dataset[i]

        # 获取数据集的形状
        y_size, x_size = pr_data.shape  # 降雨数据的维度，即列数和行数

        # 创建一个新的 TIFF 文件
        output_filename = os.path.join(output_folder, f'{time}.tiff')
        print(output_filename)
        driver = gdal.GetDriverByName('GTiff')
        dataset_output = driver.Create(output_filename, x_size, y_size, 1, gdal.GDT_Float32)

        # 设定地理坐标（写入降雨数据的地理起点与方向）
        # # 0：起点X坐标; 1：X向分辨率; 2：X向仿射变换角度; 3：起点Y坐标; 4：Y向仿射变换角度; 5：Y向分辨率
        # # 注意1：官方文档称geotransform参数应该以左上角为基准，但查阅lon_data和lats数据后发现pr数据的顺序是从地图左下角开始排列的，
        # #         如果盲目根据左上角坐标设定参数，写入pr数据后得到的tif文件会上下颠倒
        # # 注意2：我们的数据应该是栅格中心的坐标，因此地图顶点的坐标应该要偏移半个像元尺寸（像元中心配准）
        lon_resolution = (lon_dataset[-1] - lon_dataset[0]) / (len(lon_dataset) - 1)
        lat_resolution = (lat_dataset[-1] - lat_dataset[0]) / (len(lat_dataset) - 1)
        geotransform_para = [lon_dataset[0] - lon_resolution / 2, lon_resolution, 0, lat_dataset[0] - lat_resolution / 2,
                             0, lat_resolution]  # 不管降雨数据的起点在地图哪个角，都可以通用
        dataset_output.SetGeoTransform(geotransform_para)

        # 设定投影
        srs = osr.SpatialReference()
        srs.ImportFromEPSG(4326)  # 这是 WSG84 地理坐标系统的EPSG编号
        dataset_output.SetProjection(srs.ExportToWkt())

        # 写入数据
        dataset_output.GetRasterBand(1).WriteArray(pr_data)

        # 保存文件
        dataset_output.FlushCache()

总结

在不熟悉的领域开展分析时，还是需要耐心学习相关库的底层逻辑，不能盲目照搬。由于这个错误很难发现，笔者决定在网络上分享一下自己的思考和经验。由于笔者并非GIS科班出身，难免有错误之处，欢迎指出。

ZJU-安中DZH