0

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.

2

2 Answers 2

1

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)
1
  • 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.
    – Ralph Tee
    Commented May 12, 2020 at 6:02
0

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

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.