1

I've been looking at the following solution to a similar problem Trilateration using 3 latitude and longitude points, and 3 distances However with my problem the 3 points are spread over a large distance and the distance is a distance-over-ground (ignoring altitude). So the curvature of the earth comes into play.

Example of one of the data points is...[lat,long,km]

51.505348978413,-0.11270052661132,75.639168 51.4845086249581,-3.16947466601561,173.809152 55.9568033803799,-3.20266968478392,465.100416

This can be visualised at http://www.freemaptools.com/radius-around-point.htm by pasting it into the section "CSV Upload - [latitude,longitude,radius(km)] per line"

Using the solution in Trilateration using 3 latitude and longitude points, and 3 distances I get the result..

51.8247121391 -0.802763718095

which isn't bad but it's about 14 miles out. I'm really wanting this down to ~5 miles of error.

I changed the ECEF calculation to one from here http://www.mathworks.co.uk/matlabcentral/fileexchange/7942-covert-lat-lon-alt-to-ecef-cartesian/content/lla2ecef.m and used an altitude of 0.

However, this is way way off the mark with...

56.2518456429 -0.813853131422

Intermediate results for xA,yA,zA...etc. are

lla2ecef

xA = [3978.79751502] yA = [-7.82628596545] zA = [5342.89431232]

xB = [3974.53105834] yB = [-220.086727248] zB = [5342.88794881]

xC = [3573.7981748] yC = [-199.973379186] zC = [5344.22411569]

original method

xA = [3965.56758195] yA = [-7.8002627162] zA = [4986.36680448]

xB = [3961.32002974] yB = [-219.355176278] zB = [4984.92406448]

xC = [3561.02859617] yC = [-199.258852046] zC = [5279.1109334]

So... the Z values for lla2ecef method are quite different. I'm guessing this is what's causing the issue, something to do with fixing the altitude to 0 and the fact my distances are distance-over-ground.

Am I going to need to do some projections instead?

My brain is about melted for today, can anyone shed light on ways to tackle this problem?

Improve this question
6
  • Using the method at gis.stackexchange.com/questions/40660/… (but replacing the Euclidean Norm by the ellipsoidal GeoDistance), I obtain the solution (52.0292, -0.8107). Is that close enough?
    – whuber
    Jun 11, 2013 at 15:54
  • 101.873m NEE from Dansteed Way, Milton Keynes - reverse geocoded gives you the street google.co.uk/…
    – Mapperz ♦
    Jun 11, 2013 at 16:24
  • whuber - When you say replacing the Euclidean Norm by GeoDistance I assume you are using the Mathematica implementation? I'll probably be porting the solution to C++ so would like to solve the problem in generic way.
    – YMDW
    Jun 11, 2013 at 16:50
  • The Mathematica code is generic, as explained in that solution. You can use your own implementation of ellipsoidal distance and your own nonlinear fitting procedure. Both are widely available on many platforms. You can also adopt an approach advocated by @cffk (check his answers) in which an initial guess of the solution is used as the center of an oblique equidistant projection, the points are reprojected, and the guess is updated using Euclidean calculations. You will want to update the guess using the least-squares method described in my solution.
    – whuber
    Jun 11, 2013 at 17:13
  • I'm a little lost but it gives me some direction, thanks.
    – YMDW
    Jun 11, 2013 at 17:38

1 Answer 1

Reset to default
1

Using the Non Linear Fitting method as provided by whuber in (Trilateration algorithm for n amount of points) I used Matlab to solve the problem. 

P = [51.505348978413,-0.11270052661132;
    51.4845086249581,-3.16947466601561;
    55.9568033803799,-3.20266968478392];

R = [75.639168, 173.809152, 465.100416];
R = R*1000; %convert km to m

opts = statset('MaxIter',1000);
beta0 = [53.374, -1.394];
BETA=nlinfit(P,R,'myfun',beta0,opts)
  • P are the 3 trilateration points (lat, lon)
  • R are the ranges from my sample point to points in P distance-over-ground
  • beta0 are the initial guesses for a position that meets the condition defined by P and R

myfun is in a separate file

function result = myfun(beta, X)

ellipsoid = [6378137, 8.1819190842622e-2];
result(1) = distance(51.505348978413, -0.11270052661132, beta(1), beta(2), ellipsoid);
result(2) = distance(51.4845086249581, -3.16947466601561, beta(1), beta(2), ellipsoid);
result(3) = distance(55.9568033803799, -3.20266968478392, beta(1), beta(2), ellipsoid);

I'm not sure what X does in this case, but myfun needs to have two inputs (beta & X). The initial guess returns three distances to my points P. result is compared to R until a maximum of 1000 iterations (or some error is reduced below a limit - not sure what this limit is or where to set it yet).

Once the values for beta have been found that results in the same values of R the final values are stored in BETA

Improve this answer

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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