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

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]

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]

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..

Answer 1

[David Franco]

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)

radius = 6371; % Earth radius

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!