How to create a heatmap

Here I propose to create a heatmap (or frequency of events in a given area).

In my case, I have GPS coordinates that I want to map them on a pre-existing grid.

  • Area of interest : 235E : 250E / 30N : 45N (these values match a grid points)
  • Grid-x resolution : 0.625 deg
  • Grid-y resolution : 0.5 deg

Thanks to xyz2grd and the -An flag, I can directly “count” the number of events in each grid … convenient !
I also replaced the NaN values by "0"s (-di) and asked the function to match the coordinates exactly to my grid (+e flag)



# For the colorbar later
clr_txt="-Bxa+lCount -By"
gmt xyz2grd gps.txt -I0.625/0.5+e $region $projection -An -di0

Next, I got GPS coordinates of the borders of my area of interest (here the Great Basin, Nevada) that I want to clip my region on … in format 0:360 instead of -180:180. The -C0 flag after the file allows me to do the operation on longitude only and keep both longitude and latitude in the output. Had it be placed before, only longitude would have remained.

gmt math basin.txt -C0 360 ADD = gbas.txt

Now the plot :
0. Create a colormap

  1. Clip the basin
  2. Plot the heatmap
  3. Plot the basin
    • I shifted the gridlines to match the heatmap (-Bafg+<phase>)
  4. add a colorbar
gmt begin test01 png
	gmt makecpt -Clapaz -T0/150/10

	gmt clip $region $projection gbas.txt
	gmt grdimage -C -B+t"# events per area"
	gmt clip -C
	gmt plot -Wthinner,black gbas.txt -L -Bxa5f2g0.625-0.3125 -Bya5f2g0.5-0.25 --MAP_GRID_PEN_PRIMARY=faint,grey
	gmt colorbar ${clr_sty} ${clr_txt} -Np
gmt end show

Tadaa !


Very nice! You may want to check out gmtbinstats for applications like this as well.

1 Like

Damn, the name of it looks scary !

Beautiful map; thanks for sharing PlanetGus!

Thanks @Andreas :slight_smile: [padding too]


I have a doubt. I think that if you use degrees for the Increment (-I) then you will have smaller cells (in km2) as you approach to the poles. Is this correct? Does this happen if km/m are use instead?

Yes, that is correct. if you have an area with large latitudinal range you should project your data using an equal-area projection and then do the binning and plotting of the projected data since those equal-sized bins will represent the same physical area. Alternatively, you can work in degrees as above but you can normalize for the variable area and present values like number-of-something/km^2. See gmt grdmath AREA function for creating such weights.


I think the last example using count doesn’t exist ?

Yes, for most of the man page examples there are a few that uses actual files on the server (@) and then a bunch of imaginary examples. We hope to eliminate all imaginary examples but it is a big job.

What an idea to create such a versatile software !