Convert Long / lat to UTM (WGS 84)

[Eugene Dongmo]

ood evening good people.
I want to build an app to convert Longitude / Latitude into UTM, then use it to measure surface areas. I have used the algorithm developed by Brenor Brophy. But the result obtained is having a big gap according to the others software. Let say, the gap is too much. I will like to have a pricision in centimeters as the author the code said.
Please just look at it and tell me where is the error
Thanks you for all your work.

{
double Easting;
double Northing;
int Zone;
char Letter;
private Deg2UTM(double Lat,double Lon)
{
Zone= (int) Math.floor(Lon/6+31);

    Easting=0.5*Math.log((1+Math.cos(Lat*Math.PI/180)*Math.sin(Lon*Math.PI/180-(6*Zone-183)*Math.PI/180))/
        (1-Math.cos(Lat*Math.PI/180)*Math.sin(Lon*Math.PI/180-(6*Zone-183)*Math.PI/180)))*0.9996*6399593.62
    /Math.pow((1+Math.pow(0.0820944379, 2)*Math.pow(Math.cos(Lat*Math.PI/180), 2)), 0.5)*
    
    (1+ Math.pow(0.0820944379,2)/2*Math.pow((0.5*Math.log((1+Math.cos(Lat*Math.PI/180)*Math.sin(Lon*Math.PI/180-
        (6*Zone-183)*Math.PI/180))/(1-Math.cos(Lat*Math.PI/180)*Math.sin(Lon*Math.PI/180-
        (6*Zone-183)*Math.PI/180)))),2)*Math.pow(Math.cos(Lat*Math.PI/180),2)/3)

    +500000;
    
    Easting=Math.round(Easting*100)*0.01;

    
    Northing = (Math.atan(Math.tan(Lat*Math.PI/180)/Math.cos((Lon*Math.PI/180-(6*Zone -183)*Math.PI/180)))-
    Lat*Math.PI/180)*0.9996*6399593.625/Math.sqrt(1+0.006739496742*Math.pow(Math.cos(Lat*Math.PI/180),2))*
    (1+0.006739496742/2*Math.pow(0.5*Math.log((1+Math.cos(Lat*Math.PI/180)*Math.sin((Lon*Math.PI/180-(6*Zone -183)*
    Math.PI/180)))/(1-Math.cos(Lat*Math.PI/180)*Math.sin((Lon*Math.PI/180-(6*Zone -183)*Math.PI/180)))),2)*
    Math.pow(Math.cos(Lat*Math.PI/180),2))+0.9996*6399593.625*(Lat*Math.PI/180-0.005054622556*
    (Lat*Math.PI/180+Math.sin(2*Lat*Math.PI/180)/2)+4.258201531e-05*(3*(Lat*Math.PI/180+Math.sin(2*Lat*Math.PI/180)/2)+
    Math.sin(2*Lat*Math.PI/180)*Math.pow(Math.cos(Lat*Math.PI/180),2))/
    4-1.674057895e-07*(5*(3*(Lat*Math.PI/180+Math.sin(2*Lat*Math.PI/180)/2)+Math.sin(2*Lat*Math.PI/180)*
    Math.pow(Math.cos(Lat*Math.PI/180),2))/4+Math.sin(2*Lat*Math.PI/180)*Math.pow(Math.cos(Lat*Math.PI/180),2)*
    Math.pow(Math.cos(Lat*Math.PI/180),2))/3);
    
    if (Letter<'M')
        Northing = Northing + 10000000;
    Northing=Math.round(Northing*100)*0.01;
}

}

Here is the aia file, Please I really need help

[SteveJG]

You did not provide an example of a 'correct' conversion versus the conversion you achieved using the Math blocks.  How is anyone to know what you believe is a correct conversion of latitude/longitude to UTM?  Provide an example of a correct conversion and the AI2 result to help anyone who wants to check results.  Unfortunately, the Math blocks are almost impossible to read when using AI2.  One way to debug them is to use DoIt   and / or to break the equation down into small pieces.   You may need to do that.

1) Did you use the correct constants for  Equatorial Radius and  square of eccentricity.  For WGS 84 you should be using 6378137 and  0.00669438 I believe ( http://2cyr.com/emaps-bg.php?q=source ) .  These blocks in your code might not be correct:

2)You converted a php algorithm to App Inventor 2 Blocks.  When working with geometric functions, you may have to convert radians to degrees or degrees to radians to get the algorithm to work.  I am not sure the trig functions in php work the same way as they do in App Inventor 2.  You could experiment.  For instance, the TAN function in App Inventor converts degrees to radians 

3)  You might need to adjust the 'estimate' of the Earth's ellipsoid ... you use 6378137 .  Other values are commonly used       https://en.wikipedia.org/wiki/Earth_ellipsoid   .  The value might be different in the programs you compare your App Inventor results with however 6378137 should be the value to use

4) You are taking your coordinates from the LocationSensor.  A GPS receiver only provides latitude and longitude to FIVE decimal places. The latitude and longitude are provided by a NEMA sentence transmitted by the GPS satellites and it only reports five decimal places  That means one can only resolve  0.00001 of a degree of latitude  which is about 3 feet (or about 94 cm) if your positional information comes from the GPS.   You might try entering coordinate data that contains more than five decimal places accuracy instead of using the values provided by the LocationSensor and see what happens.

5) You could use a javascript UTM/lat/lon conversion algorithm within App Inventor to do the calculation and avoid using the Math blocks

https://github.com/TimothyGu/utm

or 

https://www.movable-type.co.uk/scripts/latlong-utm-mgrs.html 

You would do something like with the javascript file  https://groups.google.com/forum/#!msg/mitappinventortest/0oVTthDFChY/GmIBOEH_EQAJ;context-place=searchin/mitappinventortest/javascript%7Csort:date

or   https://appinventor.mit.edu/explore/ai2/webview-javascript.html 

Also, you might read a previous discussion about using App Inventor to convert coordinate systems  https://groups.google.com/forum/#!msg/mitappinventortest/K-bbAZmFbSo/RPudd1gmBgAJ;context-place=categories/mitappinventortest  

Hope something above helps you.

Regards,

Steve

[ABG]

The Right Click block option to switch to External View

can be used to show the nesting structure differently,

to help untangle long formulas.

ABG

[Eugene Dongmo]

Thank Steeve for all your work.

This is an example :

Decimal degree : Long = 9.25026 Lat = 4.15636

My conversion UTM = X = 501457 Y = 460127 32N

Other applications: X= 527775 Y = 459415 32N

The second conversion is good according to my GIS / Mapping knowledge

[SteveJG]

OK..that is helpful.

Here is what I see

I get slightly different values of x and y using 9.25026 and 4.15636 longitude/latitude in my Android 8.1 Tablet.

If I use an online UTM converter, I get this:

which is 52775 and 458310 .  These values are  not what you see using other applications. X= 527775 Y = 459415 

Either 

1)  your algorithm has issues and is not robust- how do you know  Brenor Brophy's algorithm calculates UTM correctly in its original version (not your translation to Blocks)?

2)  you made a mistake copying and converting to Blocks  ... use the DoIt debugger.    You may have to do the calculation manually and compare it with what the Math blocks show.step by step

3) App Inventor's Math Blocks treat Math procedures differently or round  differently or use fewer significant figures in the calculation or the system Android math routines the Math blocks use are slightly different in some Android distributions or this is an issue with the calculations possible using different cpu  or......

4) You might be able to do the calculations using the free Math parser https://puravidaapps.com/math.php  instead of using Blocks.

You  have to debug this your self Eugene unless someone decides to try to solve the puzzle.  When you find out what is going wrong, please let us know.  The discrepancy of calculate UTMs could be an issue with how the Math blocks work or an input error on your translation to Blocks.

--Steve

[Eugene Dongmo]

You're great Mr Steve. Thanks for all your works.

With this algoritm that you gave me the link (https://github.com/TimothyGu/utm), I am able to get X value (approximate good value) but Y value remains a problem, I can"t get the right value.

[Eugene Dongmo]

I found ne of the problems : the blocks Sin, Cos, Tan are not returning the good values. They are returning values in DEG instead of RAD. But my difficulty here is that, I don't know to get value in RAD.

The values were first converted into RAD (°xPI/180)

But when I use the blocks SIN, COS or TAN, the values return in DEG.

Please help me to solve this

[ABG]

Right, AI2 trig functions eat degrees, not radians.

So don't feed them radians.

See attached, where I use separate global variables with degrees values.

ABG

[Eugene Dongmo]

Thanks ABG, but this is my problem, I followed your instructions, but still having different values.

Is there a possibility to return values using RAD mode like a calculator? The calculator is proposing a good answer.

Thanks for your help

[ABG]

Data moves through blocks from right to left.

There is no MODE concept in AI2, just some blocks

that can do the multiplication/division of pi/180 (whatever)

you can use to feed trig functions.

sin and cos MUST be fed degrees.

Once sin and cos emit a value, that value is no longer an angle, so it

is an abomination upon mathematics to apply radian conversion to it.

If you want to show radians, make a copy to convert for display.

You must use the Companion's Do It facility in the blocks editor

to trace values from right to left.

Do it until your blocks look like a pin cushion.

ABG

[Eugene Dongmo]

This another demonstration

[SteveJG]

App Inventor requires the value in the bracket () of a trig function to be presented in degrees, for example  sin(degrees) or cos(degrees)

When you multiply  degrees by 3.14/180 you convert the degrees to radians  .........1deg = pi/180 = 0.0174 radians  You are 

 using a radian in the bracket instead of degrees.   The result of  sin(radian) is not what you want but that is what you are calculating.

Also you are calculating a  ln ... the algorithm you are converting from uses a log function  (similar but different).

[Eugene Dongmo]

I removed all pi/180 and my X value is correct, but the Y value still to be corrected, 

I really appreciated the way you took your time to help me. Many thanks for all your contributions.

Please If someone has an aia on how to correct the Y value, please share with me so I can make the calculations complete

Thanks one more

[ABG]

https://community.thunkable.com/t/convert-longitude-latitude-to-utm-wgs-84/92616/13

[Eugene Dongmo]

Thanks you for all your help

In the course of my research, I saw this algorithm, it is the one I am using now. 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

[SteveJG]

Hi Eugene... does the new Java algorithm work when converted App Inventor's Math Blocks?   Would you like to share the aia and Block images with us here or are you still experimenting?

-- Steve

[Eugene Dongmo]

I am still experimenting. Seeing the way the algo is, the calculation will be easier

[Peter Mathijssen]

Eugene had also a lengthy discussion here with the same question.

https://community.thunkable.com/t/convert-longitude-latitude-to-utm-wgs-84/92616