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

2020 8 19
1 个回答
回答已采纳
更新时间：2020 8 21
7 次查看（30 天）

0 个投票

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?

2 个评论
显示 无 隐藏 无

Bruno Luong
Bruno Luong 2020-8-19
编辑：Bruno Luong 2020-8-19
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
Sai Prasanth 2020-8-20
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/

请先登录，再进行评论。

请先登录，再回答此问题。

 采纳的回答

Rik
Rik 2020-8-19

1 个投票

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

6 个评论
显示 4更早的评论 隐藏 4更早的评论

Sai Prasanth
Sai Prasanth 2020-8-20
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
Sai Prasanth 2020-8-20
I am also researching on whether I can use rangesearch and avoid all of this. Please let me know what you (all) think.
Rik
Rik 2020-8-20
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
Sai Prasanth 2020-8-21
编辑：Sai Prasanth 2020-8-21
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
Rik 2020-8-21
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
Sai Prasanth 2020-8-21
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

请先登录，再进行评论。

更多回答（0 个）

类别

在 帮助中心 和 File Exchange 中查找有关 Matrix Indexing 的更多信息

产品

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by