wgrib2: rpn
Introduction
The rpn option runs a reverse polish notation
(RPN)
calculator. Having a builtin calculator is quite handy. We use it to convert
units (ex. geopotential to geopotential meters, accumulations to rates),
compute simple quantities (net flux from downward and upward fluxes),
and even compute the plant hardiness index from the 2 m temperatures. The goal
of the calculator is to reduce the need to write simple fortran programs that
handle grib2 files.
The "hardware" of the rpn calculator consists of 20 registers
and a stack (10 entries deep). (Wgrib2 prior to 2.0.6 has 10 registers.) Each register and stack entry is
an array which is the same size as the current grid size being
processed. So the main differences between the
rpn calculator and a storebought calculator
is that you are working with arrays and the top of the stack is initialized
with the array of grid values from the current grib message.
The 20 registers (0..19) are initialized when wgrib2 is started.
The stack is different, it gets initialized every time
the rpn option is run. First the stack is cleared
and the grid values (originally from currently processed grib message)
is pushed onto the top of the stack. At the end of the
the rpn option, the top of the stack is saved
in the grid values unless the stack has been emptied. To save the
calculations, you can save them in a register or write them out
by grib_out, bin, ieee, text, etc.
Callable wgrib2, wgrib2api
Callable Wgrib2 and wgrib2api use the rpn registers to transfer
gridded data between the calling program and the wgrib2 subroutine. A Fortran or C
program can read and write the wgrib2's RPN registers. For example, if a program
wants to write a grib2 file, it would first place the grid values into a register.
The calling program would then call the wgrib2 subroutine with instructions to
read a template, replace the grid values with the register values, and then
write a grib message.
Uses
 change of units when importing data (gribifying data)
 computations: ex, U,V > wind speed, wind direction, potential temperature
 merging data
 complex masking of data
 changing units before writing text/ieee files
 removing extreme data values
 finding min and max values
 finding the globally averaged precipitation
 comparing fields
Usage
rpn "A:B:C:..."
A,B,C,.. = number, rpn function, or rpn operator
Operators and Functions:
Pop X, Push Fn(X)
 abs: absolute value
 acos: arc cos, [0, pi] radians
 alt_x_scan: changes alternate x scanning to regular x scanning and vice versa
 asin: arc sin, result is [pi/2, pi/2] radians
 atan: arc tan, result is [pi/2, pi/2] radians, see atan2
 abs: absolute value
 ceil: smallest integer >= X
 cos: cosine, argument in radians
 exp: e^X
 floor: largest integer <= X
 ln: natural logorithm
 sin: sine, argument in radians
 smth9g: smth9 for global fields
 smth9r: snth9 for regional fields
 sq: X*X
 sqrt: square root
 tan: tangent, argument in radians
 xave: for nxny grids, averages in the x direction (normally zonal mean)
 xdev: for nxny grids, remove x average (normally deviation from zonal mean)
 yrev: for nxny grids, changes we:sn to we:ns and vice versa
 1/x: 1/X
Pop Y, Pop X, push Fn(X,Y)
 +: push X+Y
 : push XY
 *: push X*Y
 /: push X/Y
 <: push X < Y (1/0 if true/false)
 <=: push X <= Y (1/0 if true/false)
 ==: push X == Y (1/0 if true/false)
 !=: push X != Y (1/0 if true/false)
 >=: push X >= Y (1/0 if true/false)
 >: push X > Y (1/0 if true/false)
 atan2: push arctan(X/Y), result is [pi, pi] radians, see atan
 pow: push X**Y (X^Y)
 mask: if (Y != 0) push(X) else push(UNDEFINED)
 max: push max(X,Y)
 merge: if (Y != UNDEFINED) push(Y) else push(X)
 min: push min(X,Y)
Note: an operation involving an UNDEFINED is UNDEFINED
Stack Operators:
 clr, clear the stack (stack is emptied)
 dup, duplicate the top of the stack
 pop, remove the top of the stack
 exc/swap, exchange the top 2 stack entries
Register Operators: (note: CW2 v2.0.6+ uses registers 7,8,9 prior versions 0,1,2)
 clr_I, clear register I, I=0,1..,9 (19 for v2.0.6+)
 rcl_I, push register I on top of stack, I=0,1..,9 (19 for v2.0.6+)
 sto_I, save top of stack in register I, I=0,1..,9 (19 for v2.0.6+)
 rcl_lat, push latitudes onto the top of the stack (degrees)
 rcl_lon, push longitudes onto the top of the stack (degrees)
Variables and Constants: put on the top of the stack
 number number = floating point or integer number like 0, 10.1, 1.23e4
 days_in_ref_month number of days in the month for the reference date (conversion between monthly acc. and rates)
 days_in_verf_month number of days in the month for the verification time (conversion betwee monthly acc. and rates)
 pi 3.1415....
 rand random number uniformly distributed between 0 and 1, each grid point has a different random number
Printing Operators:
 print_corr, write cosine weighted spatial correlation, data[TOP] data[TOP1]
 print_max, print_min, data[TOP]
 print_rms, write cosine weighted RMS, data[TOP]data[TOP1]
 print_diff, write cosine weighted difference, data[TOP1]data[TOP]
 print_ave, write cosine weighted average, grid_ave(data[TOP]*cos(grid))/grid_ave(cos(grid))
 print_wt_ave, write weighted average grid_ave, data[TOP] is the weighting
print grid_ave(data[TOP]*data[TOP1])/grid_ave(data[TOP])
Examples
The standard units of grib temperature is K but you want the text output in Celcius.
$ wgrib2 a.grb match ":TMP:850 mb:" rpn "273.15:" text C.dat
Fahrenheit is easy too (F = (K273.15)*9/5+32).
$ wgrib2 a.grb match ":TMP:850 mb:" rpn "273.15::9:*:5:/:32:+" text F.dat
Suppose you want to limit the relative humidity values to 100. This example only
affect the RH fields. All submessages will be converted into messages.
$ wgrib2 a.grb if ":RH:" rpn "100:min" fi grib_out out.grb not_if ":RH:" grib out.grb
Write out the 500 mb wind speed.
$ wgrib2 a.grb match ":[UV]grd:500 mb:" \
if ":UGRD:" rpn "sto_1" fi \
if ":VGRD:" rpn "sto_2" fi \
if_reg 1:2 \
rpn "rcl_1:sq:rcl_2:sq:+:sqrt:clr_1:clr_2" \
set_var WIND \
grib_out out.grb
line 1: only process the U and V at 500 mb
line 2: store U 500mb in register 1
line 3: store V 500mb in register 2
line 4: if (register 1 and register 2 have values then
line 5: calculate the wind speed: sqrt(reg_1**2 + reg_2**2)
line 6: set variable time to WIND (wind speed)
line 7: write out the WIND data to a grib file
Note: this is a very simple script and that doesn't check the matching
date code, grid type, etc.
Note: there are options to calculate wind speed and wind direction
Suppose someone made a mistake and the latent heat flux (LHTFL) had the wrong sign. RPN to the rescue.
$ wgrib2 a.grb match ":LHTFL:" rpn "1:*" grib_out new_lhtfl.grb
You could fix the entire file by
$ wgrib2 a.grb if ":LHTFL:" rpn "1:*" fi grib_out new.grb
It would be faster if you only compressed the LHTFL fields. (grib uses the
original compressed data and grib_out uses the "data" register.)
$ wgrib2 a.grb set_grib_type jpeg \
not_if ":LHTFL:" grib new.grb if ":LHTFL:" grib_out new.grb
If both the latent and sensible heat fluxes needed a sign reversal, you could do,
$ wgrib2 a.grb if ":(LHTFLSHTFL):" rpn "1:*" fi grib_out new.grb
If you want to set certain values to undefined, you define a mask and then
apply the mask. In this example, values below 20 are set to undefined.
$ wgrib2 a.grb rpn "dup:20:>=:mask" grib_out set_grib_type c3 new.grb
The RPN calculator is used:
dup the data is duplicated
20 20 is pushed on the stack
>= test data >= 20, top of stack is 1/0 depending on test >= 20
mask apply mask to the data
set_grib_type c3 sets the grib compression to complex3
grib_out new.grb writes a grib message using the decoded data
Don't forget to enclose the argument to rpn in quotes because the shell can do unexpect things.
Printing operators
print_corr write cosine weighted spatial correlation R(TOP1), R(TOP)
print_max write max(R(TOP)) to stdout
print_min write min(R(TOP)) to stdout
print_rms write cosine weighted RMS(R(TOP1)R(TOP))
Another example: Merging
Suppose that we have a nested model, we have a low resolution TMP2m from the
the outer model and a highresolution TMP2m from the nested model. Now we
want a field that uses the TMP2m in the nestedmodel domain and the TMP2m from
the outer model elsewhere. To do this, you need to convert both fields to a
common grid. Then you use "rpn merge". Make sure that both domains are
contained in the common grid as this is a requirement of the interpolation library.
wgrib2 OUTER.T2m new_grid_winds earth new_grid A B C A1.grb
wgrib2 NESTED.T2m new_grid_winds earth new_grid A B C A2.grb
wgrib2 A2.grb rpn sto_1 import_grib A1.GRB rpn "rcl_1:merge" \
grib_out MERGED.T2m
Another example: dewpoint, totaltotal index
An example of calculating the dewpoint and
totaltotal index is more involved. Using an online infix to postfix (reverse polish)
calculator is helpful.
Yet another example: global precipitation
The model has the precipitation in the variable PRATE
which has units of mm/sec (assuming 1 gm H2O = 1cc). Suppose
I wanted the globally averaged precip in the nonmetric unit
of mm/day. It's one command away
wgrib2 gfsfile match PRATE s rpn "86400:*" stats
Comments
Warning: Reverse Polish notation can cause headaches if you try something
too complicated. An infix > postfix calculator is the suggested
remedy.
The rpn option is a piece of easy to
understand and modify code (RPN.c). If you want to add
a specialized function (ex. wind chill calculation), you many
consider adding it to the RPN calculator. Why an RPN calculator?
Well, wgrib2 is heavily influenced by the stack language Forth.
It's only natural that the calculator would be based on reverse
Polish notation.
See also:
if,
if_reg,
fi,
grib_out,
rpn_rcl,
rpn_sto,
