Translate Cartesian coordinates to WGS84, preferably in pyproj?

Question

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.

Answer 1

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)

Answer 2

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.

Answer 3

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.