Any faster way to calculate sphere distance between two tables of points?

Hi guys,

I have been stuck recently. I'm trying to calculate distance between every point in table A and B, based on their longitudes and latitudes. However, my current solution is quite time-consuming and I think there must be a better way to solve it. With my limited Matlab knowledge, I haven't figured out a new implementation. The tables are shown as below, TableA (103376x6) and TableB (32x3):

I want to get the distance of each point of tableA to all the points in tableB. So the resulting matrix, dist, is a 103376 x 32 matrix. What I have done was to use two for loops, which is so inefficient with huge amount of data. So, how can I do to make it run faster? Appreciate your suggestions and comments.

for i=1:height(A)
    for j=1:height(B)
	dist(i,j) = distance(A.Lat(i),A.Lon(i),B.Lat(j),B.Lon(j))/180*pi*6378.137;
    end  
end

[dpb]

I don't have the TB so can't do anything directly, but... I don't suppose there's anything like pdist2 available? I didn't see anything in the 'See Also' section but the TB doc's are often lacking in the detail the main product gets it seems...

All that comes to mind is that distance is vectorized internally to accept vectors for its inputs. If you built the the complete list of points and passed it instead of the explicit loop you could take advantage of that altho depending upon how it is implemented it may just be the loop anyways...but, that would bypass the overhead of calling the function so many times and that ought to be observable, at least. Whether it would be a sizable fraction of the overall time I suspect is unlikely, but "ya never know!".

Sometimes taking time is the price of trying to deal with large datasets...

The one thing that isn't shown in your above code snippet that could help measurably even with the loop is if you have not preallocated dist so you're doing dynamic reallocation. With array of this size, that would definitely be a noticeable hit.

dist=zeros(height(A),height(B));        % preallocate output array
for i=1:height(A)
  for j=1:height(B)
    dist(i,j) = distance(A.Lat(i),A.Lon(i),B.Lat(j),B.Lon(j));
  end  
end
dist=dist/180*pi*6378.137;

Also you can move the scaling calculation to the end and do it with vectorized operation that might make a tic quicker...

[QIAO WANG]

Thank you for the comment. I just tried in this way but it seems both ways consume same amount of time. I agree with you that it may be the price to deal with large amount of data. Anyway, I'll keep trying to find a better method.

[Jess Lovering]

Can you make 4 matrices that are each 103376 x 32 matrix (one for each of the lat1, lon1, lat2, lon2) and perform the distance calculation once? You could use the repmat command to replicate the lon in TableA 32 times and then repeat that for the lat. Then transpose the lat and lon arrays in TableB to be one row and 32 columns then repmat that 103376 times. Then you would have 4 matrices that are all 103376 x 32 that you could input to the distance command. Not sure what the time difference would be but I imagine it would be much faster.

[QIAO WANG]

That's a brilliant suggestion. Thank you Jessica. I spent some time trying to understand your point but once I understood it, it didn't take me too long to generate these matrices. And your expectation was right. It is now much much faster, just 1 sec.