Convert Longitude Latitude to UTM (WGS 84)


#1

Good 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 gape according to the other softwares. Let say, the gap is too much.
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

GPS_UTM_2.aia (15.7 KB)


#3

Just looking at it, my first reaction is “(*^%! that is one big nasty equation”.

How accurate is the input value? By this, I mean if you enter a Longitude and Latitude, what is the smallest difference delta that would return a different result when entering Longitude+delta, and when entering Latitude+delta?

You are getting the latitude and longitude value from the LocationSensor, which is presumed to be accurate enough (and if it is not, that is pretty much the end of it, you are stuck with it). But you are saving them in text fields, as opposed to save them in numerical variables. And that may be the source of problem: the text field will drop decimals to fit in the window, and it a poor representation of the real numeric value.
What I suggest is that you save latitude and longitude in global variables, and use THOSE to populate the text field (for user to see), and use THOSE instead of the text in the calculation (for the computer to get super good data, not the “for human formatted” version).
In a similar fashion, you are saving “ZLab”, used for additional calculations, in a text string (label) that is not even displayed (its Visible property is unchecked in the designer). Again, this value should be kept as a number, in a local variable.

Also, for your sake, mine and that of many trillions of electrons, consider simplifying the equations.
You are (re)calculation the segment “cos(Lat*pi/180)” 6 times in the evaluation of “Xlabel” and the segment “cos(Lat * D2R)” 8 times in the second (and there are several other snippets that are multi duplicates).
Those values will not change in the middle of the equation, how about computing them once, assigning them to a properly named local variable, and using those instead in the large equation?
You may even find that a few of those are combined in the same way in a few places, and could also be pre-calculated with the result used at several places.

Last comment: I do not understand your intention with this “Math.round(Easting * 100) * 0.01”.
The way it is implemented is "Math.round { (Easting * 100) * 0.01 } ", so the double multiplication is doing absolutely nothing to the result.


#4

Thanks you Sir for this first and detail answer, yeah the last comment does not help, I just wanted to check if the problem is not there.
At the level of keeping the Long / Lat in variable, I understand, thanks, I will try and get back to you.

Thanks you ones more


#5

I just did the adjustments but the result is the same.
Here is the adjusted one, please can you go trhough the formula and check if the implementation which is the problem? If the various are put in the right way.

GPS_UTM_2.aia (11.9 KB)


#6

I had to do some research as I never heard of Universal Transverse Mercator coordinate before.

So, imagine how little I know of Brenor Brophy’s work…

Beyond the implementation of the math (by the way, you could have used local variables instead of global ones, not that it matters much due to the comparatively small size of the app), what it is that you call “a big gap”?

I assume that you want to use your application to go to the 4 corners of a field, obtain an X and Y coordinates from those, so that you could derive how large the field is. Is this correct?

In this context, what are you calling a “big gap”?

I just started a small GPS monitoring app and am looking at it right now. The accuracy that it reports, with 18 satellites in view and 8 in use, is usually stated to be +/- 10 m. Over a 5 minute period, I have seen it go anywhere between 9 and 24 m, as satellites move in and out.

With this in mind, a 20 m by 20 m backyard would reportedly be measured at anything between 100 square metre right up to 900 m^2, depending on which side the error goes. And since that value fluctuates all the time, you can’t even rely on this being predictable.

(Actually, I just saw it jump right up – degrade! – to 49 m accuracy. It is now coming back down to 12, but those spikes can really mess up any measurement.)

If this is the kind of gap you have in mind, I am afraid that it the limit of what the GPS can do.
The location sensor can output its accuracy, in “m”. Perhaps you should output that on the primary display of your app so that you would know how much trust you can have on the X and Y location values at every moment.

If you intend it to measure the area of a farmer’s field – to estimate how much wheat can be grown – of a large vacant plot around a projected large shopping center – to estimate how much asphalt will be required to pave the parking lot; then the location sensor will probably give a result within a few % point of what an actual measure done by a surveyor would return.
If you intend it to estimate how large the front lawn of a house is so that a kid mowing grass for neighbors would know how much to charge for the service, then the error is likely to be larger than the answer itself. I am afraid that you would simply have reached the limit of what that tool can measure.


#7

Thank you ones more CBVG for all.

  • UTM is a coordinate system like Degree Minute Second or Decimal Degree (proposed by location sensor). We can transform one system to get another one. However, UTM is showing the coordinates in meters while the others formats show the coordinates in degree. With UTM, it is easy to have an accurate surface surface area than using degree to calculate.

  • What I call ‘Big gap’ here is that, comparing with another apps and even my GIS software, the value is not the same, and when I measure the difference, I some he gap of 27 km. So if you collect a point at a particular zone, the point will be located at 27 km away from where you’re now. That is incorrect, I have done some checking and realise that the formular was correct, but I don’t know what is happening with my blocks (instead of calculating using the RAD mode, the block are using the DEG mode (like on the calculator) but I have done the conversion at the begining) .
    So the accuracy here does not have a problem, is the maths blocks that are not calculating well.

That’s my problem.
Thanks


#8

OK, I think I see one problem.
In the calculation, you are evaluating “SinLong” and “CosLong” and both of those depend on the value of “Zlab” , which you evaluate after it being used.

The thunkable math documentation states that the trigonometric functions work in degree. And that the latitude and longitude values returned by the location sensor are also in degree. But all those multiplications by Pi/180, that is to convert degree to radians.

I suggest to remove all conversions, and perhaps the calculation will be accurate.


#9

I have found a problem, the Sin, Cos and Tan blocks are not returning the values in RAD, but they return the value in Degree. That’s maybe the problem of the gap, but I don’t know how to fix this problem. Please Help. The Long and Lat were already been converted in RAD, but whe I apply the Sin, Cos or Tan blocks return the value wrongly


#10

You solve the problem by NOT using LatRad and LongRad, but by using LatDeg and LongDeg.
Do NOT convert degrees to radian, the sin, cos and tan functions expect their input to be in degree.


#11

Thanks CBVG, 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


#12

OK, you seem to be quite confused with trigonometry.

First of all, do you know how much 31.24719 radian is?
It is 1790 degree. That is nearly 5 times full circle. What mathematical value is there to compute values that are over a full circle? If you need to know the distance between London and Paris, do you include a path that pass through New York, then Beijing, then around once more to Buenos Aires, then Melbourne and (why not?) one more time through Madrid, Los Angeles and Moscow?

You put your calculator to RAD mode. That means it will compute the sine of 31.24719 RAD, and give a result that is a pure, non dimensional number. You do NOT convert the result of a sine result to radian, it makes no sense. The result of the sine or cosine function is NOT AN ANGLE. It is a proportion. You multiply it by a distance (which could be the radius of the Earth) which is NOT expressed in radian.

If you want to duplicate the (wrong and irrelevant) result of your calculator, you would need to put

image

That will give you EXACTLY what you get on the calculator; but that it is the WRONG operation.
What you REALLY need is this:

image

because the Longitude and Latitude are ALREADY in degree, and that the ‘sin’ and ‘cos’ function expect the angle to be expressed in degree, and return a result that is a pure, non dimensional number, which does not need ANY conversion whatsoever.

I have the notion that you think that the “convert degrees to radians” somehow apply to the functions “sin” and “cos”, telling them to expect an input in a different unit. That is not what it does. On the calculator, you are not changing the values, you are changing the MODE of calculation. By sliding the selector to RAD, you are changing the function called when you press the button from sinD to sinR.
What “convert degrees to radians” in Thunkable does is take whatever number it gets from the right into that number divided by 57.2957795 on the left. It does not affect anything on the right.

Remove ALL CONVERSIONS from RAD to degree everywhere in your app. Forget that there even is such a thing as RAD in the world. Everything that the functions take as input and return as output (and that includes the Location Sensor) are in degree already in Thunkable.

What you should get everywhere is “sin(31.24719)” returning 0.51873, including in your calculator, which should be in DEG mode, because 31.24719 in an angle in degree, NOT RAD.


#13

Thanks Sir.
Please maybe I am not clear,
The formula I am using impose the use of RAD as input. This is the contrary when I use sin, cos or tan in Thunkable. The functions return the value in degree as it is stated in the help pannel.
But, I have used the calculator, inputting the value (using RAD mode) and get the accurate solution. When I conceive the algorithm in Thunkable, the result is false (and it is the same if I put the calculator in DEG mode). That’s where my problem is.

This are prooves :

This is the java code used in Thunkable

// developed by Brenor Brophy - comments removed to save space
define (“meter2nm”, (1/1852));
define (“nm2meter”, 1852);
class gPoint
{
var $ellipsoid = array(//Ellipsoid name, Equatorial Radius, square of eccentricity
“Airy” =>array (6377563, 0.00667054),
“Australian National” =>array (6378160, 0.006694542),
“Bessel 1841” =>array (6377397, 0.006674372),
“Bessel 1841 Nambia” =>array (6377484, 0.006674372),
“Clarke 1866” =>array (6378206, 0.006768658),
“Clarke 1880” =>array (6378249, 0.006803511),
“Everest” =>array (6377276, 0.006637847),
“Fischer 1960 Mercury” =>array (6378166, 0.006693422),
“Fischer 1968” =>array (6378150, 0.006693422),
“GRS 1967” =>array (6378160, 0.006694605),
“GRS 1980” =>array (6378137, 0.00669438),
“Helmert 1906” =>array (6378200, 0.006693422),
“Hough” =>array (6378270, 0.00672267),
“International” =>array (6378388, 0.00672267),
“Krassovsky” =>array (6378245, 0.006693422),
“Modified Airy” =>array (6377340, 0.00667054),
“Modified Everest” =>array (6377304, 0.006637847),
“Modified Fischer 1960” =>array (6378155, 0.006693422),
“South American 1969” =>array (6378160, 0.006694542),
“WGS 60” =>array (6378165, 0.006693422),
“WGS 66” =>array (6378145, 0.006694542),
“WGS 72” =>array (6378135, 0.006694318),
“WGS 84” =>array (6378137, 0.00669438));

var $a;										// Equatorial Radius
var	$e2;									// Square of eccentricity
var	$datum;									// Selected datum
var	$Xp, $Yp;								// X,Y pixel location
var $lat, $long;							// Latitude & Longitude of the point
var $utmNorthing, $utmEasting, $utmZone;	// UTM Coordinates of the point
var	$lccNorthing, $lccEasting;				// Lambert coordinates of the point
var $falseNorthing, $falseEasting;			// Origin coordinates for Lambert Projection
var $latOfOrigin;							// For Lambert Projection
var $longOfOrigin;							// For Lambert Projection
var $firstStdParallel;						// For lambert Projection
var $secondStdParallel;						// For lambert Projection

function gPoint($datum='WGS 84')			// Default datum is WGS 84
{
	$this->a = $this->ellipsoid[$datum][0];		// Set datum Equatorial Radius
	$this->e2 = $this->ellipsoid[$datum][1];	// Set datum Square of eccentricity
	$this->datum = $datum;						// Save the datum
}
function setXY($x, $y)
{
	$this->Xp = $x; $this->Yp = $y;
}
function Xp() { return $this->Xp; }
function Yp() { return $this->Yp; }
function setLongLat($long, $lat)
{
	$this->long = $long; $this->lat = $lat;
}
function Lat() { return $this->lat; }
function Long() { return $this->long; }
function printLatLong() { printf("Latitude: %1.5f Longitude: %1.5f",$this->lat, $this->long); }
function setUTM($easting, $northing, $zone='')	// Zone is optional
{
	$this->utmNorthing = $northing;
	$this->utmEasting = $easting;
	$this->utmZone = $zone;
}
function N() { return $this->utmNorthing; }
function E() { return $this->utmEasting; }
function Z() { return $this->utmZone; }
function printUTM() { print( "Northing: ".(int)$this->utmNorthing.", Easting: ".(int)$this->utmEasting.", Zone: ".$this->utmZone); }
function setLambert($easting, $northing)
{
	$this->lccNorthing = $northing;
	$this->lccEasting = $easting;
}
function lccN() { return $this->lccNorthing; }
function lccE() { return $this->lccEasting; }
function printLambert() { print( "Northing: ".(int)$this->lccNorthing.", Easting: ".(int)$this->lccEasting); }

function convertLLtoTM($LongOrigin = NULL)
{
	$k0 = 0.9996;
	$falseEasting = 0.0;

	//Make sure the longitude is between -180.00 .. 179.9
	$LongTemp = ($this->long+180)-(integer)(($this->long+180)/360)*360-180; // -180.00 .. 179.9;
	$LatRad = deg2rad($this->lat);
	$LongRad = deg2rad($LongTemp);

	if (!$LongOrigin)
	{ // Do a standard UTM conversion - so findout what zone the point is in
		$ZoneNumber = (integer)(($LongTemp + 180)/6) + 1;
		// Special zone for South Norway
		if( $this->lat >= 56.0 && $this->lat < 64.0 && $LongTemp >= 3.0 && $LongTemp < 12.0 ) // Fixed 1.1
			$ZoneNumber = 32;
		// Special zones for Svalbard
		if( $this->lat >= 72.0 && $this->lat < 84.0 ) 
		{
			if(      $LongTemp >= 0.0  && $LongTemp <  9.0 ) $ZoneNumber = 31;
			else if( $LongTemp >= 9.0  && $LongTemp < 21.0 ) $ZoneNumber = 33;
			else if( $LongTemp >= 21.0 && $LongTemp < 33.0 ) $ZoneNumber = 35;
			else if( $LongTemp >= 33.0 && $LongTemp < 42.0 ) $ZoneNumber = 37;
		}
		
		$LongOrigin = ($ZoneNumber - 1)*6 - 180 + 3;  //+3 puts origin in middle of zone
		//compute the UTM Zone from the latitude and longitude
		$this->utmZone = sprintf("%d%s", $ZoneNumber, $this->UTMLetterDesignator());
		// We also need to set the false Easting value adjust the UTM easting coordinate
		$falseEasting = 500000.0
		;
	}

	$LongOriginRad = deg2rad($LongOrigin);

	$eccPrimeSquared = ($this->e2)/(1-$this->e2);

	$N = $this->a/sqrt(1-$this->e2*sin($LatRad)*sin($LatRad));
	$T = tan($LatRad)*tan($LatRad);
	$C = $eccPrimeSquared*cos($LatRad)*cos($LatRad);
	$A = cos($LatRad)*($LongRad-$LongOriginRad);

	$M = $this->a*((1	- $this->e2/4-3*$this->e2*$this->e2/64-5*$this->e2*$this->e2*$this->e2/256)*$LatRad 
						- (3*$this->e2/8+3*$this->e2*$this->e2/32+45*$this->e2*$this->e2*$this->e2/1024)*sin(2*$LatRad)
						+ (15*$this->e2*$this->e2/256 + 45*$this->e2*$this->e2*$this->e2/1024)*sin(4*$LatRad)
						- (35*$this->e2*$this->e2*$this->e2/3072)*sin(6*$LatRad));

	$this->utmEasting = ($k0*$N*($A+(1-$T+$C)*$A*$A*$A/6
					+ (5-18*$T+$T*$T+72*$C-58*$eccPrimeSquared)*$A*$A*$A*$A*$A/120)
					+ $falseEasting);

	$this->utmNorthing = ($k0*($M+$N*tan($LatRad)*($A*$A/2+(5-$T+9*$C+4*$C*$C)*$A*$A*$A*$A/24
				 + (61-58*$T+$T*$T+600*$C-330*$eccPrimeSquared)*$A*$A*$A*$A*$A*$A/720)));
	if($this->lat < 0)
		$this->utmNorthing += 10000000.0; //10000000 meter offset for southern hemisphere
}


#14

OK, one last time.

Because it is getting tedious.

Let’s focus only one one small part of the equation.

cos(4.15636 * Pi/180)

What does this do?
Basically, you are taking 4.15636 degree and converting it to RAD (multiply by Pi, divide by 180).
You can check this if you want.
4.15636 degree is equal to 0.0725421 RAD.
OK so far?

Now, you are in RAD mode, so you are taking the cosine of 0.0725421 RAD. Which is 0.99736997.

OK, so now what happens if we take the cosine of 4.15636 degree, in degree mode?
I get the answer of 0.99736997.

See where I am aiming at?

cosineR(0.0725421 RAD) = cosineD(4.15636 degree), because 0.0725421 RAD = 4.15636 degree.

So on the right side, when in degree mode, you DON’T multiply by Pi and divide by 180 ANYWHERE !!!

NOWHERE!

Because you do not need to convert the angles that are in degree to RAD for the sine and cosine functions, if you are using sineD and cosD, which is what you have when you are in DEG mode, and what you have in Thunkable. Because Thunkable DOES NOT have sineR and cosR functions.
That means you have to feed them degrees. Which means DO NOT convert degree to RAD.

Thunkable is NOT Java. Java has sineR and cosR. Only.
Thunkable has sineD and cosD. Only.

Your equation then becomes

0.5 * ln ( (1 + cos (4.15636) * sin(9.25026 - whatever I don’t see because it is cut off…

Look, it is simple. Wherever you see " * Pi/180 " remove it. EVERYWHERE. Go to DEG mode in your calculator, and do the math.
Your result WILL BE 527761.0009


#15

Thanks very much.
This time it is very clear and I got all the points.
I think, you have shown a good way to do this.
I will let you know if the answer is satisfactory

Thanks one more


#16

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 CBVG. Many thanks for all your contributions.

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

Thanks one more


#17

The only place where the Y calculation differs from the in principle is that Y is making use of ATAN function. ATAN returns a result in RADIAN if you are programming in Java, but in degree if you are using Thunkable. That is where the conversion is inverted, if the result is supposed to be expressed in something that has to be consistent.


#18

Do you means the use ATAN us correct here ?
Beacause the output value is too high

1


#19

ATAN returns a value in angle units. If the function is ATAN in degree, then the result will not be the same as the one expected from ATAN in radian.

What is the expectation for the result to be in?


#20

In fact, the output value of the X (eastern) and the Y (northing) that we are talking here should be in meters (m). It is some how the distance from the equator to the point . But the input data to be converted are coordinates in degree (Lat Long). Since I removed all the conversion to Radian (pi/180), the results of the values to be used in the fonction are in degree.


#21

Hi if you get it working i would love to get an example, iv been loking for it s long time