Compute the minimum distance between each element in a vector with another larger array

I have an array say, MW that contains the latitudes and longitudes of ~1300 pixels. I also have a reference text file, say IR that contains an enormously large number of pixels with latitudes and longitudes and the pixels here are associated with unique IDs. My goal is to attribute each MW pixel to a unique ID found in the IR text file if possible. I have a function called distance_latlon (function d = distance_latlon(lat1,lat2,lon1,lon2)) that computes the distances between any two points given their lat and lon. Essentially, I would like to compute for each MW pixel, the distance to all of the IR pixels and find the IR pixel within the closest distance and attribute the ID. If the MW pixel does not have a single IR pixel within a certain distance threshold, then I can simply assign a 0.

How do I accomplish this task without the extremely expensive way that is, going through each MW pixel in a for loop, running the distance function between each MW pixel and all the IR pixels in the text file in another long for loop and finding the minimum distance? What is the least expensive way of doing this?

[Bruno Luong]

What is the size of IR (something more precise than "enormously large number")? If you provide the size of MW, it would be nice if you provide the size of IR as well.

What returns distance_latlon(lat1,lat2,lon1,lon2)? Is it the geodesic distance on sphere?

[Sai Prasanth]

It was initially of the size of 100s of millions (but I was able to subset it and bring it down to 76,000). And yes, it is the distance between two points given their lat/lon on a sphere (Heaversine formula) described here: https://www.geeksforgeeks.org/program-distance-two-points-earth/

Answer 1

[Rik]

I don't think there is a way around doing some calculation for each pair.

Is there a way to vectorize distance_latlon? If not, I would suggest a pre-selection, so you only need to do the expensive calculation a few times:

N_small=5;
[lat2,lon2]=MyFun(IR);
for n_MW=1:1300
    [lat1,lon1]=MyFun(MW(n_MW));
    quickdist=hypot(lat1-lat2,lon1-lon2);
    [~,idx1]=sort(quickdist);
    lat2_small=lat2(idx(1:N_small));
    lon2_small=lon2(idx(1:N_small));
    for n_small=1:N_small
        exactdist(n_small)=distance_latlon(lat1,lat2_small,lon1,lon1_small);
    end
    [min_exactdist,idx2]=min(exactdist);
    rowID=idx1(idx2);
end

[Sai Prasanth]

Thank you. Can you please elaborate on what you mean by vectorizing distance_latlon? Do you mean, can it take the inputs as lat-lon vectors? Sure enough, but lat_1 and lon_1 (derived from MW) will be one value and lat_2 and lon_2 (derived from IR) will be a huge array. Perhaps, I can conform the dimensions of lat_1 and lon_2 to the size of lat_2 and lon_2 and conduct the operation? The below function is the version that does it for one pair of lat, lon.

function d = distance_latlon(lat1,lat2,lon1,lon2)

% This function calculates the distance between two points presented in
% longitudes and latitudes (x1,y1 or lon1,lat1) and (x2,y2 or lon2,lat2)
% Here, the input must be in degrees and decimals. 

d = 0; 
radEarth = 6371; % In km

if(lon1-lon2 < -180 || lon1-lon2 > 180)
disp(['Add or subtract 360, your lon1-lon2 is: ',num2str(lon1-lon2)])
return
end

term1 = cos(pi/360 .* (lat1 + lat2)) .* pi/180 .* (lon1 - lon2); 
term2 = pi/180 * (lat1 - lat2); 

d = radEarth .* sqrt(term1.^2 + term2.^2); 

disp(['The distance in km is: ',num2str(d)])

end

[Sai Prasanth]

I am also researching on whether I can use rangesearch and avoid all of this. Please let me know what you (all) think.

[Rik]

As long as one pair is scalar, the function you showed will allow you to enter a vector of coordinates. That is what I meant by vectorizing the code.

[Sai Prasanth]

I understand! In this case, do I just run the MW loop and each time, feed the vector of distances each time and skip all the Quickdist etc.? The code might look something like the described one below. Is this the best way to do it?

[lat2,lon2]=MyFun(IR);

for n_MW=1:1300
  [lat1,lon1]=MyFun(MW(n_MW));   
   
   conformed_lat1 = conform_dims(size(lat2),lat1,2) % lat2 and lon2 are vectors
   conformed_lon1 = conform_dims(size(lon2),lon1,2)
   
    distfromMW  = distance_latlon(conformed_lat1,lat2,conformed_lon1,lon2);
    % distfrom MW, the output is also a vector now
    
    [mindist,index_mindist] = min(distfromMW)
    
    if mindist <= distance_threshold
    ID_MW(n_MW) = IR_file_data(index_mindist,end); 
    % The ID that needs to be assigned is in the last column of the IR text file, so I
    % am picking the row corresponding to the minimum distance and the last column. 
    else 
    ID_MW(n_MW) = 0; 
    end % end if
    
 end

I

[Rik]

Sure, the part with quickdist was based on the assumption that running it directly was too slow.

That function does seem to support vector input, so I would suggest you try it to see if it suits your need.

[Sai Prasanth]

After you pointed me towards thinking about vector inputs, I think using the combination of Matlab's distance and rangeselect functions are the quickest way to get this done.

a = -90;
b = 90;
rand_lat1 = (b-a).*rand(1300,1) + a;
rand_lat2 = (b-a).*rand(76000,1) + a;

a = -180;
b = 180;
rand_lon1 = (b-a).*rand(1300,1) + a;
rand_lon2 = (b-a).*rand(76000,1) + a;

MWpixels = [rand_lat1 rand_lon1]; % Rep. the MW pixels
IRpixels = [rand_lat2 rand_lon2]; % Rep. the IR pixels

% Objective: To find the IR pixels that are close to the MW pixel. 
tic
% rangesearch(X,Y,r): Finds all the X points that are within distance r of the Y points
[Idx,D] = rangesearch(IRpixels,MWpixels,1000,'Distance',@distfun); 
toc

%%%%%%%%%% Distance Function %%%%%%%%%%%
function D = distfun(vec_latlon_1,vec_latlon_2)

% The first column represents latitudes and the second column represents
% longitudes for both inputs

rad_Earth = 6371; % In km
D = distance(vec_latlon_1,vec_latlon_2,[rad_Earth 0]); 

end