How to use pdist2 command for distance calculation?

23 views (last 30 days)
Ivan Mich on 12 Oct 2023
Commented: Rik on 13 Oct 2023
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?

Torsten on 12 Oct 2023
If you implement your own distance function for the haversine distance, you can use pdist2.

Rik on 12 Oct 2023
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
3 CommentsShow 1 older commentHide 1 older comment
Stephen23 on 13 Oct 2023
Edited: Stephen23 on 13 Oct 2023
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 on 13 Oct 2023
% 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