# Examples showing the usage of grdmath **URL:** \ **Category:** Showcase\ **Created:** [June 1, 2025, 1:55am UTC](https://forum.generic-mapping-tools.org/t/examples-showing-the-usage-of-grdmath/5977 "2025-06-01T01:55:33Z")\ **Posts on this page:** 10\ **Page:** 1
**Author:** ![timhume](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/timhume/32/296_2.png) [@timhume](https://forum.generic-mapping-tools.org/u/timhume)\ **Post date:** [June 1, 2025, 1:55am UTC](https://forum.generic-mapping-tools.org/t/examples-showing-the-usage-of-grdmath/5977/1 "2025-06-01T01:55:33Z")
After a [recent question](http://forum.generic-mapping-tools.org/t/adding-columns-using-gmt-math/5937/9) on the use of [gmt math](https://docs.generic-mapping-tools.org/dev/gmtmath.html) I decided to investigate the [grdmath](https://docs.generic-mapping-tools.org/dev/grdmath.html) module which is analogous to [math](https://docs.generic-mapping-tools.org/dev/gmtmath.html), but for operating on grids instead. And what I found was very pleasing, so here is a simple example to demonstrate the sorts of things it can do. The example presented is computing relative humidity of the air from grids of temperature and [dew-point temperature](https://en.wikipedia.org/wiki/Dew_point). Some meteorological models only provide the temperature and dew-point temperature, so if one wants relative humidity, it is necessary to compute it from these two related variables. The way to compute relative humidity is to use this formula: rh = 100\frac{e\_{sat}(T\_{d})}{e\_{sat}(T)}, where rh is the relative humidity in %, T\_{d} is the dew-point temperature in ℃, T is the temperature in ℃ and e\_{sat}(T) is the [saturation vapour pressure of water](https://en.wikipedia.org/wiki/Vapour_pressure_of_water) at temperature T. There are several equations for computing saturation vapour pressure. [Tetens’ equation](https://en.wikipedia.org/wiki/Tetens_equation) is a simple one which has an accuracy better than ±0.1% at temperatures normally encountered near the Earth’s surface (there are complications with sub-zero temperatures which we will not discuss here): e\_{sat}(T) = 6.0178\exp\left(\frac{17.27T}{T+237.3}\right), where e\_{sat}(T) has units of hPa, and T has units of ℃. Both [math](https://docs.generic-mapping-tools.org/dev/gmtmath.html) and [grdmath](https://docs.generic-mapping-tools.org/dev/grdmath.html) support macros to simplify writing equations that use common formulae. So we start by defining some useful macros: ```auto #!/usr/bin/env bash # GMT grdmath macros. cat << END > grdmath.macros K2C = 273.15 SUB : Usage: tk K2C to return temperature in Celsius given temperature in Kelvin. ESAT = STO@TC 17.27 MUL 237.3 RCL@TC ADD DIV EXP 6.1078 MUL : Usage tk ESAT to return saturation vapour pressure in hPa given temperature in Celsius using Tetens equation. END ``` The first macro converts temperatures in Kelvin (the units usually used to disseminate temperature data) to Celsius. The second macro converts temperature in Celsius to saturation vapour pressure in hPa using Tetens’ equation. I have gridded temperature and dew-point temperature data in NetCDF format from the European Centre for Medium-range Weather Forecasts. The grid extends from 180°W to 0° to 180°E, and from 90°S to 90°N, the spacing is 0,25° in both the longitude and latitude dimensions. The data are in a file called `fld.nc`, and the temperature and dew-point temperature variables are `t2m` and `d2m` respectively. Here is the [gmt grdmath](https://docs.generic-mapping-tools.org/dev/grdmath.html) command to compute the relative humidity (using the previously defined macros): ```auto gmt grdmath fld.nc?d2m K2C ESAT fld.nc?t2m K2C ESAT DIV 100 MUL = rh_world.nc ``` And plotting it as well: ```auto gmt begin rh_world PNG gmt grdmath fld.nc?d2m K2C ESAT fld.nc?t2m K2C ESAT DIV 100 MUL = rh_world.nc gmt makecpt -Chaxby -T0/100/5 gmt grdimage rh_world.nc -JW15 -Rd -BWsen -Bpxf10a60g60 -Bpyf10a30g30 gmt coast -Dl -Wthinnest,black gmt colorbar -Bpxf5a10+l"Relative humidity (%)" -DJCB+h+w12/0.3+jCT+o0/0.5 gmt end ``` ![rh_world](https://canada1.discourse-cdn.com/flex047/uploads/gmt/original/2X/3/38a25e48d1d62ddb45899977ef3e9ad3588760a1.png) So far, so good. But here is the thing which really made me happy. The region I’m interested in crosses the International Dateline (180°E or thereabouts - it has some doglegs in it). Plotting data for an area which wraps from one edge of our grid to the other edge usually presents all sorts of hassles which need coding gymnastics to solve. But [grdmath](https://docs.generic-mapping-tools.org/dev/grdmath.html) handles it seamlessly without any pain. Here is an example for the Southwest Pacific which crosses the edge of our grid: ```auto gmt begin rh_swp PNG gmt grdmath -R160/200/-30/-5 fld.nc?d2m K2C ESAT fld.nc?t2m K2C ESAT DIV 100 MUL = rh_swp.nc gmt makecpt -Chaxby -T0/100/5 gmt grdimage rh_swp.nc -JM15 -BWseN -Bpxf5a10g10 -Bpyf5a10g10 gmt coast -Dh -Wthin,black gmt colorbar -Bpxf5a10+l"Relative humidity (%)" -DJCB+h+w12/0.3+jCT+o0/0.5 gmt end ``` ![rh_swp](https://canada1.discourse-cdn.com/flex047/uploads/gmt/original/2X/a/a4fb7931f6794f0a456d95b92aab409ce90d2be4.png) There are more modern methods for doing this sort of computation - such as the [Julia interface to GMT](https://www.generic-mapping-tools.org/GMTjl_doc/). But sometimes the old ways are pretty good too. ---
**Author:** ![timhume](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/timhume/32/296_2.png) [@timhume](https://forum.generic-mapping-tools.org/u/timhume)\ **Post date:** [June 1, 2025, 6:17am UTC](https://forum.generic-mapping-tools.org/t/examples-showing-the-usage-of-grdmath/5977/2 "2025-06-01T06:17:20Z")
As alluded to in the previous post, there are other considerations when the temperature is below 0℃. For temperatures greater than or equal to 0℃, [Tetens’ equation](https://en.wikipedia.org/wiki/Tetens_equation) above is used: e\_{sat}(T) = 6.0178\exp\left(\frac{17.27T}{T+237.3}\right). For temperatures below 0℃, [a modified equation](https://en.wikipedia.org/wiki/Tetens_equation) is used: e\_{sat}(T) = 6.0178\exp\left(\frac{21.875T}{T+265.5}\right). Here is a macro which computes e\_{sat}(T) for temperatures above and below 0℃. It’s a bit messy, but it works: ```auto ESAT = STO@TC 0 LT 21.875 17.27 IFELSE STO@A POP RCL@TC 0 LT 265.5 237.3 IFELSE STO@B POP RCL@TC RCL@A MUL RCL@TC RCL@B ADD DIV EXP 6.1078 MUL : Usage tk ESAT to return saturation vapour pressure in hPa given temperature in Celsius using Tetens equation. ``` What this does is save a variable called A for the coefficient on the top of the fraction (17.27 or 21.875) and another called B for the coefficient on the bottom of the fraction (237.3 or 265.5). The value of A and B are determined with the RPN fragments: ```auto STO@TC 0 LT 21.875 17.27 IFELSE STO@A and RCL@TC 0 LT 265.5 237.3 IFELSE STO@B ``` These A and B coefficients are then used later in the equation. In practice there’s not much visible difference between the relative humidity computed using different e\_{sat}(T) formulae for above and below 0℃. One thing I discovered is that `grdmath` macros cannot refer to other macros. If this had been possible, then I could have made the macro for e\_{sat}(T) much clearer than it is. ---
**Author:** ![Joaquim](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/joaquim/32/16_2.png) [@Joaquim](https://forum.generic-mapping-tools.org/u/Joaquim)\ **Post date:** [June 1, 2025, 11:16am UTC](https://forum.generic-mapping-tools.org/t/examples-showing-the-usage-of-grdmath/5977/3 "2025-06-01T11:16:35Z")
Tim, you have to open a Weather Channel here in the Forum 😄 ---
**Author:** ![Joaquim](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/joaquim/32/16_2.png) [@Joaquim](https://forum.generic-mapping-tools.org/u/Joaquim)\ **Post date:** [June 1, 2025, 11:58am UTC](https://forum.generic-mapping-tools.org/t/examples-showing-the-usage-of-grdmath/5977/4 "2025-06-01T11:58:51Z")
And can’t resist to hijack a bit your post. With that other option, and only for the simpler case (other requires a lengthier fun) ```auto julia> ecmwf(:forecast, vars=["2t", "2d"]) [ Info: Downloading 2t and saving to 2t_step0_ifs_oper_fc_20250601_00.grd [ Info: Downloading 2d and saving to 2d_step0_ifs_oper_fc_20250601_00.grd julia> T = gmtread("2t_step0_ifs_oper_fc_20250601_00.grd"); julia> Td = gmtread("2d_step0_ifs_oper_fc_20250601_00.grd"); julia> esat(x) = 6.0178*exp.(17.27*x ./ (x .+ 237.3)) # Define a function online esat (generic function with 1 method) julia> Grh = 100*esat(Td) / esat(T); viz(Grh, proj=:guess, colorbar=true) ``` ---
**Author:** ![Esteban82](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/esteban82/32/54_2.png) [@Esteban82](https://forum.generic-mapping-tools.org/u/Esteban82)\ **Post date:** [June 1, 2025, 1:22pm UTC](https://forum.generic-mapping-tools.org/t/examples-showing-the-usage-of-grdmath/5977/5 "2025-06-01T13:22:12Z")
Nices maps @timhume!! I think that this example can be included here: [GMT Tutorials and Examples — GMT Examples documentation](https://www.generic-mapping-tools.org/gmt-examples/index.html) ---
**Author:** ![timhume](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/timhume/32/296_2.png) [@timhume](https://forum.generic-mapping-tools.org/u/timhume)\ **Post date:** [June 1, 2025, 4:47pm UTC](https://forum.generic-mapping-tools.org/t/examples-showing-the-usage-of-grdmath/5977/6 "2025-06-01T16:47:03Z")
I was hoping you’d reply with something like this 😀 ---
**Author:** ![Joaquim](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/joaquim/32/16_2.png) [@Joaquim](https://forum.generic-mapping-tools.org/u/Joaquim)\ **Post date:** [June 1, 2025, 5:08pm UTC](https://forum.generic-mapping-tools.org/t/examples-showing-the-usage-of-grdmath/5977/7 "2025-06-01T17:08:34Z")
Yep, I felt the smell of the provocation 😄 ---
**Author:** ![timhume](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/timhume/32/296_2.png) [@timhume](https://forum.generic-mapping-tools.org/u/timhume)\ **Post date:** [June 1, 2025, 9:13pm UTC](https://forum.generic-mapping-tools.org/t/examples-showing-the-usage-of-grdmath/5977/8 "2025-06-01T21:13:40Z")
That’s the beauty of GMT. Not only can you use shell scripts, but nowadays you can use the Julia or Python interfaces too. By the way, the just activated math plugin is great! ---
**Author:** ![PlanetGus](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/planetgus/32/303_2.png) [@PlanetGus](https://forum.generic-mapping-tools.org/u/PlanetGus)\ **Post date:** [June 2, 2025, 8:13am UTC](https://forum.generic-mapping-tools.org/t/examples-showing-the-usage-of-grdmath/5977/9 "2025-06-02T08:13:23Z")
Please, don’t fall for the dark side. CLI is the way ! ---
**Author:** ![mkononets](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/mkononets/32/4029_2.png) [@mkononets](https://forum.generic-mapping-tools.org/u/mkononets)\ **Post date:** [August 29, 2025, 4:00pm UTC](https://forum.generic-mapping-tools.org/t/examples-showing-the-usage-of-grdmath/5977/10 "2025-08-29T16:00:07Z")
There’s also `gdal_calc.py`: [gdal\_calc — GDAL documentation](https://gdal.org/en/stable/programs/gdal_calc.html) “Command line raster calculator with numpy syntax”. There are few examples on the docs web page, and I tried it but on small data sets in my texture shading scripts. Not sure how it works with a global grid that wraps around the globe. It uses “normal” syntax rather than the RPN in `grd math`.