3
$\begingroup$

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!

$\endgroup$
  • 1
    $\begingroup$ Beware that for most applications where a great circle matters, you will be crossing more than one UTM zone. $\endgroup$ – Semidiurnal Simon Apr 10 at 4:38
3
$\begingroup$

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°.

enter image description here

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.

$\endgroup$
  • $\begingroup$ Are you sure great circles plot as lines in UTM? I know that's true for rhumb lines but surely not for any great circle other than a meridian? Wikipedia says " straight lines on the map (rhumb lines), other than the meridians or the equator, do not correspond to great circles" $\endgroup$ – bweise Apr 11 at 8:28
  • 1
    $\begingroup$ @bweise It depends. Explain your application. If you are worried with centimeter accuracy, then straight lines in UTM are not great circles. But otherwise they are a very good approximation WITHIN a UTM zone. Mercator is NOT the same than UTM, because UTM zones are thin enough to keep distortions small, they are just 6° wide, this is about 660 km at most. $\endgroup$ – Camilo Rada Apr 11 at 13:26
  • 1
    $\begingroup$ @bweise See the edits. I added info about the errors if you use just a straight line, and a workaround to compute great circle tracks but converting to geographic. If there is a direct formula, you might have better luck asking at Geographic Information System SE. $\endgroup$ – Camilo Rada Apr 11 at 18:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.