裁剪netCDF文件无数据求助：NOAA降水数据转美国州级分区

我正在开展一项概念验证工作，目标是将NOAA提供的含1988年全球降水数据的netCDF文件裁剪至美国州级尺度，但执行裁剪操作后始终无数据输出。

使用资源：

• 美国州级地理数据：美国人口普查局2018年发布的州级shapefile
• netCDF降水数据：NOAA CPC全球降水1988年数据集

我日常主要使用C/C++/C#，并非Python用户，当前运行环境为Windows。以下是我的代码：

#!/usr/bin/python

import xarray
import geopandas
import glob
import os, sys, getopt
from shapely.geometry import mapping

def partition(shapeFile, ncDirectory, outDirectory, stateName):
    # Load the shape file that either has all the states or just one named state
    sf_all = geopandas.read_file(shapeFile)

    # Filter to our state
    sf = sf_all.loc[sf_all['NAME'] == stateName]

    # Get all the nc files in the ncDirectory so we can repartition all the models
    file_list = sorted(glob.glob(f'{ncDirectory}/*.nc'))

    for file in file_list:
        print(f'Working on: {outDirectory}/{stateName}/{os.path.basename(file)}')

        # Open the file and do some basics
        nc_file = xarray.open_dataset(file)
        nc_file.rio.write_crs("epsg:4326", inplace=True)
        nc_file.rio.set_spatial_dims(x_dim="lon", y_dim="lat", inplace=True)

        # Convert lon from 0,360 to -180,180 if needed
        if any(nc_file.coords['lon'] > 180 ):
            nc_file.coords['lon'] = (nc_file.coords['lon'] + 180) % 360 - 180
        
        # Convert lat from 0,360 to -180,180 if needed
        if any(nc_file.coords['lat'] > 180 ):
            nc_file.coords['lat'] = (nc_file.coords['lat'] + 180) % 360 - 180

        # Perform the actual clipping here based on the geometry of the {sf} variable which should just be our state
        clipped_nc = nc_file.rio.clip(sf.geometry.apply(mapping), crs="epsg:4326", drop=True)

        if os.path.exists(f'{outDirectory}/{stateName}/{os.path.basename(file)}'):
            os.removedirs(f'{outDirectory}/{stateName}/{os.path.basename(file)}')
        
        # Save the file to a new netCDF representing the state/model combination
        clipped_nc.to_netcdf(f'{outDirectory}/{stateName}/{os.path.basename(file)}')

def main(argv):
    shapeFile = ''
    stateName = ''
    ncDirectory = ''
    outDirectory = ''

    opt, args = getopt.getopt(argv,"hf:n:o:s:",["shapedir=","ncdirectory=","outputdir=", "statename="])
    for opt, arg in opt:
        if opt =='-h':
            print('repartition.py -f <shapefile> -n <ncfiledirectory> -o <outputdir> -s <statename>')
            quit(0)
        elif opt in ("-f", "--shapefile"):
            shapeFile = arg.rstrip("/")
        elif opt in ("-n", "--ncdirectory"):
            ncDirectory = arg.rstrip("/")
        elif opt in ("-o", "--outputdir"):
            outDirectory = arg.rstrip("/")
        elif opt in ("-s","--statenames"):
            stateName = arg
        
    isExists = os.path.exists(f'{outDirectory}/{stateName}')
    if not isExists:
        os.makedirs(f'{outDirectory}/{stateName}')
    
    partition(shapeFile,ncDirectory,outDirectory, stateName)

if __name__ == "__main__":
    main(sys.argv[1:])

裁剪无数据输出的核心原因包括：经度转换后未保持单调递增、CRS匹配逻辑不严谨、文件删除函数误用、冗余的纬度转换代码干扰。以下是修复后的代码及关键说明：

修复后的代码 #

#!/usr/bin/python

import xarray as xr
import geopandas as gpd
import glob
import os, sys, getopt

def partition(shapeFile, ncDirectory, outDirectory, stateName):
    # 加载shapefile并筛选目标州
    sf_all = gpd.read_file(shapeFile)
    sf = sf_all.loc[sf_all['NAME'] == stateName]
    
    # 确保shapefile与netCDF数据CRS一致（EPSG:4326）
    if sf.crs != "EPSG:4326":
        sf = sf.to_crs("EPSG:4326")

    file_list = sorted(glob.glob(f'{ncDirectory}/*.nc'))

    for file in file_list:
        out_path = f'{outDirectory}/{stateName}/{os.path.basename(file)}'
        print(f'Working on: {out_path}')

        # 打开netCDF文件并设置空间属性
        nc_file = xr.open_dataset(file)
        nc_file = nc_file.rio.set_spatial_dims(x_dim="lon", y_dim="lat").rio.write_crs("EPSG:4326")

        # 经度从0-360转换为-180-180，并保持坐标单调递增
        if nc_file.lon.max() > 180:
            nc_file = nc_file.assign_coords(lon=(nc_file.lon + 180) % 360 - 180)
            nc_file = nc_file.sortby('lon')

        # 执行裁剪（直接传入GeoDataFrame的geometry和CRS）
        clipped_nc = nc_file.rio.clip(sf.geometry, sf.crs, drop=True)

        # 删除已存在的输出文件（用os.remove而非os.removedirs）
        if os.path.exists(out_path):
            os.remove(out_path)
        
        # 保存裁剪结果
        clipped_nc.to_netcdf(out_path)
        print(f'Saved: {out_path}')

def main(argv):
    shapeFile = ''
    stateName = ''
    ncDirectory = ''
    outDirectory = ''

    opt, args = getopt.getopt(argv,"hf:n:o:s:",["shapefile=","ncdirectory=","outputdir=", "statename="])
    for opt, arg in opt:
        if opt =='-h':
            print('repartition.py -f <shapefile> -n <ncfiledirectory> -o <outputdir> -s <statename>')
            sys.exit(0)
        elif opt in ("-f", "--shapefile"):
            shapeFile = arg.rstrip("/")
        elif opt in ("-n", "--ncdirectory"):
            ncDirectory = arg.rstrip("/")
        elif opt in ("-o", "--outputdir"):
            outDirectory = arg.rstrip("/")
        elif opt in ("-s","--statename"):
            stateName = arg
        
    # 创建州级输出目录
    state_out_dir = f'{outDirectory}/{stateName}'
    if not os.path.exists(state_out_dir):
        os.makedirs(state_out_dir)
    
    partition(shapeFile, ncDirectory, outDirectory, stateName)

if __name__ == "__main__":
    main(sys.argv[1:])

关键修改说明 #

1. 修复经度排序问题：0-360转-180-180后，调用sortby('lon')确保经度坐标单调递增，这是rioxarray正确识别空间范围的前提。
2. 移除冗余纬度转换：NOAA降水数据纬度范围为-90到90，无需转换，删除该代码避免干扰。
3. 规范CRS匹配：检查并转换shapefile到EPSG:4326，确保与netCDF数据的坐标系统完全一致。
4. 简化裁剪调用：直接传入GeoDataFrame的geometry和CRS，无需手动映射，rioxarray可直接处理GeoDataFrame。
5. 修复文件删除逻辑：用os.remove删除已存在的文件，os.removedirs仅用于删除空目录，原代码会导致删除文件时报错。

额外排查步骤 #

若仍无数据输出，可执行以下检查：

• 打印sf确认是否正确筛选到目标州（注意州名大小写敏感，如"California"而非"california"）。
• 打印nc_file.lon.min()/max()和nc_file.lat.min()/max()，确认数据范围覆盖目标州。
• 可视化netCDF数据与shapefile的空间范围，确认两者存在重叠区域。

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