0

I am extracting data from rinex observation files and would like to convert the approximate position to latitude and longitude using Proj.4. The rinex standard states that these coordinates are WGS84, which I assume refers to the ellipsoid. What are the Proj4 parameters I need to do the conversion? If this is not possible is there another route I can take? Specifically I'm using the pyproj Python library. The correct coordinates should be close to (45.4343,-123.5516).

I am using RTKLib to process the observations, but want to automate downloading the base station/CORS files.

I did just find this code that performs the conversion. I've ported it to Python and it works well for my needs. Still if there is a way to get Proj4 to do this that would be ideal as I use that library for other projects.

Here is the Python I've tested:

In [29]: prj = pyproj.Proj(proj='geocent',ellps='WGS84',datum='WGS84')

In [30]: prj(-2478053.2006,-3736604.9657)
Out[30]: (1e+30, 1e+30)

In [31]: prj(-2478053.2006,-3736604.9657,inverse=True)
Out[31]: (-22.260730649602785, -33.5664935140521)

Here is an example header.

     2.00           OBSERVATION DATA    M (MIXED)           RINEX VERSION / TYPE
                                        20180413 230424 UTC PGM / RUN BY / DATE 
28jxmqr6                                                    MARKER NAME         
                                                            MARKER NUMBER       
                    off                                     OBSERVER / AGENCY   
100037              BE-GPS-3300         2.1.47              REC # / TYPE / VERS 
                                                            ANT # / TYPE        
 -2478053.2006 -3736604.9657  4521684.4805                  APPROX POSITION XYZ 
        0.0000        0.0000        0.0000                  ANTENNA: DELTA H/E/N
     1     1                                                WAVELENGTH FACT L1/2
     4    C1    L1    C2    L2                              # / TYPES OF OBSERV 
  2018     4    13    19     8   50.5000000     GPS         TIME OF FIRST OBS   
  2018     4    13    19    23   28.0000000     GPS         TIME OF LAST OBS    
                                                            END OF HEADER       
  • You have to use all 3 numbers, not 2. – mkennedy Apr 28 '18 at 1:06
  • Yes, that was the trick. – tharen Apr 28 '18 at 6:19
3

The geocentric coordinates can be converted to latitude, longitude, and altitude by using the pyproj.transform function. The Proj.__call__ method does not accept Z values, which are required for the transformation.

In [6]: p1 = pyproj.Proj(proj='geocent',ellps='WGS84',datum='WGS84',units='m',no_defs=True)

In [7]: p2 = pyproj.Proj(init='epsg:4326')

In [8]: pyproj.transform(p1,p2,-2478053.2006,-3736604.9657,4521684.4805)
Out[8]: (-123.55166435788638, 45.43449642806122, 450.48841851670295)

The output is equivalent to the following cs2cs input:

cs2cs +proj=geocent +ellps=WGS84 +datum=WGS84 +no_defs +units=m +to +init=epsg:4326 -f "%.8f"
-2478053.2006 -3736604.9657  4521684.4805
-123.55166436   45.43449643 450.48841852
0

I don't think a normal projection library is going to do the conversion for you.

While these coordinates are WGS84, they are cartesian (XYZ) coordinates, and not geodetic (Lat/Lon).

This forum post says that you can use the program GeoConv to do the conversion.

Maybe you should have a look at the RTKLib and see if it can proces the file for you?

  • If you know the axis directions, handedness, and orientation you can use trigonometry to work out the direction, and that gives you the lat-long. Probably the X axis goes through Null Island (lat=0, long=0) and the Y axis through (lat=0, long=plus or minus 90). The Z axis probably goes through the north pole.... – Spacedman Apr 27 '18 at 23:27
  • I think PROJ.4 has a "+proj=geocent" option now. See this issue which shows a cs2cs example. – mkennedy Apr 27 '18 at 23:28
  • The math is listed here, actually: linz.govt.nz/data/geodetic-system/coordinate-conversion/… – Alex Leith Apr 27 '18 at 23:29
  • @Spacedman 3D Cartesian to lat/lon/h on an ellipsoid generally requires iteration. – mkennedy Apr 27 '18 at 23:31
  • @mkennedy I am using the geocent option. I'll try the cs2cs example to see if that works. – tharen Apr 27 '18 at 23:59

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.