Get lat lon for a pos when only dx dy is known from an other point

Question

The problem I have is getting the latitude and longitude for a position when I know one origin position and the offset to the other point (in meter)

I want to use following function (or any other function that gives me a correct lat2 and lon2):

import geographiclib
geographiclib.geodesic.Geodesic.Direct(self, lat1, lon1, azi1, s12)

I have some example data:

lon1 = 11.62113333
lat1 = 55.9862

dx = -51659.25 #meter
dy = -33702.33 #meter

This is the result that I'm hoping to achieve:

azi1 = -120.95109978727244
s12 = 61691.57175978693         
lat2 = 55.69834 
lon2 = 10.77969

When I use a UTM2latlon converter I'm getting a position that is way off.. I think that the coordinate system that is used to calculate dx and dy is ESPG:5596.

Since the distances is great and the distance have taken into account for the earthsPythagoras Theorem isn't applicable to calculate s12 and azi1. Any suggestions on functions etc.?

Do you mean EPSG:5596? Maybe the GDAL or pyproj can do the conversion. See http://gis.stackexchange.com/questions/78838/how-to-convert-projected-coordinates-to-lat-lon-using-python for an example. — , Aug 18 '15 at 18:41

score 0 · Answer 1 · answered Aug 18 '15 at 22:26

You can get to a result close-ish to your expected values in two steps:

from geographiclib.geodesic import Geodesic

lon1 = 11.62113333
lat1 = 55.9862

dx = -51659.25 #meter
dy = -33702.33 #meter

tmp = Geodesic.WGS84.Direct(
    lat1, lon1,
    90, # Go East ...
    dx) #         ... (negatively, so actually West)

destination = Geodesic.WGS84.Direct(
    tmp['lat2'], tmp['lon2'],
    0,  # Go North ...
    dy) #          ... (negatively, so actually South)

lat2, lon2 = destination['lat2'], destination['lon2']
# 55.680721562111955, 10.793499275609594

But this calculation is still conceptually wrong.

If we assume the earth to be a sphere (rather than an ellipsoid), the second step (Go North/South) is fine, as you're moving on a meridian, which is a great circle. Any way along a great circle is the shortest (on-surface) path between it's start and end point, unless it goes more than halfway around the sphere. A geodesic is a shortest path on an ellipsoid, so this fits.

For the first step, however, we'd want to go dx meters East (or -dx meters West). Following a geodesic, heading West initially won't do that, because the parallels (except for the equator) aren't great circles / geodesics. So the calculated path will drift away from the parallel, which you can easily see by looking at the intermediate location:

tmp
# {'a12': -0.4645520407433718,
#  'azi1': 90.0,
#  'azi2': 89.31397924368152,
#  'lat1': 55.9862,
#  'lat2': 55.98342231490044,
#  'lon1': 11.62113333,
#  'lon2': 10.793499275609594,
#  's12': -51659.25}

The latitude there (tmp['lat2']) isn't our initial latitude and the azimut (tmp['azi2']) changed away from 90° (straight West), too!

Thus, geodesic calculations are probably not the correct way to solve your problem. Also note that the order or the steps (whether you go North/South first and then East/West, or East/West first and then North/South) will significantly affect the result when calculated like this.

Which order would be 'correct' (if any) is determined by the projection / cartesian coordinate reference system referenced by dx and dy.

score 0 · Accepted Answer · answered Sep 14 '15 at 06:16

Thanks for the try!

Unfortunately I worked with some old algorithm so the "correct answer" was a more trivial function.

def melat(lat):
    return  degrees( log ( tan ( radians( lat / 2 + 45 ) ) ))

def latit(mlt):
    return ( degrees(atan ( exp ( radians(mlt) ) )) - 45 ) * 2

def XY2LatLong(x, y):
    EQUATOR_MINUTE_LENGTH = 1851.8518519
    GeoPoint = namedtuple('GeoPoint', ('lat', 'lon'))
    Point = namedtuple('Point', ('x', 'y'))
    base = GeoPoint(55.98619572, 11.62113936)
    factor = cos( radians(base.lat) ) * EQUATOR_MINUTE_LENGTH * 60.
    return GeoPoint(
        latit(melat(base.lat) + y / factor),
        base.lon + x / factor)

Get lat lon for a pos when only dx dy is known from an other point

2 Answers2