0

I have some Cartesian coordinates:

x = 3862525.511
y = -109738.522

I need to translate them into WGS84 coordinates, preferably using pyproj. Is this possible, or have I misunderstood the nature of Cartesian coordinates? I don't have any height info.

cartesian = ?
wgs84 = pyproj.Proj(init='epsg:4326')
new_x, new_y = pyproj.transform(cartesian, wgs84, x, y)

In case it's useful, these Cartesian coordinates weren't created by me, but I'm pretty sure they were created by copying and pasting OSGB36 eastings and northings on this site.

If it can't be done in pyproj, I can use gdal.

  • I think these are 3d coordinates with (0,0,0) at the centre of the earth, the poles will be (0,0, +/-R), Null Island (lat-long 0,0) is at maybe (0,-R,R) and so on. It depends how the axes are oriented, but you do need the z coordinate. – Spacedman Jun 25 '18 at 20:01
  • Hm. Maybe I can just set it to the average UK height above sea level (162m according to Wikipedia), or rather the equivalent in Cartesian coordinates? The WGS84 coordinates don't have to be perfectly accurate, +/- 10 metres is probably ok. – Richard Jun 25 '18 at 20:13
  • Actually I think you can treat the X,Y coordinates as being on an azimuthal, zenithal projection centred on the N pole. Hmmm.... – Spacedman Jun 25 '18 at 20:57
  • I tested the coordinate pair that you provided with an azimuthal orthographic north pole projection, and it falls near the Molson Coors building in Burton-on-Trent, so the assumption makes sense. This projection could be a way to convert your points by deprojection. Or if you wish to do it mathematically, you would need to take into account the ellipsoidal shape of the Earth. – FSimardGIS Jul 6 '18 at 0:02
0

Here's something in python. Note it doesn't use the Z coordinate:

import math

R = 6360000 # earth radius

def cc_to_ll(xy):
    (x, y, z) = xy # does not use `z`
    r = math.sqrt(x**2 + y**2)
    long =180 * math.atan2(y,x)/math.pi
    lat = 180 * math.acos(r/R)/math.pi
    return (lat, long)

To test on a known point, I've take the location of the statue of Eric Morecambe in Morecambe, which is approximately at sea level:

eric_google = (54.072878,-2.8680394)

The online converter returns the following OSGB and Cartesian coordinates for these decimal lat-long coordinates:

eric_osgb = (343296.36, 464452.65)
eric_xyz = (3745951.536, -187666.857, 5141507.615)

Test script is:

print "Eric according to google ", eric_google
print "Eric in cartesian ", eric_xyz
print "Eric converted ",cc_to_ll(eric_xyz)

The code above returns:

Eric according to google  (54.072878, -2.8680394)
Eric in cartesian  (3745951.536, -187666.857, 5141507.615)
Eric converted  (53.862520806716375, -2.868039396573805)

The converted longitude is exact to 6dp, but latitude is off by a bit. This is because the code assumes a spherical earth, and a point on the surface of the sphere. Hence its dependent on the radius, R, and the ellipsoidal shape used in the EPSG 4326 coordinate system.

I don't think you can get 10m precision without a Z coordinate in the cartesian system. The Z axis is everywhere parallel to the earth's axis, and in the UK that is at an angle of about 37 degrees, so a difference in Z coordinate of 100m will result in a shift of about 30m.

I think a better solution might be obtained by working out the intersection of the Z axis-aligned point with the EPSG 4326 ellipsoid rather than a sphere.

  • Amazing, what an effort, thanks so much for trying! Hm, as you say, Eric is perfect in latitude terms, but off quite a way in longitude terms. Sad... – Richard Jun 26 '18 at 10:16
0

Here is the code for the ellipsoid (WGS84), assuming the points are on the surface of the ellipsoid, it solves the z value and the latitude value :

import math
def xy_to_ll(x, y):
    a = 6378137.0 # WGS84 Semimajor axis
    f = 1/298.257223563 # WGS84 Flatteing
    b = a*(1-f) # Semiminor axis
    z = math.sqrt(b**2*(1-((x**2+y**2)/a**2))) # Solving for the positive z value on the ellipsoid
    p = math.sqrt(x**2+y**2) # Distance from the z-axis
    Lat = math.atan(z/(p*(1-f)**2))*180/math.pi # Solving the latitude value
    Long = math.atan2(y,x)*180/math.pi
    return (Lat, Long)

Here is the result with the test point values in @Spacedman's example :

>>> xy_to_ll(3745951.536, -187666.857)
(54.072878003746759, -2.8680393965738049)

The latitude matches the original input latitude. Of course, if the points are not on the surface and were converted to cartesian XYZ but you only have the X and Y, then it is not possible to adequately solve, unless you could somehow estimate the height of the point.

0

Actually if you want to use pyproj, it is quite straightforward:

import pyproj

x = 3862525.511
y = -109738.522

osgb36 = pyproj.Proj(init='epsg:27700')
wgs84 = pyproj.Proj(init='epsg:4326')

>>> pyproj.transform(osgb36, wgs84, x, y)
(38.52757978082899, 40.955405746532314)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.