ECEF XYZ to LAT/LON

[Natalie Angel]

I have searched the whole PB forum but did not find any real answer to my problem.

I'm working on a GPS module that gives the Longitude, Latitude and Altitude in ECEF format where XYZ is given in SINT32 Little-endian in meters.

Now I have converted the X and Y and got as an example -946051 for X and 374257 for Y

When I use Atn(Y/X) I get value like 89.9998468824687 instead of 99.013540 (note that the results are examples and might not be the real answer) which is 10 degrees off my position.

I have tried several ways of converting this but without luck so if anyone have worked with ECEF (Earth-Centered, Earth-Fixed) it would be good to know how to convert ECEF to Longitude, Latitude, Altitude (LLA) as it's used by many GPS receivers.

Thank you

[Paul Dixon]

Attila,
could your problem be the use of ATN without taking care over which quadrant you're in?

Check out this discussion:
http://www.powerbasic.com/support/pb...highlight=atn2

Paul.

[John Gleason]

Here is a pdf doc describing some conversions that may help:
http://www.microem.ru/pages/u_blox/t...G1-X-00006.pdf

[Natalie Angel]

Thank you for this. I tried the example but did not get the answer I needed.

Here is a pdf doc describing some conversions that may help...

I use the same document but still...

It says Longitude = arctan(Y/X) and that's what I try to do but it's off about 10 degrees. I have test software that save the data to file in LON/LAT so I know it works, it's just my calculations which is wrong.

I receive the data in ASCII where I take out the 4 bytes, I then convert them to HEX and add them together in reverse order as they are given in Little-endian.

[John Gleason]

Attila, perhaps if you post example code showing your calculations, someone may be able help, possibly even I.

[Natalie Angel]

OK I have tested one thing...

First of all to make the code clear in HEX I needed to format the hex valuse as they containd "5" and "1" by using

Format(Hex$(Asc(Mid$(TempGetBuffer,x,1))), "00")

and as it's 4 bytes I did so with all 4 bytes (it's on another machine so too long to type as I work with 3 screens in one time)

The 4 HEX values are added by using:

lECEFx = CLng("&H" & Format(Hex$(Asc(Mid$(TempGetBuffer,16,1))), "00") & Format(Hex$(Asc(Mid$(TempGetBuffer,15,1))), "00") & Format(Hex$(Asc(Mid$(TempGetBuffer,14,1))), "00") & Format(Hex$(Asc(Mid$(TempGetBuffer,13,1))), "00"))

Could also be in a loop ... same with Y

Where the HEX in Little-endian is 7E:90:F1:FF and if I use val("&HFFF1907E") for X then I get -946050

For Y it's D7:05:5B:20 in HEX and val("&H205B05D7") gives 542836183

For Z it's 17:2F:1F:20 in HEX and val("&H201F2F17") gives 538914583 (Z is Altitude)

Now if I take 542836183/-946050 I get -573.792276306749

Now type this into http://www.analyzemath.com/Calculato...alculator.html and you get -89.90015 wich is the same as if I use

Atn(542836183/-946050) * 57.2957795130823

so it's still wrong by 10 degrees or so??? as it should be somewhere around 99.013540 even if I use

Atn(1 - (542836183/-946050)) * 57.2957795130823 to get positive number

[Arie Verheul]

Just some remarks

Hi Attila,

I am not sure if it is the cause of your problem, but still some remarks.

1. If you find a longtitude of -89.9 degrees, this should be translated as 180 - 89,9 = 90.1 degrees.
(The +X axis points towards Greenwich).

2. If the result angle is close to 90 degrees it is more accurate to swap X and Y and calculate Atn(X/Y), relative to the 90 degree plane.

3. I do not quite understand what you mean with the last line, but it is definitely wrong.
Subtracting 1 from a tangens does not produce a sensible result.

Atn(1 - (542836183/-946050)) * 57.2957795130823 to get positive number

Arie Verheul

[Michael Mattias]

lECEFx =  _
       CLng("&H" _ 
    & Format(Hex$(Asc(Mid$(TempGetBuffer,16,1))), "00") _
    & Format(Hex$(Asc(Mid$(TempGetBuffer,15,1))), "00") _ 
    & Format(Hex$(Asc(Mid$(TempGetBuffer,14,1))), "00") _ 
   & Format(Hex$(Asc(Mid$(TempGetBuffer,13,1))), "00"))

There has got to be a better way than this.

What do bytes 13-16 look like and what should result be?

Wait a minute, I just broke that up to put it on multiple lines, and that's not valid syntax. (You can't do CLNG against a string; FORMAT is not a valid intrinsic function).

[Chris Holbrook]

And what is this "format" function? Anything like format$?

[Dale Yarker]

I put your numbers into calculator at:
http://www.ngs.noaa.gov/cgi-bin/xyz_getgp.prl
and got:
N895330.63466 E1582458.64092 532558790.9471 GRS80
About 89° like your code produced. The difference is probably due to your use of longs, while NOAA uses floating type. Where did you get the value that is 10° off?

From a quick search an "SINT32 Little-endian" is exactly the same as a PB long. Instead of the quadruple conversion, just use a pointer to grab the four bytes into a long? Something like:

#dim all

function MyFunc(ECEF as string) as long
  pMeters as long pointer
  Xmeters as long
  Ymeters as long
  Zmeters as long

  pMeters = strptr(ECEF) + 12 '13th byte
  Xmeters = @pMeters

  'repeat for Y and Z with appropriate offsets to pointer

  'do the calculation

end function

Cheers,

[Natalie Angel]

Thank you Dale,

In fact the GPS module comes with test software where the data is saved to a file and also because the location is in Chiang Mai, Thailand I know it's right.

I found a website that got the code in C++ so maybe that could help?

http://www.gpstk.org/doxygen/ECEF_8hpp-source.html and
http://www.gpstk.org/doxygen/ECEF_8cpp-source.html

Anyhow the formula is in the PDF that John pointed to so if this could be converted into PB code it would be good. I work with PB/Win 8.04

The Position (ECEF) X = +/- 2^31 given in SINT32 (32-bit signed integer) given in Little-endian and in Meters.

Now my input string is; "xxxxž9.>~�ñÿ× [ /- ¢" where the “x” is just other data.

If I convert the string (cutting out the x’es) I get this two strings;

Decimal:
158:57:46:62:126:144:241:255:215:5:91:32:23:47:31:32:32:32:1:32:162

HEX:
9E-39-2E-3E-7E-90-F1-FF-D7-05-5B-20-17-2F-1F-20-20-20-01-20-A2

where
9E:39 = GPS time (in binary function)
2E:3E = TOD
7E:90:F1:FF = ECEF X
D7:05:5B:20 = ECEF Y
17:2F:1F:20 = ECEF Z
20:20 = POI
01:20 = Information (in binary format)
A2 = ?

[Arie Verheul]

What are you doing up there ?

Just another thought.

There is of course a good quick check for your data as, ignoring the ellipticty of the earth,
for any (X,Y,Z) it must approximately hold true that :

R = Sqr(X^2 + Y^2 + Z^2)

where R = 4E7 / 2Pi.

Your coordinates give R = 5.42E8 , which suggest that you must be in orbit.

Arie Verheul

  
#Dim All
#Break On
 
Function PBMain
   
    Local X,Y,Z,R,Pi As Single
   
    Pi = 3.14159
  
    R = 4E7 / (2*Pi)
  
    Print Using$("Earth radius  ##.##^^^^",R)
    Print
  
    X = -94605
    Y = 54283618
    Z = 538914583
   
    R = Sqr(X^2 + Y^2 + Z^2)
       
    Print Using$("  X           ##.##^^^^",X)
    Print Using$("  Y           ##.##^^^^",Y)
    Print Using$("  Z           ##.##^^^^",Z)
    Print
    Print Using$("Your altitude ##.##^^^^",R)
  
    WaitKey$
  
End Function

[Michael Mattias]

>And what is this "format" function? Anything like format$?

I think that may be "brand M" BASIC code. That would explain the "legality" of passing a string to the CLNG function (I think) and the 'intrinsicness' of the FORMAT [no suffixed '$'] function.

[Natalie Angel]

It was in VB format as I was not on the PC with PB loaded.

[Natalie Angel]

OK I've got the answer now and will post the DLL code in the Source Code section soon.

In fact it was not easy but you know...