OCO2 netCDF数据Python处理求助：裁剪、滤波与CSV导出问题

依赖库安装 #

先确保安装必要的处理库：

pip install xarray netCDF4 geopandas rioxarray pandas numpy

完整代码实现 #

import glob
import re
import pandas as pd
import xarray as xr
import geopandas as gpd

# 1. 路径配置
nc_file_pattern = "oco2_LtCO2_*.nc4"  # 你的NC文件匹配路径/模式
shapefile_path = "your_shapefile.shp"  # 替换为你的目标shapefile路径
output_csv_path = "oco2_xco2_monthly_avg.csv"

# 2. 加载并预处理Shapefile
shapefile = gpd.read_file(shapefile_path)
# 统一坐标系为WGS84（OCO2数据默认坐标系）
if shapefile.crs != "EPSG:4326":
    shapefile = shapefile.to_crs("EPSG:4326")

# 3. 批量处理NC文件
daily_avg_records = []
nc_files = glob.glob(nc_file_pattern)

for file in nc_files:
    # 从文件名提取日期（匹配yymmdd格式）
    date_match = re.search(r"oco2_LtCO2_(\d{6})_", file)
    if not date_match:
        print(f"跳过无效文件名：{file}")
        continue
    date_str = date_match.group(1)
    date = pd.to_datetime(date_str, format="%y%m%d")
    
    # 读取NC文件
    try:
        ds = xr.open_dataset(file)
    except Exception as e:
        print(f"读取失败 {file}: {str(e)}")
        continue
    
    # 4. 质量标记滤波：仅保留标记为0的有效数据
    if "xco2_quality_flag" not in ds:
        print(f"文件 {file} 无质量标记变量，跳过")
        continue
    ds_filtered = ds.where(ds.xco2_quality_flag == 0, drop=True)
    if ds_filtered.xco2.size == 0:
        print(f"文件 {file} 过滤后无有效数据，跳过")
        continue
    
    # 5. Shapefile区域裁剪
    ds_filtered = ds_filtered.rio.set_crs("EPSG:4326")
    try:
        ds_clipped = ds_filtered.rio.clip(shapefile.geometry, shapefile.crs, drop=True)
    except Exception as e:
        print(f"裁剪失败 {file}: {str(e)}")
        continue
    if ds_clipped.xco2.size == 0:
        print(f"文件 {file} 裁剪后无有效数据，跳过")
        continue
    
    # 6. 计算当日xco2平均值
    daily_avg = ds_clipped.xco2.mean().item()
    daily_avg_records.append({"date": date, "xco2_daily_avg": daily_avg})
    
    # 释放内存
    ds.close()

# 7. 合并日数据并计算月平均
if not daily_avg_records:
    print("无有效数据可处理")
else:
    daily_df = pd.DataFrame(daily_avg_records)
    # 按自然月分组计算均值，同时统计当月有效天数
    monthly_avg_df = daily_df.resample("M", on="date").agg(
        xco2_monthly_avg=("xco2_daily_avg", "mean"),
        valid_days=("xco2_daily_avg", "count")
    ).reset_index()
    # 格式化年月列
    monthly_avg_df["month"] = monthly_avg_df["date"].dt.strftime("%Y-%m")
    
    # 8. 导出CSV
    monthly_avg_df.to_csv(output_csv_path, index=False)
    print(f"结果已导出至 {output_csv_path}")

关键注意事项 #

• 质量标记变量：确认你的NC文件中质量标记的变量名是否为xco2_quality_flag，若为其他名称（如quality_flag），需对应修改代码。
• 坐标系匹配：若你的OCO2数据坐标系不是WGS84，需调整rio.set_crs的参数值。
• 内存优化：处理大体积数据时，单文件处理后立即关闭数据集，避免内存溢出；若单文件仍过大，可添加chunks参数分块读取。
• 文件名适配：如果你的文件名格式有变动，需调整提取日期的正则表达式。

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