earthquake density map/grid/info

Dear all,

I have been trying to use GMT to reproduce a processing step
from a 1986 report on global seismic magnitude measurements.
The author, wishing to overcome the problem that earthquakes
are not evenly distributed, divided the surface of the Earth
into 400 tessera of equal area and extracted the largest three
earthquakes from each for his analysis.  I haven't been able
to find his original program, so I thought I'd use GMT
(v. 4.5.12 on CentOS5 Linux).

I don't need my tessera to be identical to his because I'm
using a different dataset, but I do need them to be of
equal area.  My input file consists of a list of earthquakes
for which I have (ID-number, latitude, longitude, magnitude).

I adopted the following sequence, first of all to find the
number of earthquakes per grid cell (many cells are empty!).

1. mapproject (on to a Lambert Azimuthal polar projection);
2. xyz2grd with -An option to count number in each cell;
3. grdproject to return the co-ordinates to lat/long;
4. grd2xyz to get ASCII output.

These four steps are carried out separately for northern
and southern hemispheres.  The purpose of the initial
mapproject is to convert the lats/longs into Cartesian
co-ordinates on a plane in which areas on the surface of
the Earth translate into equal areas on the
plane regardless of their latitude/longitude.  This is
necessary since xyz2grd assumes a rectangular grid.

There are two problems:

1. I can use -An in xyz2grd to get the number of earthquakes
per cell, or -Au to find the largest-magnitude earthquake
in each cell, but I can't, for instance, find the largest
three earthquakes in each cell, or find the ID number of
the largest earthquake - unless I have missed something
in xyz2grd that would allow me to store an array as the
"z" value;

2. Plotting is, so far, a travesty, because grdimage
plots each grid cell as a square or rectangle centred on
the grid node, rather than fitting the grid into the
circular Lambert projection.   Getting a nice plot of
my tessera on a global map (colour-coded to represent
number-of-earthquakes) is of secondary importance to
extracting the right earthquakes to feed to the program
featured in the 1986 report, but I know it will become of
primary importance once I get the program to produce some
publishable results.

Advice welcome.

Regards,
Sheila Peacock,
Blacknest.

To unsubscribe, send the message "signoff gmt-help" to [email protected]