I've got a question about calculation of great circles in Cartesian coordinates:
Is it possible to compute points along a great circle in UTM coordinates, i.e. given two points P1(easting1, northing1) and P2 (easting2, northing2) is there a formula to get waypoints along the geodesic connecting the two (without converting to geographical coordinates and back)?
Any help is appreciated! Thanks!
The UTM coordinate system, is a kind of Transverse Mercator projection separated in longitude bands and restricted in latitud extent such that the distortions associated with the projection remain small. Also, the UTM coordinate system is conformal projection. Therefore, it preserves the angles. That means that within the UTM zones, a straight line in UTM projection is a good approximation great circle over the relatively short distances that could be found in a UTM zone. If you don't need great accuracy, you can just trace the straight line between P1 and P2.
The following plot shows the maximum deviation that a straight line have from a great circle. In this case P1 the red dot at the Equator and Longitude 3°E, and the maximum deviation is calculated for all points in the displayed box of 6°x6°.
While the differences remain small within distances of few hundred kilometers, bigger differences can be found over larger distances, specially if P1 and P2 are in different UTM zones.
If your application requires larger distances, you shouldn't use UTM projection. Alternatively, if you your application require greater accuracy over distances suitable for UTM projection, you can calculate the great circle trajectory, but I don't know any direct formulas.
What I have used, which I know is not what you are asking for, are some Matlab functions to convert to and from UTM coordinates and another that computes great circle tracks in geographic coordinates:
function [lat, lon] = greatCircle(P1,P2,N)
%GREATCIRCLE compute points along the great circle arc connecting points P1
%and P2
% P1, P2 format: [lat,lon]
% N: Number of points along the great circle arc
lat1=P1(1);
lon1=P1(2);
lat2=P2(1);
lon2=P2(2);
lon12=lon2-lon1;
a1=atan2d(cosd(lat2).*sind(lon12),(cosd(lat1)*sind(lat2)-sind(lat1)*cosd(lat1).*cosd(lon12)));
%a2=atan2d(cosd(lat1).*sind(lon12),(-cosd(lat2)*sind(lat1)+sind(lat2)*cosd(lat1).*cosd(lon12)));
a0=atand(sind(a1).*cosd(lat1)./sqrt((cosd(a1).^2)+(sind(a1).^2).*(sind(lat1).^2)));
sigma01=atan2d(tand(lat1),cosd(a1));
sigma12=atan2d(sqrt((cosd(lat1)*sind(lat2)-sind(lat1)*cosd(lat2)*cosd(lon12))^2+(cosd(lat2)*sind(lon12))^2),(sind(lat1)*sind(lat2)+cosd(lat1)*cosd(lat2)*cosd(lon12)));
sigma02=sigma01+sigma12;
lon01=atan2d(sind(a0)*sind(sigma01),cosd(sigma01));
lon0=lon1-lon01;
sigma=linspace(sigma01,sigma02,N);
lat=atan2d(cosd(a0)*sind(sigma),sqrt(cosd(sigma).^2+sind(a0)^2*sind(sigma).^2));
lon=wrapTo360(atan2d(sind(a0)*sind(sigma),cosd(sigma))+lon0);
end
And for the conversions:
function [x,y,utmzone] = deg2utm(Lat,Lon,Huso)
% -------------------------------------------------------------------------
% [x,y,utmzone] = deg2utm(Lat,Lon)
%
% Description: Function to convert lat/lon vectors into UTM coordinates (WGS84).
% Some code has been extracted from UTM.m function by Gabriel Ruiz Martinez.
%
% Inputs:
% Lat: Latitude vector. Degrees. +ddd.ddddd WGS84
% Lon: Longitude vector. Degrees. +ddd.ddddd WGS84
%
% Outputs:
% x, y , utmzone. See example
%
% Example 1:
% Lat=[40.3154333; 46.283900; 37.577833; 28.645650; 38.855550; 25.061783];
% Lon=[-3.4857166; 7.8012333; -119.95525; -17.759533; -94.7990166; 121.640266];
% [x,y,utmzone] = deg2utm(Lat,Lon);
% fprintf('%7.0f ',x)
% 458731 407653 239027 230253 343898 362850
% fprintf('%7.0f ',y)
% 4462881 5126290 4163083 3171843 4302285 2772478
% utmzone =
% 30 T
% 32 T
% 11 S
% 28 R
% 15 S
% 51 R
%
% Example 2: If you have Lat/Lon coordinates in Degrees, Minutes and Seconds
% LatDMS=[40 18 55.56; 46 17 2.04];
% LonDMS=[-3 29 8.58; 7 48 4.44];
% Lat=dms2deg(mat2dms(LatDMS)); %convert into degrees
% Lon=dms2deg(mat2dms(LonDMS)); %convert into degrees
% [x,y,utmzone] = deg2utm(Lat,Lon)
%
% Author:
% Rafael Palacios
% Universidad Pontificia Comillas
% Madrid, Spain
% Version: Apr/06, Jun/06, Aug/06, Aug/06
% Aug/06: fixed a problem (found by Rodolphe Dewarrat) related to southern
% hemisphere coordinates.
% Aug/06: corrected m-Lint warnings
%-------------------------------------------------------------------------
% Argument checking
%
error(nargchk(2, 3, nargin)); %2 arguments required
CHuso=true;
if nargin==3
CHuso=false;
end
n1=length(Lat);
n2=length(Lon);
if (n1~=n2)
error('Lat and Lon vectors should have the same length');
end
% Memory pre-allocation
%
x=zeros(n1,1);
y=zeros(n1,1);
utmzone(n1,:)='60 X';
% Main Loop
%
for i=1:n1
la=Lat(i);
lo=Lon(i);
sa = 6378137.000000 ; sb = 6356752.314245;
%e = ( ( ( sa ^ 2 ) - ( sb ^ 2 ) ) ^ 0.5 ) / sa;
e2 = ( ( ( sa ^ 2 ) - ( sb ^ 2 ) ) ^ 0.5 ) / sb;
e2cuadrada = e2 ^ 2;
c = ( sa ^ 2 ) / sb;
%alpha = ( sa - sb ) / sa; %f
%ablandamiento = 1 / alpha; % 1/f
lat = la * ( pi / 180 );
lon = lo * ( pi / 180 );
if CHuso
Huso = fix( ( lo / 6 ) + 31);
end
S = ( ( Huso * 6 ) - 183 );
deltaS = lon - ( S * ( pi / 180 ) );
if (la<-72), Letra='C';
elseif (la<-64), Letra='D';
elseif (la<-56), Letra='E';
elseif (la<-48), Letra='F';
elseif (la<-40), Letra='G';
elseif (la<-32), Letra='H';
elseif (la<-24), Letra='J';
elseif (la<-16), Letra='K';
elseif (la<-8), Letra='L';
elseif (la<0), Letra='M';
elseif (la<8), Letra='N';
elseif (la<16), Letra='P';
elseif (la<24), Letra='Q';
elseif (la<32), Letra='R';
elseif (la<40), Letra='S';
elseif (la<48), Letra='T';
elseif (la<56), Letra='U';
elseif (la<64), Letra='V';
elseif (la<72), Letra='W';
else Letra='X';
end
a = cos(lat) * sin(deltaS);
epsilon = 0.5 * log( ( 1 + a) / ( 1 - a ) );
nu = atan( tan(lat) / cos(deltaS) ) - lat;
v = ( c / ( ( 1 + ( e2cuadrada * ( cos(lat) ) ^ 2 ) ) ) ^ 0.5 ) * 0.9996;
ta = ( e2cuadrada / 2 ) * epsilon ^ 2 * ( cos(lat) ) ^ 2;
a1 = sin( 2 * lat );
a2 = a1 * ( cos(lat) ) ^ 2;
j2 = lat + ( a1 / 2 );
j4 = ( ( 3 * j2 ) + a2 ) / 4;
j6 = ( ( 5 * j4 ) + ( a2 * ( cos(lat) ) ^ 2) ) / 3;
alfa = ( 3 / 4 ) * e2cuadrada;
beta = ( 5 / 3 ) * alfa ^ 2;
gama = ( 35 / 27 ) * alfa ^ 3;
Bm = 0.9996 * c * ( lat - alfa * j2 + beta * j4 - gama * j6 );
xx = epsilon * v * ( 1 + ( ta / 3 ) ) + 500000;
yy = nu * v * ( 1 + ta ) + Bm;
if (yy<0)
yy=9999999+yy;
end
x(i)=xx;
y(i)=yy;
utmzone(i,:)=sprintf('%02d %c',Huso,Letra);
end
and
function [Lat,Lon] = utm2deg(xx,yy,utmzone)
% -------------------------------------------------------------------------
% [Lat,Lon] = utm2deg(x,y,utmzone)
%
% Description: Function to convert vectors of UTM coordinates into Lat/Lon vectors (WGS84).
% Some code has been extracted from UTMIP.m function by Gabriel Ruiz Martinez.
%
% Inputs:
% x, y , utmzone.
%
% Outputs:
% Lat: Latitude vector. Degrees. +ddd.ddddd WGS84
% Lon: Longitude vector. Degrees. +ddd.ddddd WGS84
%
% Example 1:
% x=[ 458731; 407653; 239027; 230253; 343898; 362850];
% y=[4462881; 5126290; 4163083; 3171843; 4302285; 2772478];
% utmzone=['30 T'; '32 T'; '11 S'; '28 R'; '15 S'; '51 R'];
% [Lat, Lon]=utm2deg(x,y,utmzone);
% fprintf('%11.6f ',lat)
% 40.315430 46.283902 37.577834 28.645647 38.855552 25.061780
% fprintf('%11.6f ',lon)
% -3.485713 7.801235 -119.955246 -17.759537 -94.799019 121.640266
%
% Example 2: If you need Lat/Lon coordinates in Degrees, Minutes and Seconds
% [Lat, Lon]=utm2deg(x,y,utmzone);
% LatDMS=dms2mat(deg2dms(Lat))
%LatDMS =
% 40.00 18.00 55.55
% 46.00 17.00 2.01
% 37.00 34.00 40.17
% 28.00 38.00 44.33
% 38.00 51.00 19.96
% 25.00 3.00 42.41
% LonDMS=dms2mat(deg2dms(Lon))
%LonDMS =
% -3.00 29.00 8.61
% 7.00 48.00 4.40
% -119.00 57.00 18.93
% -17.00 45.00 34.33
% -94.00 47.00 56.47
% 121.00 38.00 24.96
%
% Author:
% Rafael Palacios
% Universidad Pontificia Comillas
% Madrid, Spain
% Version: Apr/06, Jun/06, Aug/06
% Aug/06: corrected m-Lint warnings
%-------------------------------------------------------------------------
% Argument checking
%
error(nargchk(3, 3, nargin)); %3 arguments required
n1=length(xx);
n2=length(yy);
n3=size(utmzone,1);
if (n1~=n2 || n1~=n3)
error('x,y and utmzone vectors should have the same number or rows');
end
c=size(utmzone,2);
if (c~=4)
error('utmzone should be a vector of strings like "30 T"');
end
% Memory pre-allocation
%
Lat=zeros(n1,1);
Lon=zeros(n1,1);
% Main Loop
%
for i=1:n1
if (utmzone(i,4)>'X' || utmzone(i,4)<'C')
fprintf('utm2deg: Warning utmzone should be a vector of strings like "30 T", not "30 t"\n');
end
if (utmzone(i,4)>'M')
hemis='N'; % Northern hemisphere
else
hemis='S';
end
x=xx(i);
y=yy(i);
zone=str2double(utmzone(i,1:2));
sa = 6378137.000000 ; sb = 6356752.314245;
% e = ( ( ( sa ^ 2 ) - ( sb ^ 2 ) ) ^ 0.5 ) / sa;
e2 = ( ( ( sa ^ 2 ) - ( sb ^ 2 ) ) ^ 0.5 ) / sb;
e2cuadrada = e2 ^ 2;
c = ( sa ^ 2 ) / sb;
% alpha = ( sa - sb ) / sa; %f
% ablandamiento = 1 / alpha; % 1/f
X = x - 500000;
if hemis == 'S' || hemis == 's'
Y = y - 10000000;
else
Y = y;
end
S = ( ( zone * 6 ) - 183 );
lat = Y / ( 6366197.724 * 0.9996 );
v = ( c / ( ( 1 + ( e2cuadrada * ( cos(lat) ) ^ 2 ) ) ) ^ 0.5 ) * 0.9996;
a = X / v;
a1 = sin( 2 * lat );
a2 = a1 * ( cos(lat) ) ^ 2;
j2 = lat + ( a1 / 2 );
j4 = ( ( 3 * j2 ) + a2 ) / 4;
j6 = ( ( 5 * j4 ) + ( a2 * ( cos(lat) ) ^ 2) ) / 3;
alfa = ( 3 / 4 ) * e2cuadrada;
beta = ( 5 / 3 ) * alfa ^ 2;
gama = ( 35 / 27 ) * alfa ^ 3;
Bm = 0.9996 * c * ( lat - alfa * j2 + beta * j4 - gama * j6 );
b = ( Y - Bm ) / v;
Epsi = ( ( e2cuadrada * a^ 2 ) / 2 ) * ( cos(lat) )^ 2;
Eps = a * ( 1 - ( Epsi / 3 ) );
nab = ( b * ( 1 - Epsi ) ) + lat;
senoheps = ( exp(Eps) - exp(-Eps) ) / 2;
Delt = atan(senoheps / (cos(nab) ) );
TaO = atan(cos(Delt) * tan(nab));
longitude = (Delt *(180 / pi ) ) + S;
latitude = ( lat + ( 1 + e2cuadrada* (cos(lat)^ 2) - ( 3 / 2 ) * e2cuadrada * sin(lat) * cos(lat) * ( TaO - lat ) ) * ( TaO - lat ) ) * ...
(180 / pi);
Lat(i)=latitude;
Lon(i)=longitude;
end
I hope that helps.
