Gmtbinstats: cells containing points are painted NaN in grdimage

Simple case: count the number of earthquakes in each 1° bin.

  • the grid produced is correct (checked in qgis); the number of points in each cell is correct (got confused at first since default is gridline, not pixel node registration)
  • when plotting this grid with grdimage, many of cells containing quakes (black dots) are painted gray, i.e. NaN. Why is this? I tried playing with makecpt (hence the funky colorbar), but to no avail.
  • Script and data included.
gmt begin forumpost png

gmt binstats $(gmt info eq-2000-2021.txt -I1) -I1 -Cn -Tr eq-2000-2021.txt -rp
gmt makecpt -T-10/100/10 -Cjet
gmt grdimage -C -B0 -R-5/53/90/81r -JU+34/20c
gmt coast -W -Dc
gmt plot  eq-2000-2021.txt -Sc0.05c -Gblack -l"earthquake"
gmt colorbar

gmt end show

eq-2000-2021.txt (170.7 KB)

Changing -I step solves the problem. (Try -I33k).

If I were to guess, I’d say that -I1(d) near the pole behaves “as if” it was a rectangular 1º cell (~111x111km). But it’s just a guess.

Thanks @PlanetGus, but that still does not explain why the cells with values are painted gray (NaN) by grdimage?

Except if I’m right about the cell size. Once a point is counted in a cell, it is not in other overlapping ones.

In that case, -S would be better than -T (again, this is just my guess)

I don’t get that NaNification effect

julia> G = binstats("eq-2000-2021.txt", R="-24/51/50/88", I="33k", C=:n, T=:r, r=:p);
julia> imshow("", proj=:guess, coast=true, colorbar=true)

You’re right; changing from -I1 to -I33k improves the matter, but I’m still getting gray cells that should be colored. See att. figure - I’ve highlighted some of the area that have gray cells that should have colored cells.

And why the heck does -I33k improve the matter? I want -I1!

If grdimage does any reproduction to turn your grid into a clipped rectangular image then NaNs will bleed. See -nn to avoid that.

Ah, yes. I have -n+a by default.

Still getting lots of NaNs with -nn, -n+a and -nn+a (where I should not have NaNs).

Teaching moment:

  1. Replace your grdimage command with gmt grdview -C -B0 -R-5/53/90/81r -JU+34/20c -T and realise what 1x2 degree cells look like to the north, Imagine grdimage creating an equidistant image to fill that. Even try -B0g1 to see the grid.
  2. Decide you need finer resampling of the image by adding -E600 to your current grdimage command and note how it changes,

Thanks Paul for teaching moment; both of these give reasonable results.
I prefer the grdview though, as grdimage -E600 looks smeared.
I know that 1x1 degree cells will be very thin towards the pole, but that’s OK.


  • There is not a corresponding -T option for grdimage? Why not? The only difference, in appearance at least, between the grdimage and grdview command, is that the grdview versions looks like what I expect (and want) grdimage to do.

Not sure, perhaps we decided at one time to not overload grdimage with too many options. But, your near-pole cases do suggest that perhaps the standard results (your first plots) are not a good result for users, and the remedy may require users to just know too much. Pinging our southern guru @Joaquim on what he thinks about

  1. Add -T as in grdview to grdimage [pretty easy to do: just make a function out of what grdview does and call them from those two programs]. Since grdview was written to do 3-D perspective it may not easily come to mind to use for your 2-D case.
  2. Detect cases where creating a nx by ny rectangular image either needs a boosting (e.g., a -E under the hood in a way), some warning educates the users, or we decide plotting via -T is the best default?

I don’t have much experience with grdview, so as always; pardon my ignorance.

I use grdimage 99% of the time, but my impression has always been that grdview and grdimage does have a great overlap of use cases - sometimes some rarely used feature forces me to use grdview. In that regard I have thought about grdview as being a superset of grdimage. And, in the mentality of Occam’s razor, you choose the simplest tool for the job.

I still fail to see what are those cells that are wrongly painted as NaNs. Bellow is a grdimage of -I1 and the same grid displayed in Mirone. They look the same to me.

And this now is with the new -T in grdimage. Almost equal to before, except that we now see faint outline of cells.

julia> imshow("", T=true, coast=true, colorbar=true, Vd=1)
        ["grdimage  -JX14c/0 -R-24/51/50/88 -Baf -BWSen -n+a -T -C -P -K > C:\\TEMP\\", "pscoast  -J -A200/0/2 -W0.5 -Da -K -O >> C:\\TEMP\\", "psscale  -J -Baf -DJMR -C -K -O >> C:\\TEMP\\"]

Well, just a linear projection. Then of course it look’s right. Try Andreas’s projection and then the problem shows

Ah, in my projected version, I was not plotting the points (and not even paying attention to them) to compare with colors.