GMT Part 3: Working with NetCDF data by creating a synoptic weather map
For my thesis I had to examine synoptic scale circulation over
South-Africa, I wanted to recreate a classic synoptic map that
represented mean sea level pressure and also geopotentail at 500hPa.
I knew the South-African Weather Service
(SAWS) has daily
make publicly available, but the maps are not always the highest
quality as they are hand drawn - a tedious process I would imagine.
A SAWS synoptic map.
In addition to this I was also attempting to better understand GMT’s
handling of NetCDF files. I used
data in my thesis for WRF and also some objective clustering, the ERA
dataset was thus the obvious choice. Using GMT to recreate these
synoptic maps took me a couple of days, mostly the time consumption
was related to handling the NetCDF data. GMT uses Reverse Polish
Notation (RPN) for math functions which I needed to figure out as
you’ll see. It is also useful to have CDO, ncview and ncdump installed
for data inspection and manipulating. A last point is that I compiled
it all in a bash script, this is useful for various reasons, so just
learn to do it.
First off a synoptic map is pretty simple, it consists of pressure
levels and if you really want you can add some other variables such as
CAPE, temperature, precipitate water, and so forth. Synoptic mainly
refers to a scale in this context if you wondered, a 1975
by Isidoro Orlanski explains this concept of scale related to weather
and climate clearly. You’ll notice is that the projection used in the
SAWS map is not the Mercator projection I’ve used previously, but
instead a equal area
namely Lambert Azimuthal
Coincidently someone already explained the projection using
South-Africa as an
A quick drawing of the basemap and we can see th projection is
This does not mean a lot without any land referance, the -K option
tells us we want to add something later to the map. Using pscoast and
-Df and -N1 options we can add coastlines and national borders to this
That was pretty straight forward, we can see that the map is more or
less similar to the synoptic map we want to make. Now we need to add
the pressure lines to the map, this gets a little more tricky. As
mentioned, I used ERA-Interim
data in my thesis, the data is maintained by the
ECMWF and publicly available, you do however
need to sign up to use the data. You can download the
data directly from the website, it is however easier to use the ECMWF
Do this to use the script I wrote below, otherwise it will
not work and you’ll think I’m stupid. So, once you’ve got that set up
the rest will make a little more sense.
I initially tried to just run gmt grdimage without any data conversion, this
did not work. Inspecting the ordinal data with,
Shows us the msl variable is is present,
I think the data needs to be in a 2-dimensional grid, explained
in more detail. To get a GMT compliant file we need to do,
Inspecting the data shows us that z has replaced msl
We can also see that 95334 is not a pressure value, we have to divide
it by 100. It is possible to use CDO or NCO tools, but the powerful
GMT has math functions in the form of Reverse Polish Notation
convert new mslp value by a 100 we can do
Because the data covers the whole word we can make a quick plot
without to much trouble using grdimage to see if everything seems
This gives us a image of the mslp scaled and covering the whole world.
We do however want a contoured mslp and also clip it to South-Africa.
To do this we need to call the psbasemap and pscoast functions we used
earlier and then use grdcontour to contour the data.
We get a nice map with the MSLP contoured over South-Africa, much like
the original synoptic map.
We can see that the ERA created map captures similar features as
the SAWS synoptic map; note the South-Indian Ocean High Pressure and
the South-Atlantic High Pressure and the low pressure from the
Easterly Wave over the Country. To add the L and H indicators I would
use inkscape at this time, a script would be possible but that is
beyond me right now. The weather stations indicators would also theoretically be
possible and I’ll try to figure that out in the future.
The below map is the final map I created, you can see The script below
does all the above steps with the addition of the upper 500hPa
geopotential. I also added some color. In the end we end up with a new
and improved (in my subjective opinion) synoptic map as seen below.