Inconsistent conversion results

[originlake]

Hi all,

I'm facing a strange conversion issue, that the conversion results with the same coordinates and the EPSG code are getting slightly different results. Hope I can get some help here. The pyproj version is 3.1.0

I want to convert 351541.762, 503637.841 from EPSG:3168 to EPSG:32647. Below is what I did

from pyproj import Proj, Transformer
proj1 = Proj(3168)
proj2 = Proj(32647)
# Approach 1
lat, lon = proj1(351541.762, 503637.841, inverse=True)
res1 = proj2(lat, lon)
print(res1) # get  (739157.8248517661, 503381.4067311792)
# Approach 2
transformer = Transformer.from_proj(proj1, proj2)
res2 = transformer.transform(351541.762, 503637.841)
print(res2) # get (739004.0222534146, 503348.46853315743)
# res1 != res2 

As you can see, the two approaches get different results, they are close but have around 150m distance.
The first approach's result matches the conversion in epsg.io (http://epsg.io/transform#s_srs=3168&t_srs=32647&x=351541.7620000&y=503637.8410000), the 2nd one matches the conversion in QGIS(3.18).

I'm so confused with the results, not sure which one is correct.

Any suggestion would be appreciated

[snowman2]

This might provide some insight: https://pyproj4.github.io/pyproj/stable/gotchas.html#proj-not-a-generic-latitude-longitude-to-projection-converter

[originlake]

Thanks, really appreciate your help

[snowman2]

Try this: echo 351541.762 503637.841 | cs2cs 'EPSG:3168' 'EPSG:32647'

[originlake]

Thanks, this works. Kind of curious what is missing in the conversion using the Proj4 string. Can I consider it the proj4 string has limitations of representing datum shifts?

[snowman2]

When using the EPSG codes you can see the datum information:

>>> CRS("EPSG:3168")
<Projected CRS: EPSG:3168>
Name: Kertau (RSO) / RSO Malaya (m)
Axis Info [cartesian]:
- E[east]: Easting (metre)
- N[north]: Northing (metre)
Area of Use:
- name: Malaysia - West Malaysia onshore.
- bounds: (99.59, 1.21, 104.6, 6.72)
Coordinate Operation:
- name: Rectified Skew Orthomorphic Malaya Grid (metres)
- method: Hotine Oblique Mercator (variant A)
Datum: Kertau (RSO)
- Ellipsoid: Everest 1830 (RSO 1969)
- Prime Meridian: Greenwich

>>> CRS('EPSG:32647')
<Projected CRS: EPSG:32647>
Name: WGS 84 / UTM zone 47N
Axis Info [cartesian]:
- E[east]: Easting (metre)
- N[north]: Northing (metre)
Area of Use:
- name: Between 96°E and 102°E, northern hemisphere between equator and 84°N, onshore and offshore. China. Indonesia. Laos. Malaysia - West Malaysia. Mongolia. Myanmar (Burma). Russian Federation. Thailand.
- bounds: (96.0, 0.0, 102.0, 84.0)
Coordinate Operation:
- name: UTM zone 47N
- method: Transverse Mercator
Datum: World Geodetic System 1984
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

You can lose CRS information when converting to PROJ strings: https://proj.org/faq.html#what-is-the-best-format-for-describing-coordinate-reference-systems

[originlake]

Thanks again for your help!

I found that the datum shift is missing for proj.4 string of EPSG:3168 based on this discussion, after adding +towgs84=-11,851,5, it's working as expected.