different projection results by using python vs c++ implementation of Proj

[gabrielwaibel]

Problem

I am integrating the proj libray with python and c++. Using the same coordinate and the same crs conversion outputs slightly different projected values from easting and northing.

Code

My python code using pyproj version 3.5.0.

from pyproj import Transformer

# Define the transformer with source CRS as WGS84 (EPSG:4326) and target CRS as UTM Zone 33N (EPSG:32633)
transformer = Transformer.from_crs("epsg:4326", "EPSG:2056")

# Coordinates in WGS84 (latitude, longitude)
lat, lon = 47.0074167, 8.1472778

# Convert to UTM Zone 33N
x, y = transformer.transform(lat, lon)

print(f"Latitude {lat}, Longitude {lon}")
print(f"Converted coordinates: X = {x}, Y = {y}")

This outputs:

Latitude 47.0074167, Longitude 8.1472778
Converted coordinates: X = 2653889.0601541456, Y = 1206505.6025497625

My code in c++ taken from the tutorial on the website using proj version 9.3 from source.

  auto dbContext = DatabaseContext::create();
  auto authFactory = AuthorityFactory::create(dbContext, std::string());
  auto coord_op_ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0);

  auto authFactoryEPSG = AuthorityFactory::create(dbContext, "EPSG");
  auto sourceCRS = authFactoryEPSG->createCoordinateReferenceSystem("4326");
  auto targetCRS = authFactoryEPSG->createCoordinateReferenceSystem("2056");

  auto list = CoordinateOperationFactory::create()->createOperations(sourceCRS, targetCRS, coord_op_ctxt);

  assert(!list.empty());
  PJ_CONTEXT* ctx = proj_context_create();
  auto transformer = list[0]->coordinateTransformer(ctx);

  PJ_COORD c = {{
      latitude,   // latitude in degree
      longitude,  // longitude in degree
      altitude,   // z ordinate. unused
      HUGE_VAL    // time ordinate. unused
  }};
  c = transformer->transform(c);

  // Display result
  // std::cout << std::fixed << std::setprecision(16);
  // std::cout << "Easting: " << c.v[0] << std::endl;   // should be 426857.988
  // std::cout << "Northing: " << c.v[1] << std::endl;  // should be 5427937.523

  // Destroy execution context
  proj_context_destroy(ctx);

which outputs:

Latitude: 47.007416700000000 Longitude: 8.147277799999999
easting: 2653889.054882026743144, northing: 1206505.592695795698091

Question

Where does this difference of 2653889.0601541456 -> 2653889.054882026743144 and 1206505.6025497625 -> 1206505.592695795698091 come from?

The difference is too much to be some floating point precision. Also if I am using a different target crs as epsg:31467 the error between python and c++ even increases.

[rouault]

This transformation is Z dependent. By trials, I assume in your C++ code you use an altitude of ~ 430 m See:

$ echo 47.0074167 8.1472778 430 | bin/cs2cs -d 5 EPSG:4326 EPSG:2056
2653889.05481   1206505.59257 430.00000

vs

$ echo 47.0074167 8.1472778 0 | bin/cs2cs -d 5 EPSG:4326 EPSG:2056
2653889.06015   1206505.60255 0.00000

[gabrielwaibel]

ok, thanks for the fast answer! That explains it.

but why is EPSG:2056 altitude dependent? I assumed that epsg code 2056 does not have a dependency in z. On top of this transformation I would like to use another one (EPSG:3855) which actually transforms z according to a geoid model.

Should I then use the height 0.0 for the first conversion or the one from my gps module in wgs84?

[rouault]

The transformation involves a Helmert 3D transformation (in geocentric coordinates), hence the 2D position needs to be extended to 3D to perform a geographic to geocentric conversion. https://www.iogp.org/wp-content/uploads/2019/09/373-07-02.pdf mentions at paragraph "4.1.4 Geographic 3D to 2D conversion": "The reverse conversion, from 2D to 3D, is indeterminate. It is however a requirement when a geographic 2D coordinate reference system is to be transformed using a geocentric method which is 3-dimensional (see section 4.3.1). In practice an artificial ellipsoidal height is created and appended to the geographic 2D coordinate reference system to create a geographic 3D coordinate reference system referenced to the same geodetic datum. The assumed ellipsoidal height is usually either set to the gravity-related height of a position in a compound coordinate reference system, or set to zero. As long as the height chosen is within a few kilometres of sea level, no error will be induced into the horizontal position resulting from the later geocentric transformation; the vertical coordinate will however be meaningless."

"no error" here means that the difference in results is far below the accuracy of the transformation. This transformation has an advertised accuracy of 1 metre, so the differences of a few cm you observer are considered to be unsignificant.

[gabrielwaibel]

ok makes sense that it needs to depend on z.

But do you know how people perform a second geoid conversion in practice? Basically applying two epsg codes, one for x,y and one for z.

If I do a first conversion (lat, lon, alt) to (x, y) and a second geoid transformation (lat, lon, alt) to (z) this results in (x, y, z). But if I do the geoid transformation first, the overall result would be different since the height influences the first transformation.

[jjimenezshaw]

But do you know how people perform a second geoid conversion in practice?

Apply it directly, for instance, from EPSG:4326 to EPSG:2056+5729