Pyproj Azimuthal Equidistant (AEQD) transformation always gives a small rotation as compared to country-specific CRS transformations

[robwijnhoven]

I am calculating building polygons from a satellite image.
This works accurately when using a country-specific CRS.
I would like to use a customer CRS, specific to the location at hand, but find that the resulting GPS coordinates are always slightly rotated.

The process I use is:

1. Take a satellite image at the location of the building I am interested in (center of image is the known GPS position)
2. Manually drawing building outlines on the image
3. Convert pixel coordinates to local meters. The origin of this system is at the center of the image and the scale from the satellite scale (meters/pixel) is known. Result is the building outline coordinates in the local coordinate system (in meters)
4. Convert the local coordinates to GPS using a Coordinate Reference System (CRS)
5. Overlay the GPS polygons in Openstreetmap using "OpenStreetMap WKT Playground" (https://clydedacruz.github.io/openstreetmap-wkt-playground/)

What I observe is that I get a small rotational difference when using different CRS transformations. Of course, some difference is expected when using different transformations as there is no single way to 'unwrap' the earth to a surface. My assumption is that a country-specific CRS (for a small country such as the Netherlands) should be very similar to using a custom CRS at the exact GPS position (because the datums are very close and the coverage of the country-specific transformation is quite small, in the order of ~400x400km).
After playing with a couple of locations all over the world, it seems that when using the Azimuthal Equidistant (AEQD) transformation, I end up with a small rotation as compared to the country-specific EPSG.

Concretely, I have two examples here from the Netherlands and France.
As you can see from the animated GIFs, the centerpoint of the projection is exactly the same, but there is a minor difference in rotation.
The country-specific EPSGs are most accurate (or at least seem to be).

1. Netherlands (EPSG:28992) [![Example projection of buildings in Netherlands][2]][2]
   
2. France (EPSG:2154) [![Example projection of buildings in France][1]][1]
   

More details:

• AEQD CRS string: "+proj=aeqd +lat_0={latitude} +lon_0={longitude} +ellps=WGS84" units=m"
• Pyproj version used: 3.6.1, which uses PROJ 9.3.0 (sept 1st, 2023)
• Tried multiple buildings at different locations w.r.t. country-specific CRS bounds, so one building at the center (close to the datum), and at the most north-west point, and one at the very west). The behaviour is the same for all locations.
• Tried different countries as well, similar behaviour.

Conclusion: after spending a significant amount of time debugging this, I think this might be a bug in the PROJ (or pyproj) library, which could result from a numerical issue. However, I want to rule out that I missed something in the process. Also posted this in the stackexchange geo forum (https://gis.stackexchange.com/questions/470746/pyproj-azimuthal-equidistant-aeqd-transformation-always-gives-a-small-rotation)

[djhoese]

I have no idea if this is expected, but...

My suggestion is that you should make a https://stackoverflow.com/help/minimal-reproducible-example with real coordinates and show the transformation. You could probably pick a single coordinate that shows the unexpected rotation (like one of the corners of your polygon) when compared with the other transformation. You could also include the center/reference coordinate as one of these example coordinates. Once you have shown that you get the unexpected results in pyproj, convert your example to using proj on the command line. Once you prove that PROJ has the same issue then post your example as a PROJ bug report. With a real example the PROJ maintainers should be able to tell you pretty quickly whether something is going wrong or if this is expected.

[mdsumner]

I don't see the problem, 28992 is Oblique Stereographic and has a different centre to the one you chose. The numbers and the shapes will be different in different coordinate systems.

aeqd is preserving minute accuracy of distance from that nominated centre point to any other point, sterea is preserving shape (in a local sense) at the expense of area and distance