Projecting coordinates into different EPSGs gives inconsistencies

I am using 3 EPSGs: 4326 (WGS84), 3857 (Mercator) and 3035 (EU Copernicus dataset).

I am trying to easily convert coordinates from one system to any other. To do so, I wrote a `Coordinate` class with the following method, which serves to take the point's coordinates `self.x` and `self.y` given in `self.epsg` and project them into `new_epsg`:

``````   # Project point into given EPSG system
def project(self, new_epsg):
# Create transform from src epsg to dst epsg
src = osr.SpatialReference()
src.ImportFromEPSG(self.epsg)

dst = osr.SpatialReference()
dst.ImportFromEPSG(new_epsg)

transform = osr.CoordinateTransformation(src, dst)

# Create and project point
wkt = 'POINT ({} {})'.format(self.x, self.y)
point = ogr.CreateGeometryFromWkt(wkt)
point.Transform(transform)

self.epsg = new_epsg
self.x = point.GetX()
self.y = point.GetY()

return self
``````

However, when converting from an epsg to another the coordinates are always wrong.

Testing I've done:

I've tried a combination of two things:

• flipping the coordinates in input i.e. `wkt = 'POINT ({} {})'.format(self.y, self.x)` instead of `self.x, self.y`
• flipping the coordinates in output i.e. `self.x = point.GetY()` `self.y = point.GetX()`

To test my code, I've used a single point in the 3 ESPGs using https://epsg.io/. The coordinates are (x, y):

Now here are the results:

• not flipping anything: nothing works
• flipping only input: only (1) to (2) and (3) to (2) work
• flipping only output: only (2) to (1) and (2) to (3) work
• flipping both: only (1) to (3) and (3) to (1) work

So overall by flipping coordinates around I can get all cases to work but this is very case specific and I would much rather have a general one-fits-all function.

More (runnable) code is available at https://pastebin.com/vzDkC3QM

I hope I didn't just miss something stupid.

• In each case you need to consider the axis order of the crs - so in some you need to flip and in others you don't Commented Apr 11, 2020 at 15:33
• Yes, like Ian states: axis order may not be consistent between EPSG. See here pyproj4.github.io/pyproj/stable/api/transformer.html and here: proj.org/….
– GBG
Commented Apr 11, 2020 at 16:30

https://gdal.org/tutorials/osr_api_tut.html#crs-and-axis-order

If you would like a simpler solution, when creating the spatial reference objects, make sure to do:

``````crs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
``````
• To clarify - this behaviour is imposed by GDAL/OGR, and not PROJ because the former needs to express a projection in ESRI or (modified ESRI) WKT format, but a one-to-one between PROJ and WKT is not always possible. IMHO, this is an unfortunate addition. Commented May 12, 2020 at 6:02

As pointed out by Ian and GBG the issue was indeed with axes being inverted in some cases, but not all (here for 4326 and 3035).

I'll quote the links posted by GBG as they were very instructional:

Here's my updated code that takes into account inverted axes if that can help anyone: https://pastebin.com/8CFV62z8