从NC到TIFF：Python处理ERA5-Land数据的完整地理空间工作流

气象数据科学家和GIS分析师经常需要处理来自欧洲中期天气预报中心（ECMWF）的ERA5-Land数据集。这些数据通常以NetCDF格式存储，包含丰富的时空信息，但要在GIS软件中使用，往往需要转换为更通用的GeoTIFF格式。本文将带你完成从数据加载、时间维度聚合到地理空间转换的完整流程。

1. 理解ERA5-Land数据的基本结构

ERA5-Land是ECMWF提供的高分辨率地表再分析数据集，时间分辨率可达小时级别。每个NetCDF文件通常包含多个变量（如2米气温、降水等）、时间维度、经度维度和纬度维度。

典型的ERA5-Land数据结构可以通过xarray查看：

import xarray as xr

# 加载示例数据
ds = xr.open_dataset('era5_land_sample.nc')
print(ds)

输出会显示类似这样的结构：

<xarray.Dataset>
Dimensions:    (time: 24, latitude: 721, longitude: 1440)
Coordinates:
  * time       (time) datetime64[ns] 2020-01-01 ... 2020-01-01T23:00:00
  * latitude   (latitude) float64 90.0 89.75 89.5 ... -89.75 -90.0
  * longitude  (longitude) float64 0.0 0.25 0.5 0.75 ... 359.25 359.5 359.75
Data variables:
    t2m        (time, latitude, longitude) float32 ...

关键点：

• 时间维度 ：通常为UTC时间
• 空间分辨率 ：0.1度或0.25度
• 坐标顺序 ：注意latitude可能是从北到南排列

2. 时间维度聚合：从小时到日均值

原始的小时数据对于许多应用来说过于精细，计算日均值是常见需求。xarray提供了简便的方法进行时间维度的聚合：

# 计算日均值
daily_mean = ds['t2m'].resample(time='1D').mean()

# 或者更精确的24小时平均
daily_mean_24h = ds['t2m'].groupby('time.date').mean()

注意：ERA5-Land的小时数据通常从00:00到23:00，共24个时间点。确保你的聚合方法符合分析需求。

对于批量处理多个文件，可以构建如下工作流：

import os
import xarray as xr

input_folder = '/path/to/hourly_data'
output_folder = '/path/to/daily_data'

for filename in os.listdir(input_folder):
    if filename.endswith('.nc'):
        filepath = os.path.join(input_folder, filename)
        
        # 打开数据集并计算日均值
        with xr.open_dataset(filepath) as ds:
            daily_mean = ds['t2m'].resample(time='1D').mean()
            
            # 保存结果
            output_path = os.path.join(output_folder, f'daily_{filename}')
            daily_mean.to_netcdf(output_path)

3. 地理空间转换：从NetCDF到GeoTIFF

将处理后的数据转换为GeoTIFF格式需要精确的空间参考信息。rasterio库提供了强大的栅格数据处理能力：

3.1 理解地理变换(transform)

地理变换定义了栅格数据中像素坐标与地理坐标的映射关系。关键参数包括：

• 左上角坐标
• 像素宽度和高度
• 可能的旋转（通常为0）

创建transform的典型方法：

from rasterio.transform import from_origin

# 假设已知左上角坐标和分辨率
lon_min = -180.0
lat_max = 90.0
resolution = 0.1  # 度

transform = from_origin(lon_min, lat_max, resolution, resolution)

3.2 完整的NetCDF到GeoTIFF转换函数

下面是一个完整的转换函数，考虑了坐标参考系统(CRS)：

import rasterio
from rasterio.transform import from_origin
import numpy as np

def netcdf_to_geotiff(nc_file, variable_name, output_path, crs='EPSG:4326'):
    """
    将NetCDF文件中的变量转换为GeoTIFF格式
    
    参数:
        nc_file: NetCDF文件路径
        variable_name: 要转换的变量名
        output_path: 输出TIFF文件路径
        crs: 坐标参考系统，默认为WGS84
    """
    # 打开NetCDF文件
    with xr.open_dataset(nc_file) as ds:
        data = ds[variable_name]
        
        # 如果是多维数据，取第一个时间点
        if 'time' in data.dims:
            data = data.isel(time=0)
        
        # 获取地理信息
        lons = ds['longitude'].values
        lats = ds['latitude'].values
        lon_min, lon_max = lons.min(), lons.max()
        lat_min, lat_max = lats.min(), lats.max()
        
        # 计算分辨率
        x_res = (lon_max - lon_min) / (len(lons) - 1)
        y_res = (lat_max - lat_min) / (len(lats) - 1)
        
        # 创建transform
        transform = from_origin(lon_min, lat_max, x_res, -y_res)
        
        # 写入GeoTIFF
        with rasterio.open(
            output_path,
            'w',
            driver='GTiff',
            height=data.shape[0],
            width=data.shape[1],
            count=1,
            dtype=data.dtype,
            crs=crs,
            transform=transform
        ) as dst:
            dst.write(data.values, 1)

3.3 批量处理脚本

结合时间聚合和格式转换，下面是完整的批量处理脚本：

import os
import xarray as xr
import rasterio
from rasterio.transform import from_origin

def process_era5_land(input_folder, output_folder, variable='t2m'):
    """批量处理ERA5-Land数据：时间聚合+格式转换"""
    
    # 确保输出文件夹存在
    os.makedirs(output_folder, exist_ok=True)
    
    for filename in os.listdir(input_folder):
        if filename.endswith('.nc'):
            input_path = os.path.join(input_folder, filename)
            
            # 处理文件名
            output_filename = filename.replace('.nc', '.tif')
            output_path = os.path.join(output_folder, output_filename)
            
            # 1. 打开数据并计算日均值
            with xr.open_dataset(input_path) as ds:
                # 计算日均值
                daily_mean = ds[variable].resample(time='1D').mean()
                
                # 2. 转换为GeoTIFF
                lons = ds['longitude'].values
                lats = ds['latitude'].values
                lon_min, lon_max = lons.min(), lons.max()
                lat_min, lat_max = lats.min(), lats.max()
                x_res = (lon_max - lon_min) / (len(lons) - 1)
                y_res = (lat_max - lat_min) / (len(lats) - 1)
                
                transform = from_origin(lon_min, lat_max, x_res, -y_res)
                
                with rasterio.open(
                    output_path,
                    'w',
                    driver='GTiff',
                    height=len(lats),
                    width=len(lons),
                    count=1,
                    dtype=daily_mean.dtype,
                    crs='EPSG:4326',
                    transform=transform
                ) as dst:
                    # 取第一个时间点的数据
                    dst.write(daily_mean.isel(time=0).values, 1)

4. 在GIS软件中使用结果

转换后的GeoTIFF文件可以直接在QGIS或ArcGIS等软件中使用。以下是一些常见应用场景：

4.1 QGIS中的可视化

1. 加载数据 ：直接将TIFF文件拖入QGIS工作区
2. 调整样式 ：
   • 右键图层 → 属性 → 符号化
   • 选择合适的配色方案（如热力图用于温度数据）
3. 添加底图 ：使用QuickMapServices插件添加OpenStreetMap等底图

4.2 空间分析

在GIS软件中可以执行各种空间分析操作：

• 区域统计（计算某行政区的平均温度）
• 重采样到其他分辨率
• 与其他空间数据（如土地利用、高程）进行叠加分析

4.3 常见问题排查

问题现象	可能原因	解决方案
数据位置偏移	transform设置错误	检查左上角坐标和分辨率
数据方向颠倒	纬度顺序问题	检查数据是否需翻转(numpy.flipud)
投影不正确	CRS不匹配	确保使用正确的EPSG代码

5. 高级技巧与优化建议

5.1 内存优化处理大型数据集

对于大区域或长时间序列数据，可以使用dask进行分块处理：

import dask.array as da

# 分块打开NetCDF文件
ds = xr.open_dataset('large_file.nc', chunks={'time': 10})

# 计算日均值（延迟执行）
daily_mean = ds['t2m'].resample(time='1D').mean()

# 显式计算
daily_mean.compute()

5.2 多变量处理

如果需要处理多个变量，可以修改函数增加变量选择：

def process_multiple_variables(nc_file, output_folder, variables=['t2m', 'tp']):
    """处理多个变量到单独的TIFF文件"""
    with xr.open_dataset(nc_file) as ds:
        for var in variables:
            output_path = os.path.join(output_folder, f'{var}.tif')
            # ...转换逻辑...

5.3 时间序列动画制作

将每日数据转换为时间序列动画可以直观展示变化：

import matplotlib.pyplot as plt
import imageio

# 生成每日图片
images = []
for day in range(1, 32):
    # 加载当日数据并创建图片
    # ...省略数据处理代码...
    images.append(image)
    
# 保存为GIF
imageio.mimsave('animation.gif', images, fps=5)