How to use pdist2 command for distance calculation?

I would like to use the command pdist2 in order to compute distance between two geographical coordinates. I want the haversine distance to be more specific. I know the distance(lat1,lon1,lat2,lon2) command is the suitable for this but I want to use pdist2 command. Could I modify in some manner pdist2 command?

Could you please help me?

Answer 1

[Torsten]

If you implement your own distance function for the haversine distance, you can use pdist2.

Answer 2

[Rik]

I don't expect the performance to be great, but you can use the option of specifying the distance calculation function.

D = pdist2(X,Y,@CalcDist)
function D2=CalcDist(ZI,ZJ)
% calculation of distance
%
% ZI is a 1-by-n vector containing a single observation.
% ZJ is an m2-by-n matrix containing multiple observations. distfun must
% accept a matrix ZJ with an arbitrary number of observations.
% D2 is an m2-by-1 vector of distances, and D2(k) is the distance between
% observations ZI and ZJ(k,:).
D2 = zeros(size(ZJ,1),1);
for k=1:numel(D2)
    % implement this calculation with distance(lat1,lon1,lat2,lon2)
    % its documentation states: "You can specify lat1, lon1, lat2, and lon2
    % using a combination of scalars and arrays, as long as the arrays are
    % of the same size. The function expands the scalar inputs to match the
    % size of the array inputs."
end
end

[Ivan Mich]

I tried

lat1 = 1;lon1 = 2;lat2 = 3;lon2 = 4; % edit Rik: added this to run the code
X=[lat1 lon1];
Y=[lat2 lon2];
D = (pdist2(X,Y,'@haversine'))

Error using statslib.internal.pdist2
Unrecognized 'DISTANCE' argument: @haversine.

Error in pdist2 (line 158)
[varargout{1:nargout}] = statslib.internal.pdist2(X,Y,varargin{:});

function [dist] = haversine(lat1,lon1,lat2,lon2)
dist = 2 * 6372.8 * asin(sqrt(sind((lat2-lat1)/2)^2 + cosd(lat1) * cosd(lat2) * sind((lon2 - lon1)/2)^2));
end

but no use..

[Stephen23]

lat1 = 1;
lon1 = 2;
lat2 = 3;
lon2 = 4;

X = [lat1,lon1];
Y = [lat2,lon2];
D = pdist2(X,Y,@haversine)

D = 314.4918

Checking against your original code:

2 * 6372.8 * asin(sqrt(sind((lat2-lat1)/2)^2 + cosd(lat1) * cosd(lat2) * sind((lon2 - lon1)/2)^2))

ans = 314.4918

function out = haversine(m1,m2)
tmp = sind((m2-m1)./2).^2;
out = 2 * 6372.8 * asin(sqrt(tmp(:,1)+cosd(m1(:,1)).*cosd(m2(:,1)).*tmp(:,2)));
end

[Rik]

% Define the variables
lat1 = 1;lon1 = 2;lat2 = 3;lon2 = 4;
X=[lat1 lon1];
Y=[lat2 lon2];

try
    D = (pdist2(X,Y,'@haversine'))
catch ME,ThrowAsWarning(ME),end

Warning: Unrecognized 'DISTANCE' argument: @haversine.

I had expected that after more than 160 questions, you would be able to read an error message. It tells you that the argument is incorrect. That should prompt you to have a look at the examples in the documentation. That will show you that you shouldn't provide a char vector, but a function handle, so removing the quotes does the trick.

try
    D = (pdist2(X,Y,@haversine_version1))
catch ME,ThrowAsWarning(ME),end

Warning: Error evaluating distance function 'haversine_version1'.
because:
Not enough input arguments.

And now you see a different error message: the pdist2 function doesn't result in enough input arguments for your custom function. Again, looking at the documentation, you see that you need a different defnition. (see below for version2)

For this try, let's also remove the extra parentheses that don't actually do anything.

try
    D = pdist2(X,Y,@haversine_version2)
catch ME,ThrowAsWarning(ME),end

D = 314.4918

Good, this seems to work.

Here is version 1:

function [dist] = haversine_version1(lat1,lon1,lat2,lon2)
dist = 2 * 6372.8 * asin(sqrt(sind((lat2-lat1)/2)^2 + cosd(lat1) * cosd(lat2) * sind((lon2 - lon1)/2)^2));
end

And here is version 2:

function [dist] = haversine_version2(m1,m2)
lat1 = m1(:,1);
lat2 = m2(:,1);
lon1 = m1(:,2);
lon2 = m2(:,2);
dist = 2 * 6372.8 * asin(sqrt(sind((lat2-lat1)/2)^2 + cosd(lat1) * cosd(lat2) * sind((lon2 - lon1)/2)^2));
end

Now we are almost at the version Stephen suggested. We define a temporary variable that does most of the calculation already and makes the next line more compact (and let it fit on a single line). We also remove the square brackets that don't do anything.

function dist = haversine_version3(m1,m2)
tmp = sind((m2-m1)/2)^2;
dist = 2 * 6372.8 * asin(sqrt(tmp(:,1)+cosd(m1(:,1))*cosd(m2(:,1))*tmp(:,2)));
end
function out = haversine_Stephen(m1,m2)
tmp = sind((m2-m1)./2).^2;
out = 2 * 6372.8 * asin(sqrt(tmp(:,1)+cosd(m1(:,1)).*cosd(m2(:,1)).*tmp(:,2)));
end

You will notice that there are still differences between this third version and what Stephen suggested: several periods to ensure element-wise multiplication and division. That is required for array support. What's also missing is documention. That I will leave for you to write. Remember that code without comments is completely and utterly useless once it doesn't work.

function ThrowAsWarning(ME)
% This is just a helper function to show the error message as warning.
if ~isempty(ME.cause)
    warning('%s\n  because:\n%s',ME.message,ME.cause{1}.message)
else
    warning('%s',ME.message)
end
end