How to calculate the Euclidean distance beetwen all points of Latitude Longitude pairs?

53 views (last 30 days)
David Franco on 20 Apr 2020
Answered: David Franco on 21 Apr 2020
I have a 399 cities array with LON LAT coordinates (first column for the Longitudes), like the picture below.
How can I calculate the 399x399 matrix with all distances between this 399 cities?
I used pdist and squareform but the result are small number. Am I correct?
D = pdist(XY);
Z = squareform(D);
For example, the two first points (-50.3125 -23.3005; -48.9918 -24.6617) have a Euclidean distance between them of 216 km (see picture below).
Thank you!
David Franco on 21 Apr 2020
In time: I found that the values in my Z matrix is in fact the 2-norm (or Euclidean distance) between all my 399 points. But how can I convert them to km? Only with haversine formula?
Geoff Hayes on 21 Apr 2020
David - yes, Haversine is one way to get the distance between two latitude and longitude points. The Vincenty algorithm (see https://www.mathworks.com/matlabcentral/fileexchange/5379-geodetic-distance-on-wgs84-earth-ellipsoid) is another..

David Franco on 21 Apr 2020
function [d1,d2] = pos2dist(point1,point2)
% Distance:
% d1: distance in km based on Haversine formula
% d2: distance in km based on Pythagoras theorem
% Inputs:
% point1: lat lon of origin point [lat lon]
% point2: lat lon of destination point [lat lon]
%
% Outputs:
% d1: distance calculated by Haversine formula
% d2: distance calculated based on Pythagoran theorem
%
% Example 1, short distance:
% point1 = [-43 172];
% point2 = [-44 171];
% [d1 d2] = pos2dist(point1,point2)
% d1 =
% 137.365669065197 (km)
% d2 =
% 137.368179013869 (km)
%
% Example 2, longer distance:
% point1 = [-43 172];
% point2 = [20 -108];
% [d1 d2] = pos2dist(point1,point2)
% d1 =
% 10734.8931427602 (km)
% d2 =
% 31303.4535270825 (km)
lat1 = point1(1)*pi/180;
lat2 = point2(1)*pi/180;
lon1 = point1(2)*pi/180;
lon2 = point2(2)*pi/180;
deltaLat = lat2-lat1;
deltaLon = lon2-lon1;
a = sin((deltaLat)/2)^2 + cos(lat1)*cos(lat2) * sin(deltaLon/2)^2;
c = 2*atan2(sqrt(a),sqrt(1-a));
x = deltaLon*cos((lat1+lat2)/2);
y = deltaLat;
d1 = radius*c; % Haversine distance
d2 = radius*sqrt(x*x + y*y); % Pythagoran distance
end
Thanks Geoff!

Categories

Find more on Geometric Geodesy in Help Center and File Exchange

R2019b

Community Treasure Hunt

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

Start Hunting!