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.?

Answer 1

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) 

Answer 2

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.