How to create a weighted heatmap

Hi,

I would like to plot a weighted density plot based on occurence of events at a certain location.

I have a matrix with ( x-coordinate | y-coordinate | value), and I would like to use the value as a scaling factor.

Is it possible or would it be better to directly feed the result from another software then use grdimage ?

Thanks !
Guillaume

See the blockmean|median modules. They can compute bin statistics and create grids (GMT6 only).

I ended up using Matlab (accumarray). Now I’m writing the code with grdimage, will see soon :slight_smile:

Thanks !

Alright … sooo…

I have a matrix stored in a NetCDF file.

gmt grdimage -R0/10/0/10 -JX10c/10c $file?matrix -Ccolor.cpt -O -K -nn+bn >> $psFile

It works.

(1) However I’d like to switch the x-axis and the y-axis. The -p90 option makes thing a little messy, is there another way ?

(2) I don’t want the grid to be interpolated (categorical), thus the -nn+bn options, but it doesn’t work. Any idea ?

Thanks,
Guillaume

You mean plotting in portrait instead of landscape? Add -P to the grdimage command.

No no… I mean literally switch the axis.

My matrix is f([x,y]), and I want to plot it f([y,x]).

?

I think that you should use -:

If your grid is actually transposed, you will need to transpose it again. See grdedit -Et.

Thanks,

Now about the interpolation ?

This is what I got :

And this is what I try to reproduce (ignore colors and aspect ratio):
from_matlab

I also considered this example, but I’m not sure it is actually easier…

Thanks

Sorry, I dont know enough about your grid to answer. What are the grid increments? You are doing A Cartesian plot so there is no interpolation here. That smoothing must be in the data, unless we are looking at a screendump of PDF plot on macOS which likes to do that to trick us.

I indeed am using macOS right now …

The matrices are like this :
(matrix 1)
n 10degrees wide bins in X
m 1 month wide bins in Y

example:

0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ;
0 ,0 ,0 ,0 ,1.7 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ;
0 ,0 ,1.7 ,0 ,5.3 ,8.9 ,1.7 ,0 ,5.3 ,5.3 ,3.5 ,0 ,0 ;
1.7 ,0 ,3.5 ,7.1 ,10 ,3.5 ,5.3 ,1.7 ,3.5 ,3.5 ,8.9 ,8.9 ,0 ;
0 ,0 ,0 ,1.7 ,0 ,0 ,1.7 ,3.5 ,0 ,0 ,0 ,0 ,0 ;
0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ;

The full command I’m running is :

gmt grdimage -R$lon1/$lon2/0.5/12.5 -JX10c/10c $hist2D?matrix1 -Ccnt.cpt -K -O >> $psFile

Thanks again…

I fail to see the relation between the two figures. The one at top has clearly a much higher resolution. But if I take your numbers and do a quick test on Julia, I get this (not so different from your bottom fig)

julia> m = [0 0 0 0 0 0 0 0 0 0 0 0 0 ;
       0 0 0 0 1.7 0 0 0 0 0 0 0 0 ;
       0 0 1.7 0 5.3 8.9 1.7 0 5.3 5.3 3.5 0 0 ;
       1.7 0 3.5 7.1 10 3.5 5.3 1.7 3.5 3.5 8.9 8.9 0 ;
       0 0 0 1.7 0 0 1.7 3.5 0 0 0 0 0 ;
       0 0 0 0 0 0 0 0 0 0 0 0 0];

julia> imshow(m, fmt=:png, Vd=1)
        grdimage  -JX12c/0 -Baf -BWSen -R0.5/13.5/0.5/6.5 -C -P -K > C:\TMP\GMTjl_tmp.ps

Hi,
Ok so … to lower the resolution I need to use -R0.5/6.5/… instead of -R$lon1/$lon2/… and re-generate the labels manually ? (-Bpxcxlabels.txt)

The top on is interpolated by Preview and does not reflect data resolution. Convert the ps to a PNG with psconvert and you should get something close to what you expect.

What kind of sorcery is this ?
It works. So basically, if I want transparency (related to my earlier post with Ghostscript) and not interpolation, I’d rather save in PNG format than PDF. Is that correct ?
The thing is I wrote the code to have a single PS file (with multiple pages). Is there a trick to keep 1 ps and convert to multiple png ?

Apple has decided taht low-resolution images in PDF should be displayed as interpolated as that is what you want (…). It is just a Preview issue. Open in Acrobat Reader or other PDF viewer and it is fine. Transparency is always processed via PDF then to the raster and is unaffected by Preview. And yes, rasters are single page items so best to keep illustrations in single documents.

Ok, so transparency necessarily needs ghostscript 9.50 or 9.53 (hopefully), and anything other than Preview will display fine grids. GOT IT !

(just tried and confirmed) Many thanks for these insights.