Clustering XY locations based on a certain grid

Dear GMT developers and users,

I have XY locations corresponding to seismicity of earthquakes. I want to compute the total amount of seismic moment and earthquake rate by firstly clustering/grouping the events into a certain grid, say every 1deg, as shown by the figure attached. The outputs I expect are a number of grids, with each grid listing the number of earthquakes that fall within the grid, including their XY locations and magnitudes. Is it possible to do that in GMT?

Regards,
Rino Salman

Dear Rino

I have used the programs blockmean, blockmedian and blockmode in the past to perform similar tasks.

With blockmean you can use the argument -S to obtain the number of events in cell and also the sum of the z values, for example the seismic moment in order to obtain the cumulative seismic moment (the plot of the left in your figure). You can output grids directly and then use grdmath to perform operations on these.

If you want to obtain the events in the cells as ASCII tables maybe you have to use gmtselect in a loop selecting the events by region (-R). Something like:

for x in {xmin..xmax..dx} do
   for y in {ymin..ymax..dy} do
      gmt select input.file -Rx/x+dx/y/y+dy > output_x+dx.file
   done
done

Best regards

Dear Jose-Alvarez,

Thank you for the clue, I really appreciate it.

Best regards
Rino