Convert Longitude Latitude to UTM (WGS 84)


#22

Yes, but that is not the point.
If the equation is making some kind of mathematical correction fix that is a polynomial of some sort, that relies on fixed coefficients that ALREADY take into account the fact that the result would need to be converted, then those corrections factors would also need to be converted. Basically, they might have the pi/180 hidden.

For instance, you have that coefficient 6399593.625 in the original equation. What is the meaning of that number? Is it possible that it could be 366669705 already multiplied by pi/180? (I am not saying that it is, only that it could be).
What about 4.258201521e-5? Could it apply only in the context of working in RAD, and that it would need to be adjusted to be 0.00243976976 instead?

What you need to do is to process the equations in parallel, on a calculator in RAD mode (where you are confident that the result is correct, including the multiplications by pi/180) and in DEG (where you omit the pi/180). Whenever you get to an intermediate result that WILL NOT be the argument of a trigonometric function (sine or cosine), then the values should be equal. Whenever you have a result that will be the argument of a sine or cosine function, they have to be off by a factor of 57.2957795 (180/pi) between the RAD and the DEG.
And if the value is the argument of a reverse trigonometric function (ATAN) then they have to be euqal, but the RESULT of the ATAN function has to be scaled by 180/pi.

That is why you should have made yourself a favor and turn those huge equations in a bunch of smaller ones, so that you could track and compare the intermediate results.


#23

Thanks you very much for all your help CBGV.
You helped a lot. I learned also to much in mathematics formulas.
May God Bless you and the hole Community.

In the course of my research, I saw this algorithm, it is simple and detailed. They were so many difficulties to implement the other java code.

function [x,y,utmzone] = deg2utm(Lat,Lon)

% -------------------------------------------------------------------------

% [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, 2, nargin)); %2 arguments required

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 );

Huso = fix( ( lo / 6 ) + 31);

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 );

Easting = epsilon * v * ( 1 + ( ta / 3 ) ) + 500000;

Northing = 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