Re: earthquake density map/grid/info

Hi, All,

It seems that if you have the coordinates of the tessera from the earlier work you can use an algorithm to search for events in all the tessera. These events can then be sorted and you have the full list available to you. 

I have an algorithm for selection referred to above. Let me know if you need it.

Frarncis


On Mar 19, 2014, at 11:29 AM, Walter (HF) Smith <[email protected]> wrote:

> Dear Sheila,
> 
> Thanks for your faith in GMT. :)
> 
> There are a lot of ways to do this, but perhaps getting the top three is the hard part, at least as GMT is set up now.
> 
> I think you will be disappointed with using "grdimage" because the "image" operator in PostScript necessarily paints areas with rectangular boundaries.
> 
> How about this:
> 1, mapproject all your data into some equal-area projection (more about that below)
> 2, something (if not block mean/median/mode then maybe an awk) to print to each column an index, k, or pair of indecies, i and j, which correspond to x/size and y/size rounded to integers. In this step x,y are the projected equal-area coordinates and the idea is to create integers that specify what "tile" each data point belongs in.
> 3, sort the output of 2 (unix sort, see -n, -r and -k options) so that the first sort key or keys are the index k or the indices i and j, and then the last sort key is the earthquake "size" (I don't think you said what magnitude estimate or moment estimate this was, but it doesn't matter, let's call it size). By doing a reversed (-r) numerical (-n) sort, you'll get the quakes in each cell listed from largest to smallest. 
> 4, the next step is a while loop that will print out from zero up to 3 (in case a tile isn't full) records from the start of each new tile. This can again be done with awk.
> 
> I'm a dinosaur so I am familiar with awk. Others on this list will be more facile with perl or other ways of doing this.
> 
> Note that the above mainly bypasses GMT after the mapproject step. Others who are better GMT users than I (yes, you know who you are!) may be able to see a way to do it with more GMT tools. GMT's functionality has grown enormously but my ability to think in GMT is not always keeping up.
> 
> Walter
> 
> On Mar 19, 2014, at 2:11 PM, Sheila Peacock <[email protected]> wrote:
> 
>> 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]
> 
> To unsubscribe, send the message "signoff gmt-help" to [email protected]


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