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       
2
  • You have to use all 3 numbers, not 2.
    – mkennedy
    Apr 28, 2018 at 1:06
  • Yes, that was the trick.
    – tharen
    Apr 28, 2018 at 6:19

2 Answers 2

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?

9
  • 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, 2018 at 23:27
  • I think PROJ.4 has a "+proj=geocent" option now. See this issue which shows a cs2cs example.
    – mkennedy
    Apr 27, 2018 at 23:28
  • The math is listed here, actually: linz.govt.nz/data/geodetic-system/coordinate-conversion/…
    – Alex Leith
    Apr 27, 2018 at 23:29
  • @Spacedman 3D Cartesian to lat/lon/h on an ellipsoid generally requires iteration.
    – mkennedy
    Apr 27, 2018 at 23:31
  • @mkennedy I am using the geocent option. I'll try the cs2cs example to see if that works.
    – tharen
    Apr 27, 2018 at 23:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.