The following full text code and data have been released to the Hejing community . Copy the link below or read the original text to go to it, and you can fork it with one click:

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

The pcp_combine in MET can be used to calculate the cumulative precipitation. This program has a small bug (the test version is V5.2, the new version does not know whether it has been fixed or not), and an error will be reported when the forecast time exceeds 100.

pcp_combine reads two GRIB files, calculates the difference of total precipitation, and then outputs it to the netcdf file, referring to the attributes in the output file, we can use python to achieve the same function.

Note: The WRF data tested here is the Lambert projection. If it is changed to other projections, the file properties need to be adjusted accordingly.

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()



Welcome to follow, forward and contribute, Xiaobian WeChat QXBWL_001

To join the exchange group, please note your name + unit.

picture

picture