How can I calculate distance between a point in a given radius inside a geographic coordinate grid?

Hello all,

I tried to search for similar questions, and I have got different solutions, but nothing to combine all my tasks, hence I still can't find the proper approach for my problem. Some of the approaches include:

https://www.mathworks.com/matlabcentral/answers/492031-interpolation-of-scattered-earthquake-data

https://www.mathworks.com/matlabcentral/answers/423998-circles-with-particular-radius-for-latitude-and-longitude-points-plotted-on-map

I need to create a grid which contains earthquake data, namely longitude and latitude points (N & E geographic coordinates). What I need is to divide the area into a 0.5 km interval 2-D grid in the horizontal coordinates.

For each grid point, I would need to plot a circle with radius 6 km, and count how many events fall within the circle as well as calculate the distance between the center point of the circle and each event point (earthquake event).

Attached is the file containing the earthquake data needed for the analysis.

Cheers,

H.

[darova]

This is how i understand your question

Is this a correct interpretation?

[Harith Al Kubaisy]

Yes exactly, thanks for the follow up. The grid would be uniformly horizontal however.

[darova]

You can easily identificate indices of points inside circle by

% x1,y1 - circle center
% x,y - data
% R - radius of a circle
ix = (x-x1).^2 + (y-y1).^2 <= R^2

[Harith Al Kubaisy]

I will give that a go, thanks. What about the distances between the circle center, and said data points? And is there a way to modify the code to make it draw circles and calculate the needed distannces for every 0.5 km in the horizontal coordinates?

[Harith Al Kubaisy]

suppose I have this bit of code (sorry for the extra zeros)

clc

clear all

close all

grid on

long =

[57;52.1000000000000;50;52;43.5000000000000;50;53.5000000000000;52

;

43.2000000000000;54;42;43.5000000000000;52;53.5000000000000;56;56

;

52;51.5100000000000;44;52;53.5000000000000;52;52;51.5000000000000;

48;55.4000000000000;48.7000000000000;56.7000000000000;42.6000000

000000;47.3000000000000;46.7000000000000;45.8000000000000;56.700

0000000000;43;57.4000000000000;51.2000000000000;56.5000000000000

;

43.1000000000000;53.4000000000000;43.6000000000000;45.4000000000

000;53.6000000000000;53.8000000000000;51.8000000000000;53.700000

0000000;51.7260000000000;42.8000000000000;46.4000000000000;54.80

00000000000;44.5000000000000;53.9000000000000;44.3000000000000;

52.2000000000000;56.4000000000000;50.2000000000000;48.900000000

0000;53.5000000000000;48.1000000000000;51.6000000000000;43;42.800

0000000000;53.6000000000000;56.4000000000000;53.7000000000000;5

3.4000000000000;46.3100000000000;49.7000000000000;51.70000000000

00;50.3000000000000;42.5000000000000;54.1000000000000;43.6000000

000000;56.5700000000000;48.4000000000000;52;56.9000000000000;57.4

960000000000;50.9000000000000;52.3000000000000;56.4150000000000

;

50.5000000000000;53.5000000000000;56.5400000000000;56.200000000

0000;50.8000000000000;56.2500000000000;52;51.6000000000000;53.700

0000000000;54.1000000000000;56.2760000000000;46.1000000000000;56

.

7400000000000;56.8510000000000;45;53.3000000000000;48.2000000000

000;51.3000000000000;51.7000000000000;53.7000000000000];

l$t =

[16;15.7000000000000;13.5000000000000;14.5000000000000;11.80000000

00000;13.5000000000000;14.5000000000000;14;11.7000000000000;14.500

0000000000;11.5000000000000;11.5000000000000;13.5000000000000;15;1

3.5000000000000;14.5000000000000;13.5000000000000;13.50000000000

00;10.9000000000000;14;14.5000000000000;14;14;13.5000000000000;13;1

3.7000000000000;12.1000000000000;13.9000000000000;11.600000000000

0;12;12.1000000000000;12.2000000000000;13.9000000000000;11.5000000

000000;13.1000000000000;13.7000000000000;14.5000000000000;11.5000

000000000;14.5000000000000;12.1000000000000;11.9000000000000;14.10

00000000000;15.6000000000000;13.6000000000000;14.5000000000000;1

3.7920000000000;11.6000000000000;12.1000000000000;14.700000000000

0;12.1000000000000;14.6000000000000;12.2000000000000;14.500000000

0000;14.2000000000000;13.4000000000000;13.3000000000000;14.100000

0000000;12.7000000000000;13.9000000000000;12;11.5000000000000;14.2

000000000000;14.4000000000000;14.1000000000000;15.1000000000000;1

2.2000000000000;13.3000000000000;14.1000000000000;13.30000000000

00;11.8000000000000;14.2000000000000;12.1000000000000;14.32000000

00000;13;14.2000000000000;13.9000000000000;12.9890000000000;13.400

0000000000;14.7000000000000;14.4250000000000;13.1000000000000;14.

3000000000000;14.3700000000000;14.3000000000000;13.3000000000000

;

14.6300000000000;14.5000000000000;13.6000000000000;14.2000000000

000;14.6000000000000;14.4540000000000;12.1000000000000;14.1800000

000000;14.2400000000000;11.9000000000000;14.4000000000000;12.4000

000000000;13.3000000000000;14.4000000000000;14.5000000000000];

plot(long, lat, '.b', 'MarkerSize', 12);

This would plot me a sample of my study region containing 100 events. What I'd like, is an algorithm that could divide the area into a uniform grid points, with intervals of 0.5 km in the horizontal coordiatnes, and then plot a 6 km radius circle around each grid point, to calculate how many events lie within it, and how far each of them is from the circle's center, i.e. the grid point.

[darova]

Do you have a better format data? I can't use this

Please use attachemnt button for long data

[Harith Al Kubaisy]

Sorry about that, I have attached the excel sheet with longitude and latitude in the original post.

However, for simplicity purposes, you can use this smaller dataset:

long = [57;52.1;50;52;43.5;50;53.5;52;43.2]

lat = [16;15.7;13.5;14.5;11.8;13.5;14.5;14;11.7]

[darova]

• What I'd like, is an algorithm that could divide the area into a uniform grid points, with intervals of 0.5 km in the horizontal coordiatnes

You have lat and long but want to draw grid 0.5km. Can you explain more?

Is this correct?

[Harith Al Kubaisy]

Yes, I'd like the grid points to be spaced out 0.5 km from each other, which equates to 0.004 in the decimal degrees.

The figure you put is exactly what I want. Circle around every grid point, prefrably around the cell's center, but this will do also. Thank you for keeping up with my questions!

Answer 1

[darova]

Here is how i see solution for this question

clc,clear

long = [57;52.1;50;52;43.5;50;53.5;52;43.2];
lat = [16;15.7;13.5;14.5;11.8;13.5;14.5;14;11.7];
step = 1;
x = min(lat):step:max(lat);
y = min(long):step:max(long);
[X,Y] = meshgrid(x,y);              % create grid
    
R = 1;                              % circle radius

plot(lat,long,'.r')                 % plot data
hold on
for i = 1:length(x)
    viscircles([x(i),y(1)], R)      % plot circle
    ix = (lat-x(i)).^2 + (long-y(1)).^2 < R^2;
    plot(lat(ix),long(ix),'ob')     % plot data inside circle (if exists)
end
plot(X,Y,'k',X',Y','k')             % plot grid
hold off
view(2)

[Harith Al Kubaisy]

You've been extremely helpful, thank you. For the most part, this does capture my task. However, the plot produced by the code draws circles around the lowermost grid points, i.e. the five points at Y ~ 43. How would I be able to draw it for the entire grid.

[Harith Al Kubaisy]

Moreover, the point Y = 57, X = 16 is plotted outside grid.

[Harith Al Kubaisy]

Ok, a few modifications (code below) made it possible to plot the circles around every grid point. The last question would be, how to incorporate a way to count data points inside every circle, and the distance between each data point and the circle's center?

clc,clear

long = [57;52.1;50;52;43.5;50;53.5;52;43.2];

lat = [16;15.7;13.5;14.5;11.8;13.5;14.5;14;11.7];

step = 1;

x = min(long):step:max(long);

y = min(lat):step:max(lat);

[X,Y] = meshgrid(x,y); % create grid

R = 1; % circle radius

plot(long,lat,'.r','MarkerSize',6) % plot data

hold on

for j = 1:length(y)

for i =1:length(x)

viscircles([x(i),y(j)], R) % plot circle

ix = (lat-x(i)).^2 + (long-y(1)).^2 < R^2;

plot(lat(ix),long(ix),'ob') % plot data inside circle (if exists)

end

end

plot(X,Y,'k',X',Y','k') % plot grid

hold off

view(2)

axis equal

[darova]

• how to incorporate a way to count data points inside every circle, and the distance between each data point and the circle's center?

It's in these lines

    ix = (lat-x(i)).^2 + (long-y(1)).^2 < R^2;
    plot(lat(ix),long(ix),'ob')     % plot data inside circle (if exists)

Just add

if ~isempty(sum(ix))
    N(i,j) = sum(ix);               % number of points inside (i,j) circle
else
    N(i,j) = 0;
end

[Harith Al Kubaisy]

Thank you, sorry for troubling you constantly, but I don't completely understand how to extract the needed information from this.

Also, would it be possible to incorporate the haversine distance equation to calculate the distance between each grid point (Circle's center) and every data event that lies within that given circle?

[darova]

• I don't completely understand how to extract the needed information from this.

Please be more specific.

• Also, would it be possible to incorporate the haversine distance equation to calculate the distance between each grid point (Circle's center) and every data event that lies within that given circle?

I think so. What problems are you facing?

[Harith Al Kubaisy]

Sorry for my poor formulation. I hope this will do a better job.

1. I'm not able to retrieve which data points (i.e. earthquake event) lie within a given circle from the matrix results. The matrices i and j return a count value. Is there a way to identify which data event is inside the circle (i.e. giving the data points an I.D. string value to be able to identify them later on)?
2. My problem is incorporating it into automatically calculating distances between data points that lie within a circle and the center of that given circle. Hence, ignoring circles which contain no data points within them.

[darova]

See attached script. I calculated Euclidean distance

D - distance point-circle (if point is inside)

N - number of points inside (i,j) circle

• Also, would it be possible to incorporate the haversine distance equation

How to calculate it?

[Harith Al Kubaisy]

This is brilliant, you're just natural. Thank you yet again. I'm now able to get counts of data points within a given circle. Is it possible to add an identifier field, such an ID number to distuinguish the events, and be able to se that in the results? I need that for my analysis, since each earthquake (data point) has a different magnitude. This would also enable to me to identify the distance between specific entries. Since the results of D and N are random at the moment (I can't tell which is which).

[darova]

• Is it possible to add an identifier field, such an ID number to distuinguish the events, and be able to se that in the results?

Can you show how you imagine this?

[Harith Al Kubaisy]

If I currently select data on the plot, and click on it, I would get X and Y values, which are the longitude/latitude values for that event. I'm suggesting to use a 3xM matrix for the data point such (longitude, latitude, ID no.)?

[darova]

I added these lines

ID = lat*nan;
for i = 1:length(x), for j = 1:length(y)
        % code
        ID(ix) = sub2ind(size(N),i,j);  % write ID
end,end

to call it

>> D(3)         % distance of 3d point to closest earthquake
ans =
    0.6800
>> ID(3)        % ID of eartquake
ans =
    11
>> N(ID(3))     % number of points inside earthquake
ans =
     2
>> D(1)         % 'NaN' means that is outside earthquake
ans =
   NaN

sub2ind link

[Harith Al Kubaisy]

Thank you. I'm running the code with a 972x1 double file for both longitude (X) and latitude values (Y) as the input data points. Note that I'm changing the plot commmand to plot(long,lat,'.r') in my code. The code runs slow with the viscircles function, since my analysis requires a very small step size (0.009) and radius (0.05). Hence I'm likely running into memory problems which is likely causing MATLAB to shut down before it finishes calculating, or run out of charge.

When I put a comment % on the second plot circle command viscircles([x(i),y(j)], R,'edgecolor','r'); and then run the code, it finishes much quicker, but a black rectangle covers the plot area and the results seem unrealistic, as I'm getting NaN for all ID values, and N is a 794x306 double full of 0's.

Any idea on how to go farther?

One idea I have is using an external hard drive to see if the memory issue is the culprit. The code seems to have finished calculating now (with the second vicircles command still on) because the Figure window popped-up now but it is still empty, and I keep getting notifications from my MacBook that my disk is almost full.

[Harith Al Kubaisy]

I just ran the code with a smaller input matrix (51x1 double for both long and lat files), and the code finished up plotting the figure in a reasonable time, with both viscircles commands on. However, the N file is a 131x217 double matrix of 0's again, and D and ID are both a 51x1 double of NaN.

Is something of the following nature more feasible for my task?

---------

hold on

% Get the values from the matrix

point_id=Matrix1(:,1);

x=Matrix1(:,2);

y=Matrix1(:,3);

% Set the text to show up

for i=1:length(x)

text(x(i),y(i), point_id)

end

---------

Moreover, to my current understanding, the currnet code calculates the distance between circle center and the closest earthquake inside the circle. However, there would be circles that contain more than one earthquake, and I would the distance of each of them from the point.

[Harith Al Kubaisy]

My idea about the Haversine distance calculation followed this logic: since every defined grid point has X and Y coordinates, and since every earthquake has the same, wouldn't be easier to implement a loop the utilizes MATLAB's built in function for calculating distance between long1lat1 & long2lat2, and then transform it from degrees to kilometers:

deg2km(distance(lat1, lon1, lat2, lon2))

In this case, one set of the long,lat would be the grid point, and the other any earthquake that lines within the circle around that grid point.

[Harith Al Kubaisy]

Ok, I think I made the stupid mistake of not changing some necessary variables since you're plotting lat, long while I need it to be long, lat. After minor modifications, the code currently reads like this:

clc,clear
load lat4
load long4
lat = (lat4);
long = (long4);
step = 0.009;
x = min(long):step:max(long);
y = min(lat):step:max(lat);
[X,Y] = ndgrid(x,y); % create grid
R = 0.05; % circle radius
plot(long,lat,'.r') % plot data
N = X*0; % number of points inside each circle
D = long*nan; % distance to closest circle
ID = long*nan;
hold on
for i = 1:length(x)
for j = 1:length(y)
D2 = (long-x(i)).^2 + (lat-y(j)).^2;
ix = D2 < R^2; % find smaller distances
D(ix) = D2(ix); % write smaller distances
ID(ix) = sub2ind(size(N),i,j); % write ID
N(i,j) = sum(ix); % number of points inside (i,j) circle
plot(long(ix),lat(ix),'ob') % plot data inside circle (if exists)
if sum(ix) % if there are points inside
viscircles([x(i),y(j)], R,'edgeColor','g'); % plot circle
else
viscircles([x(i),y(j)], R,'edgecolor','r'); % plot circle
end
end
end
plot(X,Y,'k',X',Y','k') % plot grid
hold off
axis equal

The results of the N, ID, D and D2 matrices make sense now. However, can you please still consider my inqueiries of the previous issues about the text ID and Haversine distances.

[darova]

To make your code faster use plot instead of viscircles. BY the way, shouldn't your grid look like following?

% elev = latitude
% azim = longitude
lat = -40:5:40;             % in degrees
long = -40:5:40;            % in degrees
lat = deg2rad(lat);         % convert to radians
long = deg2rad(long);
[LONG,LAT] = meshgrid(long,lat);
[xx,yy,zz] = sph2cart(LONG,LAT,earthRadius);
t = linspace(0,2*pi,20);    % 20 points for one circle
R = earthRadius/50;         % radius of circle
[x0,y0] = pol2cart(t,R);    % 20 points for circle
X = bsxfun(@plus,xx(:),x0*0);
Y = bsxfun(@plus,yy(:),x0);
Z = bsxfun(@plus,zz(:),y0);
plot3(X',Y',Z','r')
axis equal

To plot 'green' circle inside for loop:

if sum(ix) % if there are points inside
%     viscircles([x(i),y(j)], R,'edgeColor','g'); % plot circle
    plot(x0+x(i),y0+y(i),'g')
% else
%     viscircles([x(i),y(j)], R,'edgecolor','r'); % plot circle
end

• About Haversine formula

If im correct about your mesh (lat,long should be converted to spherical system to plot)

Then distance between points will be just:

D2 = (xdata-xx(i,j)).^2 + (ydata-yy(i,j)).^2; + (zdata-zz(i,j)).^2

• NaN values

If i understood you correctly: you don't like too many NaN values. So just

ID(isnan(ID)) = [];         % remove NaN
D(isnan(D)) = [];           % remove NaN

[Harith Kubaisy]

As ever, I'm grateful for your valuble contributions. Thank you a lot. I had problems signing in to my original account, so I made a new one since I have some more issues:

1- "BY the way, shouldn't your grid look like following?"

Yes and no. I'm trying to replicate the results made by a journal paper, in which they made a regional 2-D scan over a given study area, creating a 2-D grid at 0.5 km intervals in the horizontal coordinates, with a radius of 6 km around each grid point. My analysis calls for 1 km intervals, but with the same search radius. However, the sphere plot could come in handy as well.

------------------------------------------------------------------------------------------------------

2- For the following lines

t = linspace(0,2*pi,20);    % 20 points for one circle
R = earthRadius/50;         % radius of circle

I don't understand the use of 20 points. Is it so, that you are creating arbitrary points for the sake of demonstration, which would be my earthquake data points?

As for the radius, should the value in the denominator depend on my desired radius? What is the value of earthRadius in this case?

------------------------------------------------------------------------------------------------------

3- "To plot 'green' circle inside for loop:"

How is the for loop implemented with this modified code? Would it be:

for i=length(xx)
    for j=1:length(yy)

------------------------------------------------------------------------------------------------------

4- " Then distance between points will be just:

D2 = (xdata-xx(i,j)).^2 + (ydata-yy(i,j)).^2; + (zdata-zz(i,j)).^2

• NaN values

If i understood you correctly: you don't like too many NaN values. So just "

It's not that I don't like too many NaN values, it's that they are NaN that I don't understand. These should be numerical values for distance, or am I mistaken?

Also, would this calculate the distance between EVERY data point inside the circle, and the reference grid point?

Thank you again.

[darova]

• When use plot() you can control number of points. Can draw circle in 20 points or 200

t = linspace(0,2*pi,20);    % 20 points for one circle

• It's not that I don't like too many NaN values, it's that they are NaN that I don't understand. These should be numerical values for distance, or am I mistaken?

Attach your data

• Also, would this calculate the distance between EVERY data point inside the circle, and the reference grid point?

No actually. I changed the script. Now it creates a table of distances like following

See attached script

[Harith Kubaisy]

Attached are my data files. I will try out the new script.

[darova]

Looks ok for me

[Harith Kubaisy]

Thank you, the code runs instantly now, thank you. I still have some issues, so kindly bear with me like you have been.

1. When I set the step size to 0.009, which is the decimal degrees equivalent for 1 km for the grid size, I still don't get as many grids. Also, the distribution of the earthquakes in my figure is the mirror of the image you attached:
2. Why are the axis labelled in this manner? Shouldn't be these the longitude/latitude values?
3. Suppose I take the most north-western grid in your figure, the first one at the top from the left, it seems to only have one earthquake inside the circle. Which one would that be in the tables?

[darova]

• When I set the step size to 0.009, which is the decimal degrees equivalent for 1 km for the grid size, I still don't get as many grids. Also, the distribution of the earthquakes in my figure is the mirror of the image you attached:

Try to change step to 0.009*pi/180 (in radians). I think it's to small step, you'll get too many circles. My computer even couldn't handle it

• Why are the axis labelled in this manner? Shouldn't be these the longitude/latitude values?

Agree. I found better way to visualize (worldmap function)

• Suppose I take the most north-western grid in your figure, the first one at the top from the left, it seems to only have one earthquake inside the circle. Which one would that be in the tables?

Maybe better be to pick those points. I added additional section at the end of the code (just press Ctrl+Enter to run)

See script

[Harith Kubaisy]

Thank you for the modification, the axis makes sense now.

You're right, the step size is too small, but that's the grid size we need for the case study (1 km in the horizontal coordinates). Indeed I got an error when I changed the step size to 0.009*pi/180. How would you recommend going about this issue? My complete earthquake set is even bigger (2027 earthquake events i.e. data points). Would creating a datastore and integrating a tall array make a difference?

Also, when I tried your tip of moving the cursor to line 41 and pressing ctrl+enter, nothing happened. Unless I'm missing something here. I'm able to select specific circles from the plot via the mouse cursor, and it returns a count answer, which i persume is the number of data points inside the circle. However, the ID table is still just zeros (972x1092 logical), and thus I'm unable to deduce which datapoint this is, and in which circle it is, and the distances between the datapoints and the circle center. My thinking is, it would be more intuitive, to give each grid point, and circle, a fixed label, whether it be text or numerical.

(p.s. I'm really, really grateful for all your contributions. You're making me enjoy MATLAB more and more. Thank you).

[Harith Kubaisy]

Hello again, I hope you're not fed up with our discussion. I'm trying something of the following nature, but I don't know if logic is on my side:

(Note: completenondecl.csv being a spreadsheet cointaining the full list of earthquakes (2028 events), the first column is longitude and the secodn is latitude. )

clc,clear

% load lat5.mat
% load long5.mat
ds = datastore('completenondecl.csv','TreatAsMissing','NA','MissingValue',0);
T = tall(ds);
lat0 = T(:,2);
long0 = T(:,1);
step = 0.009*pi/180;        % in degrees
long = min(long0):step:max(long0);
lat = min(lat0):step:max(lat0);
[xdata,ydata,zdata] = sph2cart(long0*pi/180,lat0*pi/180,earthRadius);
% elev = latitude
% azim = longitude
[LONG,LAT] = ndgrid(long,lat);
[xx,yy,zz] = sph2cart(LONG*pi/180,LAT*pi/180,earthRadius);
t = linspace(0,2*pi,20);    % 20 points for one circle
R = 6000;                   % radius of a circle in [m]
[x0,y0] = pol2cart(t,R/earthRadius*180/pi);

I get the following error

Error using tall/min (line 17)
Invalid data type. First argument must be numeric or logical.

Error in darova7 (line 10)
long = min(long0):step:max(long0);

[darova]

You declared wrong step

Can you attach the data?

[Harith Kubaisy]

Here you go.

[darova]

I don't have datastore. So i used importdata

Here is another approach

T = importdata('completenondecl.csv');
lat0 = T(:,2);
long0 = T(:,1);
step = 0.009;
lat  = round((lat0 -min(lat0))/step)*step+min(lat0);
long = round((long0-min(long0))/step)*step+min(long0);

t = linspace(0,2*pi,20);    % 20 points for one circle
R = 1000;                   % radius of a circle in [m]
[x0,y0] = pol2cart(t,R/earthRadius*180/pi);
LATC = bsxfun(@plus,lat(:),x0)';
LONGC = bsxfun(@plus,long(:),y0)';
%
worldmap([min(lat0) max(lat0)], [min(long0) max(long0)])
plotm(lat0,long0,'.r')
hold on
plotm(LATC,LONGC,'b')
hold off

But i don't how precise it can be because of angle

[Harith Kubaisy]

Thank you darova, but the code is drawing circles around the earthquakes themselves, not around the arbitrary grid points we want to create. Also, what are the variables x0 and y0?

[darova]

• the code is drawing circles around the earthquakes themselves

Please add this line to your code

plotm(lat,long,'.g')

You will see circle centers

[Harith Kubaisy]

Noted. See attached figure. Something is still off. The grid should have way too many small squares if 0.009 is treated as decimal degrees (1 km increment). There is no need for plotting a circle around the earthquake events themselves. Only around the grid points, to eventually calculate the kernel density distribution of the area. So what is needed is: the count of earthquakes within 6 km radius of a grid point, which are in 1 km (approx = 0.009 decimal degrees) increments in the horizontal coordinates, and the distance of each invidual earthquake from the grid point center. That's why an ID is necessary for each earthquake and distance, as the earthquakes have different sizes, and I would need to account for the magnitude of each earthquake in a later part of the analysis. Said magnitude is also a column of the data spreadsheet, namely the 6th column.

[darova]

• The grid should have way too many small squares if 0.009 is treated as decimal degrees (1 km increment)

Use meshgrid and previous scripts to plot/draw grid

• So what is needed is: the count of earthquakes within 6 km radius of a grid point

It's a lot of computations. I have an idea how to simplify it. But of couse for cost of precision

Snap each earthquake data to grid point (lost of precision)

A then just use this mask to find neighbour points/grid/earthquake (R=6)

[Harith Kubaisy]

I don't understand how you made the snap, and the transition to the second figure (the mask).

[darova]

Another example: you have lat and long as following

lat = [1 1.1 2.3 5.1 3.2];
long = [5 6.2 1.3 2.1 6.2];
    % round data
lat1 = round(lat);
long1 = round(long);
lat1 =
     1     1     2     5     3
long1 =
     5     6     1     2     6

Now data is of integer format (1,2,3,...) so those value can be positions in matrix

Z = zeros(7,7);         % matrix of zeros 7x7
Z(1,5) = 1;             % put data inside using new coordinates
Z(1,6) = 1;
Z(2,1) = 1;

Is it clearer now?

[Harith Kubaisy]

Ok I get you now, but this loss of precision is quiet costly, I don't think the analysis will make sense if this approach is used. I think it's better to keep using the previous scripts since it was going in the right direction, namely the one from 21st March. I'm sorry if I'm taking your asssitance for granted, but I'm quiet sure that a few tweaks could solve this.

So from that I want to ask: what is the smallest step size that our computers can handle with that code? I'm discussing this with my supervisor; using the smallest step size possible, while increasing the circle radius to make sure all data points are accounted for.

[darova]

• I think it's better to keep using the previous scripts since it was going in the right direction

I understand that. In the previous script pdist2 function tries to create matrix. Matrix represents all combinations of distances between points. The data you attached consists of 2028 points. Then you want grid with step 0.009 (its 743x1812 = 1 346 316 points).

So totally pdist2 tries to create matrix of size 2028 x 1 346 316 (its big). That is why it failed.

a = zeros(2028,1345573);            % my computer has difficulties to create such matrix
% 1 value of double is 8 bytes

Why do you need to calculate distance from every point to every grid point? Because you don't know how far each point form current grid.

By the way, can you explain more about kernel distribution? You want number of earthquakes for each grid point in 6km radius. What do you want to do with them?

• So from that I want to ask: what is the smallest step size that our computers can handle with that code?

It depends on rounding and size of matrix you can afford. Earlier example:

lat = [1 1.12 2.33 5.13 3.22];
long = [5 6.21 1.32 2.15 6.26];
step = 0.1;
    % create indices for matrix
ilat = round(lat/step);
ilong = round(long/step);
ilat =
    10    11    23    51    32
ilong =
    50    62    13    21    63

You know what i mean?

[Harith Kubaisy]

• Why do you need to calculate distance from every point to every grid point? Because you don't know how far each point form current grid.

I don't. I actually only need the distance between earthquakes that are inside any given circle, and that specific circle's centre, i.e. that specific grid point. That's why I got confused with these tables.

• By the way, can you explain more about kernel distribution? You want number of earthquakes for each grid point in 6km radius. What do you want to do with them?

I'm going to copy-paste a part of a journal paper that we're trying to emulate:

"For each grid point, both parameters are calculated utilising all recorded events within a 6 km radius. The parameters are calculated based on the kernel density estimation as an approach to obtain the spatial distribution through a probability density function, using the distance to weight each event from a reference point (each grid point, the common centre of its adjacent events). This circular-shape-based approach prevents any directional bias. The earthquake kernel density parameter, ρNk, is calculated by counting all the weighted events within a 6 km radius from each grid point, dividing their sum by the sampler area (π r 2 ) and normalising by the duration of the earthquake catalogue:

where N is the total number of events within the radius r, d(n) is the distance between an event n and the circle centre, σ is the standard deviation of the Gaussian function and T is the duration of the earthquake catalogue. Units are events per squared kilometre per year (events km−2 yr−1).

The M0 kernel density parameter, ρM0 k , is obtained by first calculating the seismic moment released by each event separately, using the empirical relation between M0 and ML, as obtained by Shapira and Hofstetter (1993) after converting units from dyne×cm to N ×m:

log[M0]=10+1.3ML. (2)

Secondly, each amount of energy is weighted according to the distance of the corresponding event from the circle centre (like the calculation of the earthquake kernel density). Then, we sum the weighted M0 released from all the events within a 6km radius, divide the sum by the circle area (πr2) and normalise by the duration of the catalogue:

"

Reference:

M. Sharon et al. (2020): Assessment of seismic sources and capable faults through hierarchic tectonic criteria 129

[Harith Kubaisy]

For a smaller sample (N=51) of my complete data set, this block of code manages to work reasonably fast, even with the desired small step length:

clc,clear
load lat4
load long4
lat = (lat4);
long = (long4);
step = 0.009;
x = min(long):step:max(long);
y = min(lat):step:max(lat);
[X,Y]= ndgrid(x,y); % creates grid
R = 0.05; % circle radius
plot(long,lat,'.r') % plot data
N = X*0; % no. of points inside each circle 
D = long*nan; % distance to closest circle 
ID = long*nan; 
hold on
for i= 1:length(x) 
    for j=1:length(y)
      D2 = (long-x(i)).^2 + (lat-y(j)).^2;
        ix = D2 < R^2; % find smaller distances
        D(ix) = D2(ix); % write smaller distances
        ID(ix) = sub2ind(size(N),i,j); %write ID
        N(i,j) = sum(ix); % no. of points inside (i,j) circle
        plot(long(ix),lat(ix),'ob') %plot data inside circle (if exists)
        if sum(ix) % if there are points inside
            viscircles([x(i),y(j)], R,'edgeColor','g'); % plot circle
        else
%           viscircles([x(i),y(j)], R,'edgecolor','r'); % plot circle
        end
    end
end
hold off
axis equal 

However, I'm still not able to deduce which of the resulting distances are relevant for my analysis.

So to summarize what I need:

1. Plot earthquakes on a grid with intervals of 1 km and plot a circle around each grid point with a radius of 6 km.
2. How many earthquakes are in each circle.
3. What is the distance of each earthquake from the circles from step 2.
4. Account for the magnitude of each earthquake from step 3.

That's why I was suggesting using some sort of identifier string. I need to distuingish between the data for the second part of the analysis (the seismic moment distribution) since every earthquake has a different magnitude.

Attached here also is a spreadsheet with a new column in the end, numbered from 1 to 2028 that should serve as an ID for the earthquakes.

I'm sorry if I'm being redundant here, and as ever, thank you for your contribution. You're a true star.

[darova]

So you don't actually need to store all data. You can do calculation during for loop

[Harith Kubaisy]

Bear with me, but what is sig^2 , and T?

Also the p=p/pi*R^2 , is this replacing anything?

[darova]

• Bear with me, but what is sig^2 , and T?

Haha, you tell me, those are your formulas

• Also the p=p/pi*R^2 , is this replacing anything?

Just dividing all data by circle area (didn't wanted to do this inside for loop)

More operations outside for loop - faster code

R = 6;
tic
for i = 1:1e6               % do million operations
    p = 1/(pi*R^2);
end
t1 = toc
A = pi*R^2;                 % calculate area once
tic
for i = 1:1e6               % do million operations
    p = 1/A;
end
t2 = toc
sprintf('second loop faster %f times',t1/t2)

[Harith Kubaisy]

• Haha, you tell me, those are your formulas

D'oh. I completely failed to realize that it is indeed the same formula. In that case, T which is my duration (from the year 1916 to 2019) should be defined before the loop. Right?

• Just dividing all data by circle area (didn't wanted to do this inside for loop)

Makes sense.

[darova]

• should be defined before the loop. Right?

Everything that can be outside loop - should be outside

[Harith Kubaisy]

I'm getting strange plots with the current script. I'm only taking into account the first 50 earthquakes of the full catalog.

clc,clear
% load lat4
% load long4
AA = importdata('completenondecl3.csv');
lat = AA(1:50, 2);
long = AA(1:50, 1);
step = 0.009;
x = min(long):step:max(long);
y = min(lat):step:max(lat);
[X,Y]= ndgrid(x,y); % creates grid
R = 0.05; % circle radius
plot(long,lat,'.r') % plot data
P = X*0; % no. of points inside each circle 
T = 2019-1916+1;
sig = 2;
 
hold on
for i= 1:length(x) 
    for j=1:length(y)
      D2 = (long-x(i)).^2 + (lat-y(j)).^2;
        ix = D2 < R^2; % find smaller distances
        D(ix) = D2(ix); % write smaller distances
        P(i,j) = sum(-exp(D(ix).^2/2/sig^2))/T; %write ID
%         N(i,j) = sum(ix); % no. of points inside (i,j) circle
        plot(long(ix),lat(ix),'ob') %plot data inside circle (if exists)
        if sum(ix) % if there are points inside
            viscircles([x(i),y(j)], R,'edgeColor','g'); % plot circle
        else
%           viscircles([x(i),y(j)], R,'edgecolor','r'); % plot circle
        end
    end
end
plot(X,Y,X',Y')
P = P/(pi*R^2);
imagesc(P)
hold off
axis equal 

[darova]

• You can trust me

• check units

• Look what viscircles is doing

[Harith Kubaisy]

• check units

To be honest, I'm still confused with the definition of ix in the line ( ix = D2 < R^2)

According to the original equation, it should be the distance between the earthquake and the circle's center, hence km.

• Look what viscircles is doing

I realize viscircles is plotting something. But take a look at the following plot, for which I used 51 earthquakes only for the plot. Look at the axis. The X-Y axes are the longitude/latitude values and should not reach these values.

Also, what is the purpose of the following command:

plot(long(ix),lat(ix),'ob') %plot data inside circle (if exists)

Plotting circles around the earthquakes themselves is not necessary.

(Also, thank you.)

[darova]

It was a joke actually. Look at the minus sign, it should be inside exp(). Distance is already squared (no need for 2 power)

• UNITS

• X-Y plane: specifiy axis

imagesc(x,y,P)          % no axis before, so X and Y just rows and columns

[Harith Kubaisy]

• It was a joke actually. Look at the minus sign, it should be inside exp(). Distance is already squared (no need for 2 power)

Pardon me, this task is stripping away my sense of humor. Like i said in my previous posts, I'd assume D(ix) is in degrees as well in this script. Right?

      D2 = (long-x(i)).^2 + (lat-y(j)).^2;
        ix = D2 < R^2; % find smaller distances
        D(ix) = D2(ix); % write smaller distances

I've read/wrote those lines so much now that I'm confused about them, and I must have driven you crazy by now (just in case if CoVid-19 did not).

But in all seriousness, I'm still confused about:

1. This line:

        plot(long(ix),lat(ix),'ob') %plot data inside circle (if exists)

2. Which variable defines the distance between earthquakes in a given circle, and that same circle's center (the grid point)? I'm assuming it is D(ix)? However, the ix is a 51x1 logical column of 0's only.

The plot looks like this now:

The axes makes sense again, however:

3. What do the colors represent? and why is the majority of the plot overlain by yellow?

4. The result of P is a 131x217 double. I assume this is the result of the summation equation at every grid point in the plot. The following screenshot represents some of the values in that double. They're negative.

Sorry if I'm making you crazy. Thanks.

[Harith Kubaisy]

Also, seems that the Statistics Toolbox has a function ksdensity (https://de.mathworks.com/help/stats/ksdensity.html)

Is it helpful in this quest?

[darova]

• I'd assume D(ix) is in degrees as well in this script. Right?

Yes. Degrees of power 2 ( )

• 1. This line:
• plot(long(ix),lat(ix),'ob') %plot data inside circle (if exists)

It was one of the one first scripts. I didn't know that you don't need circles around points themselves. You can just remove it. My script is only suggestion, you decide what should be in the result

• 2. Which variable defines the distance between earthquakes in a given circle, and that same circle's center (the grid point)? I'm assuming it is D(ix)? However, the ix is a 51x1 logical column of 0's only

Yes D is array of distances for each earthquake (not circle/grid) (*not distances but degrees^2, they should be converted into km to be distances). ix is indicator (remember you loop through grid points, if grid point doesn't have closest earthquake - all zeros)

You asked me previously if D array can store more distances (earthquake can be in two grid points). The answer: in this case it stores only last grid point.

• 3. What do the colors represent? and why is the majority of the plot overlain by yellow?

It is the question to you: what you are calculating?

• 4. ... The following screenshot represents some of the values in that double. They're negative.

Did you replace minus sign inside exp() function?

sum(exp(-D/2/sig^2))/T; % write ID

If it's your first time with programming and MATLAB you need some time to understand all those things. Programming is like math or playing on guitar

[Harith Kubaisy]

It is my first time working independetly on MATLAB, and like you can see, it's not so independent after all. I have basic experience with the program.

• It is the question to you: what you are calculating?

My task is calculating the earthquake kernel density and seismic moment kernel density. The end product of such an analysis is a map like this:

(The study used a 0.5 km interval grid and 6 km radius)

So I'd like to make sense of the plot generated by MATLAB since there is no scale/reference to the colors.

• Did you replace minus sign inside exp() function?

I just have. Now the plot looks like this:

[darova]

Try to scale data (to km)

scale = (pi/180*earthRadius/1000)^2;
% code
D(ix) = scale*D2(ix); % write smaller distances
P(i,j) = sum(exp(-D(ix)/2/sig^2))/T; %write ID

An idea to speed up your code: work only in specific area

[Harith Kubaisy]

Darova, I don't know who you are, but I'm so grateful that people like you exist. Thank you.

I've used the scale idea. This is the my current script:

clc,clear
load lat4
load long4
AA = importdata('completenondecl3.csv');
% lat = AA(1:50, 2);
% long = AA(1:50, 1);
lat = (lat4);
long = (long4);
step = 0.009;
x = min(long):step:max(long);
y = min(lat):step:max(lat);
[X,Y]= ndgrid(x,y); % creates grid
R = 0.05; % circle radius
plot(long,lat,'.r') % plot data
P = X*0; % no. of points inside each circle 
T = 2019-1916+1;
sig = 2;
scale = (pi/180*earthRadius/1000)^2;
hold on
for i= 1:length(x) 
    for j=1:length(y)
      D2 = (long-x(i)).^2 + (lat-y(j)).^2;
        ix = D2 < R^2; % find smaller distances
        D(ix) = scale*D2(ix); % write smaller distances
        P(i,j) = sum(exp(-D(ix)/2/sig^2))/T; %write ID
%         N(i,j) = sum(ix); % no. of points inside (i,j) circle
%         plot(long(ix),lat(ix),'ob') %plot data inside circle (if exists)
        if sum(ix) % if there are points inside
            viscircles([x(i),y(j)], R,'edgeColor','g'); % plot circle
        else
%           viscircles([x(i),y(j)], R,'edgecolor','r'); % plot circle
        end
    end
end
P = P/(pi*6.^2);
imagesc(x,y,P)
hold off
axis equal 

I agree about working with a smaller area. In the script I only chose 51 earthquakes (long4 and lat4 are 51x1 doubles), and the calculation was fast. I can break it up by areas of similar longitude/latitude values. I don't mind repeating the process 40 times : )

Also, you can see that I edited the final equation, since R = 0.05 (degrees). I chose 6, since the rest of the values are already scaled to kilometers. Is it correct going about it this way?

[darova]

Works ok?

You don't use earthquake magnitude? If you take a closer look at image, do you any difference between points? they are look similar

[Harith Kubaisy]

• You don't use earthquake magnitude?

The earthquake magnitude is taken into account for the second part of the analysis, which is the Seismic Moment Kernel Density distribution. That will wait until I'm done successfully calculating the Earthquake Kernel Density distribution.

• If you take a closer look at image, do you any difference between points? they are look similar

What points do you mean?

--------------------------------------------------------------------------

Another question:

if I call this command

>> N(23,51)

ans =

     3

So I know that in that specific grid (or the circle around that grid to be precise), I have 3 earthquakes. How can I get the individual distance of each of them from that specific grid point, along with knowing which distance is for which earthquake? This is important for the second part of the analysis. A reminder that my data spreadsheet has a column in the end which should serve as an ID for each earthquake (1 to 2028).

[Harith Kubaisy]

Or, just a way to know which three earthquakes are the ones inside the circle, and then by knowing the coordinates of the grid point, I could calculate the distance between them using Haversine formula.

[darova]

• What points do you mean?

Most of them look the same

• if I call this command

Do you have coordinates of that grid point?

[Harith Kubaisy]

• Most of them look the same

I think bigger differences would be more visible when all data are plotted.

• Do you have coordinates of that grid point?

No. That's what I'm trying to figure out. By definiton, the coordinates should be there somewhere.

[Harith Kubaisy]

Hello Darova, I hope you can help me with this still lingering issue:

I managed to get a reasonable plot after some orientation errors. My supervisor thinks the map makes sense but we need to prove the math now. So I'm attaching the current script I'm using:

clc,clear
% load latfull
% load longfull
% load longmatrix
% load latmatrix
% load coord
load long5
load lat5
% AA = importdata('completenondecl3.csv');
long = (long5);
lat = (lat5);
step = 0.009;
x = min(long):step:max(long);
y = min(lat):step:max(lat);
[X,Y]= ndgrid(x,y); % creates grid
R = 0.05; % circle radius
plot(long,lat,'.r') % plot data
P = X*0; % no. of points inside each circle 
T = 2019-1916+1;
sig = 2;
scale = (pi/180*earthRadius/1000)^2;
hold on
for i= 1:length(x) 
    for j=1:length(y)
      D2 = (long-x(i)).^2 + (lat-y(j)).^2;
        ix = D2 < R^2; % find smaller distances
        D(ix) = scale*D2(ix); % write smaller distances
        P(i,j) = sum(exp(-D(ix)/2/sig^2))/T; %write ID
        N(i,j) = sum(ix); % no. of points inside (i,j) circle
%         plot(long(ix),lat(ix),'ob') %plot data inside circle (if exists)
        if sum(ix) % if there are points inside
            viscircles([x(i),y(j)], R,'edgeColor','g'); % plot circle
        else
%           viscircles([x(i),y(j)], R,'edgecolor','r'); % plot circle
        end
    end
end
p = P/(pi*6.^2);
G = rot90(flip(p),3);
imagesc(x,y,G);
colorbar
hold off
axis equal 

For simplicity, I attached two data files with ten earthquakes only to make it easier to control and compute.

To prove the calculations correctly I have chosen an arbitrary cell which is N(21,43) which has 2 earthquakes. Now I want to solve the density formula by hand, and for that I need the distance between each one of these earthquake to the grid point, which according to what I make from the figure, is at X = 42.0016 and Y = 11.8004. So if there is just a method to know which two earthquakes are the ones inside the circle, and then computing the distance between each will be easy by Haversine formula.

I need to be able to distuinguish between these earthquakes since the next part of the analysis is magnitude dependant. I'm still having a hard time with this.

Thank you.

[darova]

Hello!

• For simplicity, I attached two data files with ten earthquakes only to make it easier to control and compute.

It's very generous of you. THank you!

• I need the distance between each one of these earthquake to the grid point

Sure. You can use pdist2 again (note: distance between points in degrees)

Also MATLAB gives you some valuable warnings

And only one advice from me

• I need to be able to distuinguish between these earthquakes since the next part of the analysis is magnitude dependant. I'm still having a hard time with this.

I'll hold my firsts for you!

[Harith Kubaisy]

Thanks for your reply. The duplicates are two separate events that occured at the same coordinates. That's part of the problem (not distinguishing events). Thanks for the pdist2 suggestion, somehow I missed it from earlier, or did not understand how to use it. I will give the calculation a go soon and see how it goes.

• And only one advice from me

I'm assuming you mean that using the plot command would be better here?

[darova]

Yes, definitely

[Harith Kubaisy]

• Sure. You can use pdist2 again (note: distance between points in degrees)

Alright, but in doing the following command

pdist2([long lat],[X(21,43) Y(21,43)]) 

We get a list of distances between every earthquake, and this grid point. I want only the distances for the two earthquakes that are within the circle of N(21,43). Unless I'm missing out on something again..

[darova]

Sure. Use logical operators (or find)

D = pdist2([long lat],[X(21,43) Y(21,43)]);     % all distances
D(D<0.05)  % distances inside circle

[Harith Kubaisy]

Thank you a million! This definitely solves a major issue. Alright.

For cell N(21,43), which has two earthquakes within the circle, who are both 0.0016 degrees away from the grid point, which approximates to d(n) = 0.1776 km. T = 104 years. r = 6 km. sigma = 2.

Solving by hand for the eathquake kernel density of cell N(21,43) gives me a value of about 8.468e-05, while the array p from the script returns it as 1.6932e-04.

There is something off.

Scratch that. I forgot to multiply by two. My calculator reads 1.6936e-04 now. Thank you.

[Harith Kubaisy]

Hello Darova, I need to visualize the final density values (p) on a proper, projected geographic map. Is there a way to do this on MATLAB with the current script? If not, what is the best way to export the data to be able to visualize it on ArcMap 10.4?

Answer 2

[Cris LaPierre]

You can compute distances between two lat/lon coorinates using the Haversine formula. You can google it to find the equation. You can see an example of its implementation in MATLAB around the 6:50 mark in this video.

[Harith Al Kubaisy]

Thanks for the response. That I'm aware of, but this partially captures the task I'm facing. I need to create circles around arbitrary grid points, which are uniformly spaced at 0.5 km (0.09 decimal degrees) intervals in the horizontal coordinates. How would I go ahead establishing this grid and circles automatically, not having to write down coordinates for each circle?

[Harith Al Kubaisy]

I don't have the lat/lon coordinates for the grid points, because they are not defined yet. This is part of the question I'm asking: defining the grid points.