Python：TROPOMI CH4甲烷数据NC格式批量转换tif

相较于先前常规的读取nc数据方法不太一样，因为TROPOMI又进行了分组，经纬度以及柱浓度数据都在PRODUCT文件夹里面。

我的方法是第一步先将NC文件转成tif数据格式，但是直接转会有无效值以及空值

空值为：

  :_FillValue = 9.96921E36f; // float

故第二部在转换后还需将tif格式的文件中的空值过滤。完整代码如下：

import netCDF4 as nc
import numpy as np
import rasterio
from rasterio.transform import from_origin
import os
from glob import glob

# 输入目录路径
input_directory = '输入文件路径'
# 输出目录路径
output_directory = '输出文件路径'

# 确保输出目录存在
os.makedirs(output_directory, exist_ok=True)

# 获取所有NetCDF文件路径
nc_files = glob(os.path.join(input_directory, '*.nc'))

for file_path in nc_files:
    # 获取文件名（不包含扩展名）
    file_name = os.path.splitext(os.path.basename(file_path))[0]

    # 读取NetCDF文件
    dataset = nc.Dataset(file_path)

    # 提取PRODUCT组中的变量
    product = dataset['PRODUCT']
    methane_mixing_ratio = product.variables['methane_mixing_ratio'][:]
    lat = product.variables['latitude'][:]
    lon = product.variables['longitude'][:]

    # 获取无效值的填充值
    fill_value = product.variables['methane_mixing_ratio']._FillValue

    # 重新调整为一维数组
    valid_lat = lat.flatten()
    valid_lon = lon.flatten()
    valid_ch4 = methane_mixing_ratio.flatten()

    # 设置栅格的分辨率
    resolution = 0.1  # 可以根据需要调整分辨率

    # 计算栅格的边界
    min_lon, max_lon = np.min(valid_lon), np.max(valid_lon)
    min_lat, max_lat = np.min(valid_lat), np.max(valid_lat)

    # 计算栅格的大小
    width = int((max_lon - min_lon) / resolution) + 1
    height = int((max_lat - min_lat) / resolution) + 1

    # 创建空栅格
    raster_data = np.full((height, width), np.nan)

    # 计算每个点的栅格索引
    x_indices = ((valid_lon - min_lon) / resolution).astype(int)
    y_indices = ((max_lat - valid_lat) / resolution).astype(int)

    # 将值填入栅格
    raster_data[y_indices, x_indices] = valid_ch4

    # 定义栅格的仿射变换
    transform = from_origin(min_lon, max_lat, resolution, resolution)

    # 输出为临时栅格文件
    temp_output_file = os.path.join(output_directory, f'temp_{file_name}.tif')
    with rasterio.open(
        temp_output_file,
        'w',
        driver='GTiff',
        height=height,
        width=width,
        count=1,
        dtype=raster_data.dtype,
        crs='EPSG:4326',
        transform=transform,
    ) as dst:
        dst.write(raster_data, 1)

    # 读取临时栅格文件并过滤填充值
    with rasterio.open(temp_output_file) as src:
        data = src.read(1)
        data[data == fill_value] = np.nan

        # 保存处理后的栅格文件
        output_file = os.path.join(output_directory, f'output_{file_name}.tif')
        with rasterio.open(
            output_file,
            'w',
            driver='GTiff',
            height=src.height,
            width=src.width,
            count=1,
            dtype=data.dtype,
            crs=src.crs,
            transform=src.transform,
        ) as dst:
            dst.write(data, 1)

    # 删除临时文件
    os.remove(temp_output_file)

    print(f'Filtered raster file saved as {output_file}')

print('All files processed.')

PS：如果有（肯定有）直接在NC直接过滤空值的方法，请@我，谢谢。。。