How to create a heatmap

PlanetGus · March 20, 2021, 5:59pm

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)

lonW=235
lonE=250
latS=30
latN=45
region="-R${lonW}/${lonE}/${latS}/${latN}"

projection="-JX10c"

# For the colorbar later
clr_sty="-DjCL+ef+w8c/0.3c"
clr_txt="-Bxa+lCount -By"

gmt xyz2grd gps.txt -Ggridded.nc -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

Clip the basin
Plot the heatmap
Plot the basin
- I shifted the gridlines to match the heatmap (-Bafg+<phase>)
add a colorbar

gmt begin test01 png
	gmt makecpt -Clapaz -T0/150/10

	gmt clip $region $projection gbas.txt
	gmt grdimage  gridded.nc -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 !

pwessel · March 20, 2021, 6:10pm

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

PlanetGus · March 20, 2021, 6:12pm

Damn, the name of it looks scary !

Andreas · March 21, 2021, 9:31am

Beautiful map; thanks for sharing PlanetGus!

PlanetGus · March 21, 2021, 3:30pm

Thanks @Andreas [padding too]

Esteban82 · March 25, 2021, 3:52am

Nice!!

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?

pwessel · March 25, 2021, 4:15am

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.

PlanetGus · April 28, 2021, 4:02pm

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

pwessel · April 28, 2021, 9:03pm

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.

PlanetGus · April 28, 2021, 10:10pm

What an idea to create such a versatile software !

PlanetGus · March 8, 2023, 10:23am

UPDATE

And to help @annafiglioli (see her thread for data)

#! /bin/bash

## PARAMETERS
data_input="prova_GMT.txt"
data_gridded="prova_GMT_gridded.nc"

lon_start=11.50894
lon_stop=15.77485
lat_start=35.27988
lat_stop=39.50000
depth_min=0
depth_max=400

grid_step="10k"

map_width=5
profile_height=3

## STYLE
# colorbar later
clr_sty="-DJCR+ef+w4c/0.2c"
clr_txt="-Bxa+lCount -By"

## GMT SCRIPT
gmt begin prova png  

    gmt makecpt -Coleron -T0/100/5 -H -D > heat.cpt 

    gmt gmtbinstats -R${lon_start}/${lon_stop}/${lat_start}/${lat_stop} \
                    ${data_input} -G${data_gridded} \
                    -Cn -I${grid_step}+e -Tr -E0

    gmt coast -JM${map_width}c -G -Bafg -B+glightblue+t"Sicilia"
        gmt grdimage ${data_gridded} -Cheat.cpt
		gmt coast -Wthin -Di
    gmt coast -Q
	gmt colorbar ${clr_sty} ${clr_txt} -Np

	gmt plot -R${lon_start}/${lon_stop}/${depth_min}/${depth_max} -JX${map_width}d/-${profile_height} \
				${data_input} -i0,2 -Sc0.1c -Glightgray -Wfaint \
				-Bxafg -Byafg+l"Depth (m)" -Y-4c

gmt end show

Ideally, the clipping with coast would use DCW (iso-3166) but, in this case, Sicilia (IT-82) is not part of the available code. Only “IT” or “italy” is.

Esteban82 · March 8, 2023, 12:00pm

But you could use dcw-collections:

tag: SCLI Sicily Island
region: 12.421/15.653/36.65/38.303

PlanetGus · March 8, 2023, 12:04pm

Ah!
We definitely need a table with all the code available.
I tried stuff like IT.SI and so on but none worked.
As long name (e.g. -Etexas) works, why -Esicily did not?

Esteban82 · March 8, 2023, 5:26pm

You should try -ESCLI or -E"Sicily Island". The longitude and latitude values are related to the geogrpahic feature (the island) and NOT the political region. In this case, the difference are the minor island like Lampedusa and others.

PlanetGus · March 8, 2023, 7:36pm

Just saying, something in the doc ‘round here would be useful. The same way we have octal codes and symbols or built-in cpts

Esteban82 · March 8, 2023, 7:54pm

Just a list? or do you want also maps?

PlanetGus · March 8, 2023, 8:08pm

Just a list so we can do a quick search in the page.

Alphabetically is fine, per continent eventually? I don’t know how deep in the tree we can go without too much work.

Option A :

a fullname or code
b fullname or code
…

Option B :

Africa
- a fullname or code
- b fullname or code
- …
Asia
…

Option C :

Continent 1
- Country i
  - region a
  - region b
  - …
- Country ii
  …

I have a preference for Option B. It’s ordered enough with less work and is less controversial too (e.g. Taïwan, Cachemire, Nepal, Donbass, … love politics )

Esteban82 · March 8, 2023, 8:30pm

The list would be about 250 records. It seems to me a bit long.

Link:

github.com

GenericMappingTools/dcw-gmt/blob/master/dcw-collections.txt

# Data from United Nations geoscheme M49 (https://en.wikipedia.org/wiki/UN_M49)
tag: UN002 UN Africa
list: DZ,EG,LY,MA,SD,TN,EH,IO,BI,KM,DJ,ER,ET,TF,KE,MG,MW,MU,YT,MZ,RE,RW,SC,SO,SS,TZ,UG,ZM,ZW,AO,CM,CF,TD,CG,CD,GQ,GA,ST,BW,SZ,LS,NA,ZA,BJ,BF,CV,CI,GM,GH,GN,GW,LR,ML,MR,NE,NG,SH,SN,SL,TG
tag: UN015 UN Northern Africa
list: DZ,EG,LY,MA,SD,TN,EH
tag: UN202 UN Sub-Saharan Africa
list: IO,BI,KM,DJ,ER,ET,TF,KE,MG,MW,MU,YT,MZ,RE,RW,SC,SO,SS,TZ,UG,ZM,ZW,AO,CM,CF,TD,CG,CD,GQ,GA,ST,BW,SZ,LS,NA,ZA,BJ,BF,CV,CI,GM,GH,GN,GW,LR,ML,MR,NE,NG,SH,SN,SL,TG
tag: UN014 UN Eastern Africa
list: IO,BI,KM,DJ,ER,ET,TF,KE,MG,MW,MU,YT,MZ,RE,RW,SC,SO,SS,TZ,UG,ZM,ZW
tag: UN017 UN Middle Africa
list: AO,CM,CF,TD,CG,CD,GQ,GA,ST
tag: UN018 UN Southern Africa
list: BW,SZ,LS,NA,ZA
tag: UN011 UN Western Africa
list: BJ,BF,CV,CI,GM,GH,GN,GW,LR,ML,MR,NE,NG,SH,SN,SL,TG
tag: UN019 UN Americas
list: AI,AG,AW,BS,BB,BQ,KY,CU,CW,DM,DO,GD,GP,HT,JM,MQ,MS,PR,BL,KN,LC,MF,VC,SX,TT,TC,VG,VI,BZ,CR,SV,GT,HN,MX,NI,PA,AR,BO,BV,BR,CL,CO,EC,FK,GF,GY,PY,PE,GS,SR,UY,VE,BM,CA,GL,PM,US
tag: UN419 UN Latin America and the Caribbean
list: AI,AG,AW,BS,BB,BQ,KY,CU,CW,DM,DO,GD,GP,HT,JM,MQ,MS,PR,BL,KN,LC,MF,VC,SX,TT,TC,VG,VI,BZ,CR,SV,GT,HN,MX,NI,PA,AR,BO,BV,BR,CL,CO,EC,FK,GF,GY,PY,PE,GS,SR,UY,VE
tag: UN029 UN Caribbean

This file has been truncated. show original

PlanetGus · March 8, 2023, 8:57pm

It actually works for me.
It’s a little ugly to read through, but quick search works.

Though, real names would be more useful.

For example, in the Lesser Antilles category, MQ refers to Martinique.
But if I don’t know the code, I’ll need to resort to Wikipedia.

I think it is important that each available code was referenced with its fullname

(Maybe we should go to github instead)

Raisha · April 14, 2024, 11:48pm

Hi, is there any clue on how to do exactly this, but with a smoothing applied instead of patchy density between the blocks?