Interpolation of scattered earthquake data

Hello all, I am currently stuck trying to write a mat lab code that can interpolate my scattered earthquake data and produce a file that contains a continuous grid of data. I have an excel file that contains 3 columns. The first column is Longitude (Lo), the 2nd column is Latitude (Lt) while the third is earthquake magnitude (Mg). Please how can I go about this?. I am trying to do a latitude and longitude grid spacing of one degree each and export the resulting grid as either an excel or CSV file.

I look forward to your response.

[Akira Agata]

Looking at your data, in some area, there are many data points with almost the same Longitude/Lattitude but with different Magnitude. So I think it's not good to apply interpolation using all these data points. Is it possible to apply some pre-processing for your data before interpolation? (such as averaging the Magnitude for each grid square area?)

[Michael Olaniran]

Thank you Akira for your quick repsonse. So I am trying to sum up all the magnitudes that occur which a particular cell grid. Let say the area of my grid cell is 1 degree by 1 degree, I want to sum up all earthquakes that occur within that grid cell and progressively do that for all the grid cells within the boundary limit (Lattitude 10 -60, Longitude 60-140) of my data

Answer 1

[Star Strider]

There are several options. First, see if the griddata function will do what you want. (The Mapping Toolbox may have more applicable functions, however I do not have it, so I have no experience with its functions.)

EDIT — (19 Nov 2019 at 3:51)

[D,S] = xlsread('Earthquakedata.xlsx');

Lon = D(:,1);
Lat = D(:,2);
Mag = D(:,3);

LonX = [min(Lon) max(Lon)];
LatX = [min(Lat) max(Lat)];

LonV = floor(LonX(1)):ceil(LonX(2));                        % One-Degree Resolution
LatV = floor(LatX(1)):ceil(LatX(2));                        % One-Degree Resolution

[LonM,LatM] = meshgrid(LonV, LatV);
MagM = griddata(Lon, Lat, Mag, LonM, LatM);

figure
surfc(LonM, LatM, MagM)
grid on
xlabel('Longitude')
ylabel('Latitude')
zlabel('Magnitude')

The griddata call produces this Warning:

Warning: Duplicate data points have been detected and removed - corresponding values have been
averaged. 
> In griddata>useScatteredInterp (line 185)
  In griddata (line 126)

So griddata does the ‘heavy lifting’ on its own, without any need to pre-process the data.

[Michael Olaniran]

Thank you @Star strider. This helps alot . I am trying the understand the 1 degree resolution line you wrote. If I wanted to do like 0.5 resolution how can I do that?

[Star Strider]

My pleasure.

To make it 0.5° resolution, simply change the step size from 1 (implied if not specifically stated) to 0.5:

LonV = floor(LonX(1)):0.5:ceil(LonX(2));                    % One-Half-Degree Resolution
LatV = floor(LatX(1)):0.5:ceil(LatX(2));                    % One-Half-Degree Resolution

There are other ways to create these vectors, although this is likely the easiest in this situation.

[Michael Olaniran]

Thanks @Star strider . I just thought of an idea I am playing around with- inserting the old magnitude data in the new o.5 resolution grid without interpolating the magnitude. So the old magnitude will maintain their longitude and Lattitude locations while the new grid points will be Replaced by 1s.So you eventually get the old data and 1s rather than interpolated old data.Please how can I go about this?

[Star Strider]

As always, my pleasure.

I do not understand what you want to do. In any event, I would not replace the new grid points by 1s, I would use NaN since NaN values will not plot. This would leave blanks for the NaN values, which is apparently what you want. To do that, first create a NaN matrix with the appropriate dimensions, then put the known data at the correct ‘Lon’ and ‘Lat’ points. You would have to round these to integers if you want to use them as matrix subscripts, so to get fractional degree resolution, multiply all of them by 10, or more if you wanted even greater resolution. Note that according to the Warning that griddata produced, there are also duplicate coordinates, so you will have to devise a way to deal with those.

However, I am not certain that you can plot the data without interpolating them, regardless, especially on any sort of regular grid. The reason is that they are essentialy random, not gridded.

This is straightforward to demonstrate that graphically, using the ‘Lon’ and ‘Lat’ vectors from my earlier code, and data that are designed to be gridded:

figure
plot(Lon, Lat)
grid
title('Earthquake Coordinates')

[X,Y] = ndgrid(0.5:0.1:1.5);
figure
plot(X(:), Y(:))
grid
title('Gridded Coordinates')

To emphasize that, do this plot as well:

figure
stem3(Lon, Lat, Mag)
grid on

You can also use plot3, however I prefer stem3 here because it more accurately depicts the dimensions of each point.