Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

Given a list of entities (or GeoPositions) and their distances from a particular location, how do I find that location's GeoPosition?

dist = Map[{#, GeoDistance[$GeoLocation, #]}&, GeoNearest["University", $GeoLocation, 5]]

enter image description here

How do I find this point?

GeoGraphics[{Thick, GeoMarker[], GeoCircle @@@ dist}, ImageSize -> 500]

enter image description here

If you want to try with this data:

{{Entity["University", "VanderCookCollegeOfMusic149639"],   Quantity[1.407900938586476, "Miles"]},
 {Entity["University", "IllinoisInstituteOfTechnology145725"], Quantity[1.5607147441441935, "Miles"]},
 {Entity["University", "UniversityOfIllinoisAtChicago145600"],   Quantity[1.6865509085430637, "Miles"]},
 {Entity["University", "ShimerCollege148849"], Quantity[1.7240912464169733, "Miles"]},
 {Entity["University", "IllinoisCollegeOfOptometry145628"],   Quantity[1.7283592978609772, "Miles"]}}
share|improve this question

1 Answer 1

up vote 4 down vote accepted

I'm definitely not proud of this, but here's a brute force method that is pretty accurate.

First let's define some functions, cheating and using the built in GeoDistance function. We'll also try to minimize the error between the distances we're using as inputs and the distances we're calculating based on our guess {latitude, longitude}:

geodist[pos1 : {_?NumericQ, _?NumericQ}, pos2 : {_?NumericQ, _?NumericQ}] := GeoDistance[pos1, pos2]

findpt[latlons : {{_, _} ..}, distances_?ListQ] := Module[{solpos}, 
  solpos = {$lat, $lon} /. 
  NMinimize[
   Total[
    Abs@QuantityMagnitude[
     (geodist[{$lat, $lon}, #] & /@ latlons) - distances]], 
   {{$lat, Max[Min[latlons[[All, 1]]] - 10, -90], Min[Max[latlons[[All, 1]]] + 10, 90]},
    {$lon, Max[Min[latlons[[All, 2]]] - 10, -180], Min[Max[latlons[[All, 2]]] + 10, 180]}}, 
 Method -> "NelderMead"][[2]];
GeoPosition[solpos]
]

The only thing to check is to make sure the distances calculated by GeoDistance are in the same units as the input distances.

So let's test:

inputlatlon = {39.72, -104.81};
dist = Map[{#, GeoDistance[GeoPosition[inputlatlon], #]} &, 
  GeoNearest["University", GeoPosition[inputlatlon], 5]]
GeoGraphics[{Thick, GeoMarker[GeoPosition[inputlatlon]], 
GeoCircle @@@ dist}, ImageSize -> 500]

enter image description here

findpt[QuantityMagnitude[LatitudeLongitude /@ dist[[All, 1]]], dist[[All, 2]]]

GeoPosition[{39.72, -104.81}]

GeoDistance[%, GeoPosition[inputlatlon]]

Quantity[2.62342*10^-6, "Feet"]

Hmmm, not too bad.

UnitConvert[%, "Angstroms"]

Quantity[7996.18, "Angstroms"]

Damn, 8000 Angstroms...

share|improve this answer
1  
Thanks! :) I was holding out to see if someone might come up with a symbolic solution, but it seems like this is the way to go. Angstroms... Mathematica's geodesy is almost as precise as the NSA's spy satellites. –  mfvonh yesterday

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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