利用python实现pcp_combine计算累计降水的功能

以下全文代码和数据均已发布至和鲸社区，复制下面链接或者阅读原文前往，可一键fork跑通：

https://www.heywhale.com/mw/project/6288e81ab5f5d68e30a3315d

MET中的pcp_combine可以用于计算累计降水，这个程序有一个小BUG（测试版本为V5.2，新版本不知道有没有修复），当预报时效超过100时会报错。

pcp_combine读入两个GRIB文件，计算总降水的差，再输出到netcdf文件，参考输出文件中的属性，我们完全可以用python来实现相同的功能。

注意：这里测试的WRF数据为兰伯特投影，如果换成其他投影，文件属性需要作出相应的调整。

from eccodes import *import xarray as xr from datetime import timezoneimport netCDF4 as ncimport numpy as np
fn1 = "FCST_d01_2022-05-09_06_00_00"ds1 = xr.open_dataset(fn1, engine='cfgrib', backend_kwargs={'filter_by_keys': {'typeOfLevel': 'surface',}, 'indexpath':""}) 
fn2 = "FCST_d01_2022-05-09_12_00_00"ds2 = xr.open_dataset(fn2, engine='cfgrib', backend_kwargs={'filter_by_keys': {'typeOfLevel': 'surface',}, 'indexpath':""}) 
nx = ds2.tp.GRIB_Nx ny = ds2.tp.GRIB_Nyfhr1 = 6fhr2 = 12acc  = fhr2 - fhr1
outfile = "test_out.nc"fw=nc.Dataset(outfile, 'w', format="NETCDF3_CLASSIC")
fw.MET_version = "V5.2"fw.Projection  = ds2.tp.GRIB_gridDefinitionDescriptionfw.scale_lat_1 = ds2.tp.GRIB_Latin1InDegreesfw.scale_lat_2 = ds2.tp.GRIB_Latin2InDegreesfw.lat_pin     = ds2.tp.GRIB_latitudeOfFirstGridPointInDegrees fw.lon_pin     = ds2.tp.GRIB_longitudeOfFirstGridPointInDegreesfw.x_pin       = ds2.tp.GRIB_latitudeOfSouthernPoleInDegrees fw.y_pin       = ds2.tp.GRIB_longitudeOfSouthernPoleInDegrees 
fw.lon_orient  = ds2.tp.GRIB_LoVInDegreesfw.d_km        = ds2.tp.GRIB_DxInMetres / 1000.fw.r_km        = 6371.200000fw.nx          = nxfw.ny          = f"{ny} grid_points" 
fw.createDimension('lat',ny)fw.createDimension('lon',nx)
lat   = fw.createVariable('lat', 'f4', ('lat','lon'))lon   = fw.createVariable('lon', 'f4', ('lat','lon'))APCP  = fw.createVariable(f"APCP_{acc:02d}","f4",("lat","lon",),fill_value=-9999)
lon.long_name = "longitude"lon.units     = "degrees_east"lon.standard_name="longitude"
lat.long_name ="latitude"lat.units     = "degrees_north"lat.standard_name="latitude"
init  = ds2.time.data.astype('datetime64[s]').tolist().replace(tzinfo=timezone.utc)valid = ds2.valid_time.data.astype('datetime64[s]').tolist().replace(tzinfo=timezone.utc)
APCP.setncattr("name", f"APCP_{acc:02d}") #保留关键字，必须要setncattr设置APCP.long_name = "Total precipitation"APCP.level     = f"A{acc}"APCP.units     = "kg/m^2"APCP.init_time     = init.strftime("%Y%m%d_%H%M%S")APCP.init_time_ut  = int(init.timestamp())APCP.valid_time    = valid.strftime("%Y%m%d_%H%M%S")APCP.valid_time_ut = int(valid.timestamp())APCP.accum_time    = f"{acc:02}0000"APCP.accum_time_sec= int(acc*3600
lat[:] = ds2.latitude.data[:]lon[:] = ds2.longitude.data[:]APCP[:] = ds2.tp.data[:] - ds1.tp.data[:]fw.close()