Indexing nearest value to my coordinate

Hello everyone, I have a list "Stations" (59x2) with coordinates (lon, lat).

I also have a .nc file that contains longitude & latitude (both 621x837) and a describing variable "v".

I now want to find the indices of the longitude and latitude in the .nc file, that are closest/nearest to each coordinate in my list "Stations" because they do not always match exactly.

The nc file is unfortunately way too large to attach to this question.

My code looks like this:

%% Data import and information
filename = 'D:\xxx\ccc\vvv\1999.nc';

%% Step 1 - Variables: Load latitude and longitude arrays
lon_data = ncread(filename, 'lon');
lat_data = ncread(filename, 'lat');
t_data = ncread(filename, 'time');

%% Step 2: Define the coordinates of the 59 points of interest
Stations = [14.0355, 61.2536;
    ... % 58 other coordinates that I can add if needed (otherwise it makes the code confusing)
    ];

%% Step 3: Find the indices of the points of interest
lat_indices = zeros(59, 1);
lon_indices = zeros(59, 1);
eps = 0.0005;

for i = 1:59
    [~, lat_idx] = abs(lat_data - Stations(i, 2));
    [~, lon_idx] = abs(lon_data - Stations(i, 1));
    lat_indices(i) = find((lat_idx < eps), 1);
    lon_indices(i) = find((lon_idx < eps), 1);
end

The code, however, returns the following error:

Error using abs

Too many output arguments.

Error in Indexing_lat_lon (line 83)

[~, lat_idx] = abs(lat_data - Stations(i, 2));

Answer 1

[Torsten]

lon_data = rand(621,837);
lat_data = rand(621,837);
stations = rand(59,2);
lon_indices = zeros(59, 1);
lat_indices = zeros(59, 1);
for i = 1:59
    [~,lon_indices(i)] = min(abs(lon_data-stations(i,1)),[],"all","linear");
    [~,lat_indices(i)] = min(abs(lat_data-stations(i,2)),[],"all","linear");
end
[row_lon,col_lon] = ind2sub([621 837],lon_indices)

row_lon = 59×1

434 365 272 327 150 334 11 98 91 546

col_lon = 59×1

150 441 239 534 214 518 83 667 590 484

[row_lat,col_lat] = ind2sub([621 837],lat_indices)

row_lat = 59×1

370 291 378 195 381 9 372 577 473 186

col_lat = 59×1

544 789 218 24 61 390 692 115 250 825

It's strange that you want the min for lon and lat separately. Shouldn't you match one index pair (ilon,ilat) common for both nc-files ?

Shouldn't it be

for i = 1:59
    [~,lonlat_indices(i)] = min((lon_data-stations(i,1)).^2+(lat_data-stations(i,2)),[],"all","linear");
end
[row_lonlat,col_lonlat] = ind2sub([621 837],lonlat_indices)

row_lonlat = 1×59

409 598 21 321 328 421 2 446 230 328 256 266 561 421 162 262 118 12 598 309 309 21 561 468 342 201 162 266 356 21

col_lonlat = 1×59

157 598 680 815 763 567 227 416 219 763 570 240 646 567 348 288 300 163 598 642 642 680 646 716 564 789 74 731 276 680

Answer 2

[Bruno Luong]

It is odd that you want to find the index of nc file data closest to the Station.

It just sounds more natural to look for the opposite, index of the Station closest to the nc data.

Note that the distance here is euclidian on (lon,lat) which is NOT the geodesic distance.

lon_data = rand(10,10);
lat_data = rand(10,10);
Stations = rand(59,2);

DT = delaunayTriangulation(lon_data(:), lat_data(:));
idx = DT.nearestNeighbor(Stations)

idx = 59×1

79 11 22 3 30 68 30 13 60 81

% alternatively if you have the Statistics and Machine Learning Toolbox 
%idx = knnsearch([lon_data(:), lat_data(:)], Stations)

figure
h1 = plot(lon_data, lat_data, '+g');
hold on
h2 = plot(Stations(:,1),Stations(:,2), 'or');
Lon = [Stations(:,1),lon_data(idx)]'; Lon(3,:) = NaN;
Lat = [Stations(:,2),lat_data(idx)]'; Lat(3,:) = NaN;
h3 = plot(Lon, Lat, 'b');
legend([h1(1) h2(1) h3(1)], 'nc data', 'stations', 'nearest nc')
axis equal