Вот некоторые фрагменты кода, написанные на Matlab, которые могут быть полезны.Если это полезно, вы должны конвертировать их в Python.Подход грубой силы и не самый элегантный.Тем не менее, я попытался включить приближенные преобразования системы координат, которые учитывают форму Земли как эллипсоид.Все может немного упроститься, если считать Землю сферой.В качестве альтернативы, чтобы повысить точность (хотя это весьма вероятно на незначительную величину), можно локально аппроксимировать поверхность эллипсоида поверхностью сферы (сферы, которая лучше всего приближает эллипсоид в данной точке) и использовать сферический вместо евклидовагеометрия.
Возможны некоторые опечатки или ошибки, но, возможно, вы сможете получить представление о координатах, преобразованиях и методе.
С помощью следующих двух функций вы можете преобразовать в:
геодезические (то есть долгота широта) координаты вблизи точки long_lat0 = [long0, lat0] с евклидовыми координатами, которые являются линейным приближением первого порядка фактических, истинных геодезических координат на эллипсоиде Земли WGS84
и наоборот, вы можете преобразовать обратно из евклидовых координат в геодезические длинные широты
long_lat0 = [long0, lat0]; % a point from dataset 2
long_lat % the n x 2 matrix of points from dataset 1 (or a chunk of it)
%center of approximate Euclidean coordinate system is point long_lat0
% with long_lat coordinates and the scaling coefficient
% a of longitude and b of latitude,
% which equalizes longitude and latitude distance at point long_lat0, is
function [x, a, b] = convert_to_local_Eucl(long_lat, long_lat0)
% long_lat0 = [long_0, lat_0] is the origin of the local coordinate system
% long_lat = [long_1, lat_1;
% long_2, lat_2;
% ............
% long_n, lat_n] is an n x 2 array of points in lat and long coordinates
% on the Earth's ellipsoid
% x = [x_1, y_1;
% x_2, y_2;
% ..........
% x_n, y_n]
% is the n x 2 matrix of Euclidean coordinates with origin the point long_lat0
% a is a number, correction factor of longitude coordinate
% b is a number, correction factor of latitude
R = 6378137.0 %in meters;
e_2 = ( R^2 - (6356752.314245)^2 ) / R^2;
a = R * (1-e_2) * cosd(long_lat0(2)) / (1 - e_2*sind(long_lat0(2))^(1/2)); % dlong
b = R * (1-e_2) / (1 - e_2*sind(long_lat0(2))^(3/2); %dlat
% a and b are correcting/rescaling coefficients
% that correct the longitude-latitude coordinates of all points
% near point long_lat0 in geodetic coordinates of WGS84.
x = long_lat .- long_lat0; % subtract the long_lat0 from the coordinates of every point
% from the list long_lat, i.e. for each j = 1...n
% x(j, 1) = long_lat(j, 1) - long_lat0(1);
% x(j, 2) = long_lat(j, 2) - long_lat0(2);
x = [ a * x(:,1), b * x(:, 2)];
% multiply the first column of coordinates by the scaling factor a and
% multiply the second column of coordinates by the scaling factor b
% these coordinates are first order linear Euclidean approximation
% of the real geodetic coordinates of WGS84.
% Near the point long_lat0
% the error is negligible, especially within a couple of kilometers.
% The farther you go from that point, the error slowly increases,
% but then it doesn't matter since such points are not the closest anyway.
end
function long_lat = convert_to_long_lat(x, long_lat0, a, b)
% from Euclidean coordinates x = [x(1), x(2)] of a point near long_lat0 go back to
% long_lat = [long, lat] coordinates of that points. a and b are the scaling
% coefficients at point long_lat0
long_lat = [long_lat0(1) + x(1)/a, long_lat0(2) + x(2)/b];
end
Для каждой точки long_lat0 = [long0, lat0] из набора данных 2 начните с преобразования изгеодезический long-lat для аппроксимации евклидовых координат в long_lat0 весь (или часть) long_lat список набора данных 1 секунда и третий столбец:
x = convert_2_local_Eucl(long_lat, long_lat0);
Затем вычислите величины (то есть длины) всех 2Dвекторы-строки x(j,:) = [x(j,1), x(j,1)] из набора данных x
magnitudes = norm(x); %you have to either find this function or write one yourself
После этого найдите индекс и минимум элемента из x:
[j, min] = min(magnitudes);
Тогдадля двух пар: x1 = x(j,:) and x2 = x(j+1,:) и x1 = x(j,:) and x2 = x(j-1,:) используйте следующую функцию для вычисления ближайшей точки:
function [dist, long_lat] = dist_point_to_reference(x1, x2, long_lat0, a, b)
% calculates the shortest distance dist from the point long_lat0
% to the closest point on the segment between x1 and x2
% and then obtain the long_lat coordinates of this closest point
dist = dot(x1, x1) * dot(x2 - x1, x2 - x1) - dot(x1, x2 - x1)^2 ; % dot is dot product
dist = sqrt( dist / ( dot(x2 - x1, x2 - x1)^2) );
% dist is the distance from the point at the origin [0, 0]
% to the straight Euclidean interval between
% the points x1 = [x1(1), x1(2)] and x2 = [x2(1), x2(2)]
if dot(x1, x2 - x1) > 0 % if the height of the triangle is outside, on the side of x1
dist = sqrt( dot(x1, x1) );
long_lat = x1;
elseif dot(x2, x1 - x2) > 0 % if the height of the triangle is outside, on the side of x2
dist = sqrt( dot(x2, x2) );
long_lat = x1;
else
long_lat(1) = - x2(2) + x1(2);
long_lat(2) = x2(1) - x1(1);
long_lat = long_lat / sqrt(dot(long_lat, long_lat));
long_lat = - dot(x1, long_lat) * long_lat; % despite the name, these are Eucldean coordinates
end
long_lat = convert_to_long_lat(long_lat, a, b); % finally, geodetic coordinates
end