Interpolation of scattered earthquake data

2 views (last 30 days)
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.
  2 Comments
Akira Agata
Akira Agata on 20 Nov 2019
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?)
scatter3.png
Michael Olaniran
Michael Olaniran on 20 Nov 2019
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

Sign in to comment.

Accepted Answer

Star Strider
Star Strider on 20 Nov 2019
Edited: Star Strider on 20 Nov 2019
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')
Interpolation of scattered earthquake data - 2019 11 19.png
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.
  4 Comments
Michael Olaniran
Michael Olaniran on 23 Nov 2019
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
Star Strider on 23 Nov 2019
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.

Sign in to comment.

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!