Python利用shp文件裁剪netcdf文件

第1步：引入库函数

import netCDF4                # 用于读取NetCDF4格式数据（气象数据常用格式）
import numpy as np            # 数值计算库，处理数组数据
from osgeo import gdal, osr, ogr  # GDAL地理数据处理库，用于栅格化矢量、坐标转换等
import matplotlib.pyplot as plt   # 绘图库，用于数据可视化

第2步：掩膜生成函数 makeMask

def makeMask(lon, lat, res, shapefile):# 打开矢量文件（Shapefile）source_ds = ogr.Open(shapefile)  # 使用GDAL打开Shapefile，失败返回Nonesource_layer = source_ds.GetLayer()  # 获取矢量文件的图层（如国界多边形）# 在内存中创建空白栅格mem_ds = gdal.GetDriverByName('MEM').Create('', lon.size, lat.size, gdal.GDT_Byte)# 参数解释：# '' : 内存栅格无文件名# lon.size : 栅格列数（经度方向）# lat.size : 栅格行数（纬度方向）# gdal.GDT_Byte : 像素类型为字节（0-255）# 设置地理变换参数（定义栅格空间位置和分辨率）mem_ds.SetGeoTransform((lon.min(), res, 0, lat.max(), 0, -res))# 参数格式：(左上角经度, 经度分辨率, 旋转项, 左上角纬度, 旋转项, 纬度分辨率)# 纬度分辨率为负值，因为纬度从上到下递减band = mem_ds.GetRasterBand(1)  # 获取栅格第一个波段（单波段掩膜）# 栅格化矢量图层（将多边形区域填充为1）gdal.RasterizeLayer(mem_ds, [1], source_layer, burn_values=[1])# [1] : 目标波段列表（此处为第一个波段）# burn_values=[1] : 将矢量覆盖区域像素值设为1array = band.ReadAsArray()  # 将栅格数据转为NumPy数组mem_ds = None  # 释放内存栅格（GDAL对象需手动清理）band = None    # 释放波段对象return array   # 返回二值掩膜数组（1=区域内，0=区域外）

第3步：主程序

# 设置Matplotlib图像分辨率
plt.figure(dpi=600)  # 创建画布并设置DPI为600（高清图像）# 定义文件路径
shapefile = "E:/BaiduNetdiskDownload/cut_nc/vn_shp/VNM_adm0.shp"  # 越南国界Shapefile
ncs = "E:/BaiduNetdiskDownload/cut_nc/ds_tmys_mmsVN_1976-2005_pr_historical_GFDL-ESM2M.nc"  # NetCDF降水数据文件# 读取NetCDF文件
nc = netCDF4.Dataset(ncs, 'r')  # 以只读模式打开NetCDF文件
pr = nc.variables['pr'][:]       # 提取降水变量数据（可能为3D：时间×纬度×经度）# 显示原始降水数据
plt.imshow(pr)  # 若pr为3D数组，此处会报错（需指定时间层，如pr[0]）
plt.show()      # 显示图像# 提取经纬度数据
lons = nc.variables['lon'][:]  # 经度数组（1D）
lats = nc.variables['lat'][:]  # 纬度数组（1D）# 计算像元尺寸（潜在问题点）
cellsize = lons[:][1] - lons[:][0]  # 等价于lons[1] - lons[0]，但写法冗余
# 正确写法：cellsize = lons[1] - lons[0]# 生成掩膜
mask = makeMask(lons, lats, cellsize, shapefile)  # 调用函数生成掩膜数组# 显示掩膜
plt.imshow(mask)  # 显示二值掩膜（黑白区域）
plt.show()# 应用掩膜到降水数据
precip = np.ma.masked_where(mask == 0, pr)  # 使用NumPy掩膜模块，将mask=0的区域设为无效
# 潜在问题：若pr为3D（时间×纬度×经度），需指定时间层，如pr[0]# 显示掩膜后的降水数据
plt.imshow(precip)
plt.show()# 输出统计指标
print(np.min(precip), np.mean(precip), np.max(precip))  # 输出最小值、均值、最大值

绘图结果如下：

注意：
（1）上面的 pr 类型是array of float类型；
（2）mask的范围是以netcdf文件的范围和坐标范围为准的；

例如在利用中国shp文件对世界范围Netcdf文件裁剪时代码修改部分如下：

# 路径设置
plt.figure(dpi=600)
shapefile = "E:/Selecting_ET_with_ETC_method/matfiles/ChinaShp/ChinaBoundary_WGS1984.shp"
ncs= "E:/DataDealing/ET_DOLCE/DOLCE_v3_1980.nc"# 读入NC文件
nc = netCDF4.Dataset(ncs,'r')# 得到降水数据
pr = nc.variables['hfls'][:].mean(axis=0) # 变量名修改，并对其时间维度进行平均
pr = np.flip(pr,0) # 需要注意范围是否需要镜像# 展示降水数据
plt.imshow(pr)
plt.show()

上面的中国例子中，shp是中国范围，但nc文件时全球的，裁剪结果是世界范围内！