您的位置：首页 > Web前端 > JavaScript

Javascript 地图经度纬度之间的距离的算法(轉)

2012-01-12 12:39 447 查看

Calculatedistance,bearingandmorebetweenLatitude/Longitudepoints

http://www.movable-type.co.uk/scripts/latlong.html

Thispagepresentsavarietyofcalculationsforlatitude/longitudepoints,withtheformul?andcodefragmentsforimplementingthem.

Alltheseformul?areforcalculationsonthebasisofasphericalearth(ignoringellipsoidaleffects)–whichisaccurateenough*formostpurposes…[Infact,theearthisveryslightlyellipsoidal;usingasphericalmodelgiveserrorstypicallyupto0.3%–seenotesforfurtherdetails].

Entertheco-ordinatesintothetextboxestotryoutthecalculations.Avarietyofformatsareaccepted,principally:

Distance

Thisusesthe‘haversine’formulatocalculatethegreat-circledistancebetweentwopoints–thatis,theshortestdistanceovertheearth’ssurface–givingan‘as-the-crow-flies’distancebetweenthepoints(ignoringanyhills,ofcourse!).

Haversineformula:

a=sin²(Δlat/2)+cos(lat1).cos(lat2).sin²(Δlong/2)
c=2.atan2(√a,√(1−a))
d=R.c
whereRisearth’sradius(meanradius=6,371km);
notethatanglesneedtobeinradianstopasstotrigfunctions!
JavaScript:
varR=6371;//kmvardLat=(lat2-lat1).toRad();vardLon=(lon2-lon1).toRad();varlat1=lat1.toRad();varlat2=lat2.toRad();vara=Math.sin(dLat/2)*Math.sin(dLat/2)+Math.sin(dLon/2)*Math.sin(dLon/2)*Math.cos(lat1)*Math.cos(lat2);varc=2*Math.atan2(Math.sqrt(a),Math.sqrt(1-a));vard=R*c;

Thehaversineformula1‘remainsparticularlywell-conditionedfornumericalcomputationevenatsmalldistances’–unlikecalculationsbasedonthesphericallawofcosines.The‘versedsine’is1-cosθ,andthe‘half-versed-sine’(1-cosθ)/2=sin²(θ/2)asusedabove.ItwaspublishedbyRWSinnottinSkyandTelescope,1984,thoughknownaboutformuchlongerbynavigators.(Forthecurious,cistheangulardistanceinradians,andaisthesquareofhalfthechordlengthbetweenthepoints).A(surprisinglymarginal)performanceimprovementcanbeobtained,ofcourse,byfactoringoutthetermswhichgetsquared.

SphericalLawofCosines

Infact,whenSinnottpublishedthehaversineformula,computationalprecisionwaslimited.Nowadays,JavaScript(andmostmoderncomputers&languages)useIEEE75464-bitfloating-pointnumbers,whichprovide15significantfiguresofprecision.Withthisprecision,thesimplesphericallawofcosinesformula(cosc=cosacosb+sinasinbcosC)giveswell-conditionedresultsdowntodistancesassmallasaround1metre.(Notethatthegeodeticformofthelawofcosinesisrearrangedfromthecanonicalonesothatthelatitudecanbeuseddirectly,ratherthanthecolatitude).

Thismakesthesimplerlawofcosinesareasonable1-linealternativetothehaversineformulaformanypurposes.Thechoicemaybedrivenbycodingcontext,availabletrigfunctions(indifferentlanguages),etc.

Sphericallaw
ofcosines:		d=acos(sin(lat1).sin(lat2)+cos(lat1).cos(lat2).cos(long2−long1)).R
JavaScript:
varR=6371;//kmvard=Math.acos(Math.sin(lat1)*Math.sin(lat2)+Math.cos(lat1)*Math.cos(lat2)*Math.cos(lon2-lon1))*R;

Excel:=ACOS(SIN(lat1)*SIN(lat2)+COS(lat1)*COS(lat2)*COS(lon2-lon1))*6371
(Notethathereandinallsubsequentcodefragments,forsimplicityIdonotshowconversionsfromdegreestoradians;seebelowforcompleteversions).

Equirectangularapproximation

Ifperformanceisanissueandaccuracylessimportant,forsmalldistancesPythagoras’theoremcanbeusedonanequirectangularprojection:*

x=Δlon.cos(lat)
y=Δlat
d=R.√x²+y²
JavaScript:
varx=(lon2-lon1)*Math.cos((lat1+lat2)/2);vary=(lat2-lat1);vard=Math.sqrt(x*x+y*y)*R;

(lat/loninradians!)
Thisusesjustonetrigandonesqrtfunction–asagainsthalf-a-dozentrigfunctionsforcoslaw,and7trigs+2sqrtsforhaversine.Accuracyissomewhatcomplex:alongmeridianstherearenoerrors,otherwisetheydependondistance,bearing,andlatitude,butaresmallenoughformanypurposes*(andoftentrivialcomparedwiththesphericalapproximationitself).

Bearing

Ingeneral,yourcurrentheadingwillvaryasyoufollowagreatcirclepath(orthodrome);thefinalheadingwilldifferfromtheinitialheadingbyvaryingdegreesaccordingtodistanceandlatitude(ifyouweretogofromsay35°N,45°E(Baghdad)to35°N,135°E(Osaka),youwouldstartonaheadingof60°andenduponaheadingof120°!).

Thisformulaisfortheinitialbearing(sometimesreferredtoasforwardazimuth)whichiffollowedinastraightlinealongagreat-circlearcwilltakeyoufromthestartpointtotheendpoint:1

Formula:θ=atan2(sin(Δlong).cos(lat2),
cos(lat1).sin(lat2)−sin(lat1).cos(lat2).cos(Δlong))
JavaScript:
vary=Math.sin(dLon)*Math.cos(lat2);varx=Math.cos(lat1)*Math.sin(lat2)-Math.sin(lat1)*Math.cos(lat2)*Math.cos(dLon);varbrng=Math.atan2(y,x).toDeg();

Excel:=ATAN2(COS(lat1)*SIN(lat2)-SIN(lat1)*COS(lat2)*COS(lon2-lon1),
SIN(lon2-lon1)*COS(lat2))
*NotethatExcelreversestheargumentstoATAN2–seenotesbelow
Sinceatan2returnsvaluesintherange-π...+π(thatis,-180°...+180°),tonormalisetheresulttoacompassbearing(intherange0°...360°,with-vevaluestransformedintotherange180°...360°),converttodegreesandthenuse(θ+360)%360,where%ismodulo.

Forfinalbearing,simplytaketheinitialbearingfromtheendpointtothestartpointandreverseit(usingθ=(θ+180)%360).

Midpoint

Thisisthehalf-waypointalongagreatcirclepathbetweenthetwopoints.1

Formula:Bx=cos(lat2).cos(Δlong)
By=cos(lat2).sin(Δlong)
latm=atan2(sin(lat1)+sin(lat2),√((cos(lat1)+Bx)²+By²))
lonm=lon1+atan2(By,cos(lat1)+Bx)
JavaScript:
varBx=Math.cos(lat2)*Math.cos(dLon);varBy=Math.cos(lat2)*Math.sin(dLon);varlat3=Math.atan2(Math.sin(lat1)+Math.sin(lat2),Math.sqrt((Math.cos(lat1)+Bx)*(Math.cos(lat1)+Bx)+By*By));varlon3=lon1+Math.atan2(By,Math.cos(lat1)+Bx);

Justastheinitialbearingmayvaryfromthefinalbearing,themidpointmaynotbelocatedhalf-waybetweenlatitudes/longitudes;themidpointbetween35°N,45°Eand35°N,135°Eisaround45°N,90°E.

Destinationpointgivendistanceandbearingfromstartpoint

Givenastartpoint,initialbearing,anddistance,thiswillcalculatethedestinationpointandfinalbearingtravellingalonga(shortestdistance)greatcirclearc.

Formula:lat2=asin(sin(lat1)*cos(d/R)+cos(lat1)*sin(d/R)*cos(θ))
lon2=lon1+atan2(sin(θ)*sin(d/R)*cos(lat1),cos(d/R)−sin(lat1)*sin(lat2))
θisthebearing(inradians,clockwisefromnorth);
d/Ristheangulardistance(inradians),wheredisthedistancetravelledandRistheearth’sradius
JavaScript:
varlat2=Math.asin(Math.sin(lat1)*Math.cos(d/R)+Math.cos(lat1)*Math.sin(d/R)*Math.cos(brng));varlon2=lon1+Math.atan2(Math.sin(brng)*Math.sin(d/R)*Math.cos(lat1),Math.cos(d/R)-Math.sin(lat1)*Math.sin(lat2));

Excel:lat2:=ASIN(SIN(lat1)*COS(d/R)+COS(lat1)*SIN(d/R)*COS(brng))
lon2:=lon1+ATAN2(COS(d/R)-SIN(lat1)*SIN(lat2),SIN(brng)*SIN(d/R)*COS(lat1))
Forfinalbearing,simplytaketheinitialbearingfromtheendpointtothestartpointandreverseit(usingθ=(θ+180)%360).

Intersectionoftwopathsgivenstartpointsandbearings

Thisisarathermorecomplexcalculationthanmostothersonthispage,butI'vebeenaskedforitanumberoftimes.SeebelowfortheJavaScript.

Formula:d12=2.asin(√(sin²(Δlat/2)+cos(lat1).cos(lat2).sin²(Δlon/2)))
φ1=acos(sin(lat2)−sin(lat1).cos(d12)/sin(d12).cos(lat1))
φ2=acos(sin(lat1)−sin(lat2).cos(d12)/sin(d12).cos(lat2))

ifsin(lon2−lon1)>0
θ12=φ1,θ21=2.π−φ2
else
θ12=2.π−φ1,θ21=φ2

α1=(θ1−θ12+π)%2.π−π
α2=(θ21−θ2+π)%2.π−π

α1=|α1|
α2=|α2|

α3=acos(−cos(α1).cos(α2)+sin(α1).sin(α2).cos(d12))
d13=atan2(sin(d12).sin(α1).sin(α2),cos(α2)+cos(α1).cos(α3))
lat3=asin(sin(lat1).cos(d13)+cos(lat1).sin(d13).cos(θ1))
Δlon13=atan2(sin(θ1).sin(d13).cos(lat1),cos(d13)−sin(lat1).sin(lat3))
lon3=(lon1+Δlon13+π)%2.π−π

wherelat1,lon1,θ1:1stpoint&bearing
lat2,lon2,θ2:2ndpoint&bearing
lat3,lon3:intersectionpoint

%=mod,||=abs

note–ifsin(α1)=0andsin(α2)=0:infinitesolutions
ifsin(α1).sin(α2)<0:ambiguoussolution
thisformulationisnotalwayswell-conditionedformeridionalorequatoriallines
Notethiscanalsobesolvedusingvectorsratherthantrigonometry:

Foreachpointφ,λ(lat=φ,lon=λ),wecandefineaunitvectorpointingtoitfromthecentreoftheearth:u{x,y,z}=[cosφ·cosλ,cosφ·sinλ,sinφ](takingx=0º,y=90º,z=north–notethattheseformulædependonconventionusedfordirectionsandhandedness)

Andforanygreatcircledefinedbytwopoints,wecandefineaunitvectorNnormaltotheplaneofthecircle:N(u1,u2)=(u1×u2)/||u1×u2||where×isthevectorcrossproduct,and||u||thenorm(lengthofthevector)

Thevectorrepresentingtheintersectionofthetwogreatcirclesisthenui=±N(N(u1,u2),N(u3,u4))

WecanthengetthelatitudeandlongitudeofPibyφ=atan2(uz,sqrt(ux²+uy²)),λ=atan2(uy,ux)

Theantipodalintersectionpointis(-φ,λ+π)

Cross-trackdistance

Here’sanewone:I’vesometimesbeenaskedaboutdistanceofapointfromagreat-circlepath(sometimescalledcrosstrackerror).

Formula:dxt=asin(sin(d13/R)*sin(θ13−θ12))*R
whered13isdistancefromstartpointtothirdpoint
θ13is(initial)bearingfromstartpointtothirdpoint
θ12is(initial)bearingfromstartpointtoendpoint
Ristheearth’sradius
JavaScript:
vardXt=Math.asin(Math.sin(d13/R)*Math.sin(brng13-brng12))*R;

Here,thegreat-circlepathisidentifiedbyastartpointandanendpoint–dependingonwhatinitialdatayou’reworkingfrom,youcanusetheformulæabovetoobtaintherelevantdistanceandbearings.Thesignofdxttellsyouwhichsideofthepaththethirdpointison.

Thealong-trackdistance,fromthestartpointtotheclosestpointonthepathtothethirdpoint,is

Formula:dat=acos(cos(d13/R)/cos(dxt/R))*R
whered13isdistancefromstartpointtothirdpoint
dxtiscross-trackdistance
Ristheearth’sradius
JavaScript:
vardAt=Math.acos(Math.cos(d13/R)/Math.cos(dXt/R))*R;



Closestpointtothepoles
And:‘Clairaut’sformula’willgiveyouthemaximumlatitudeofagreatcirclepath,givenabearingandlatitudeonthegreatcircle:

Formula:latmax=acos(abs(sin(θ)*cos(lat)))
JavaScript:
varlatMax=Math.acos(Math.abs(Math.sin(brng)*Math.cos(lat)));



Rhumblines
A‘rhumbline’(orloxodrome)isapathofconstantbearing,whichcrossesallmeridiansatthesameangle.

Sailorsusedto(andsometimesstill)navigatealongrhumblinessinceitiseasiertofollowaconstantcompassbearingthantobecontinuallyadjustingthebearing,asisneededtofollowagreatcircle.RhumblinesarestraightlinesonaMercatorProjectionmap(alsohelpfulfornavigation).

Rhumblinesaregenerallylongerthangreat-circle(orthodrome)routes.Forinstance,LondontoNewYorkis4%longeralongarhumblinethanalongagreatcircle–importantforaviationfuel,butnotparticularlytosailingvessels.NewYorktoBeijing–closetothemostextremeexamplepossible(thoughnotsailable!)–is30%longeralongarhumbline.

Distance/bearing
Theseformulægivethedistanceand(constant)bearingbetweentwopoints.

Formula:Δφ=ln(tan(lat2/2+π/4)/tan(lat1/2+π/4))[=the‘stretched’latitudedifference]
ifE:Wline,q=cos(lat1)
otherwise,q=Δlat/Δφ
d=√(Δlat²+q².Δlon²).R[pythagoras]
θ=atan2(Δlon,Δφ)
wherelnisnaturallog,Δlonistakingshortestroute(<180º),andRistheearth’sradius
JavaScript:
vardPhi=Math.log(Math.tan(lat2/2+Math.PI/4)/Math.tan(lat1/2+Math.PI/4));varq=(!isNaN(dLat/dPhi))?dLat/dPhi:Math.cos(lat1);//E-WlinegivesdPhi=0
//ifdLonover180°takeshorterrhumbacross180°meridian:if(Math.abs(dLon)>Math.PI){
dLon=dLon>0?-(2*Math.PI-dLon):(2*Math.PI+dLon);}vard=Math.sqrt(dLat*dLat+q*q*dLon*dLon)*R;varbrng=Math.atan2(dLon,dPhi);

Destination
Givenastartpointandadistancedalongconstantbearingθ,thiswillcalculatethedestinationpoint.Ifyoumaintainaconstantbearingalongarhumbline,youwillgraduallyspiralintowardsoneofthepoles.

Formula:α=d/R(angulardistance)
lat2=lat1+α.cos(θ)
Δφ=ln(tan(lat2/2+π/4)/tan(lat1/2+π/4))[=the‘stretched’latitudedifference]
ifE:Wlineq=cos(lat1)
otherwiseq=Δlat/Δφ
Δlon=α.sin(θ)/q
lon2=(lon1+Δlon+π)%2.π−π
wherelnisnaturallogand%ismodulo,Δlonistakingshortestroute(<180°),andRistheearth’sradius
JavaScript:
vardLat=d*Math.cos(brng);varlat2=lat1+dLat;vardPhi=Math.log(Math.tan(lat2/2+Math.PI/4)/Math.tan(lat1/2+Math.PI/4));varq=(!isNaN(dLat/dPhi))?dLat/dPhi:Math.cos(lat1);//E-WlinegivesdPhi=0vardLon=d*Math.sin(brng)/q;//checkforsomedaftbuggergoingpastthepole,normaliselatitudeifsoif(Math.abs(lat2)>Math.PI/2)lat2=lat2>0?Math.PI-lat2:-(Math.PI-lat2);
lon2=(lon1+dLon+Math.PI)%(2*Math.PI)-Math.PI;

Mid-point
Thisformulaforcalculatingthe‘loxodromicmidpoint’,thepointhalf-wayalongarhumblinebetweentwopoints,isduetoRobertHillandCliveTooth1(thxAxel!).

Formula:latm=(lat1+lat2)/2
f1=tan(π/4+lat1/2)
f2=tan(π/4+lat2/2)
fm=tan(π/4+latm/2)
lonm=[(lon2−lon1).ln(fm)+lon1.ln(f2)−lon2.ln(f1)]/ln(f2/f1)
wherelnisnaturallog
JavaScript:
varlat3=(lat1+lat2)/2;varf1=Math.tan(Math.PI/4+lat1/2);varf2=Math.tan(Math.PI/4+lat2/2);varf3=Math.tan(Math.PI/4+lat3/2);varlon3=((lon2-lon1)*Math.log(f3)+lon1*Math.log(f2)-lon2*Math.log(f1))/Math.log(f2/f1);


Usingthescriptsinwebpages
Usingthesescriptsinwebpageswouldbesomethinglikethefollowing:

<scriptsrc="latlon.js">/*Latitude/Longitudeformulae*/</script>
<scriptsrc="geo.js">/*Geodesyrepresentationconversions*/</script>
...
<form>
Lat1:<inputtype="text"name="lat1"id="lat1">Lon1:<inputtype="text"name="lon1"id="lon1">
Lat2:<inputtype="text"name="lat2"id="lat2">Lon2:<inputtype="text"name="lon2"id="lon2">
<buttononClick="varp1=newLatLon(Geo.parseDMS(f.lat1.value),Geo.parseDMS(f.lon1.value));
varp2=newLatLon(Geo.parseDMS(f.lat2.value),Geo.parseDMS(f.lon2.value));
alert(p1.distanceTo(p2));">Calculatedistance</button>
</form>


IfyouusejQuery,thecodecanbeseparatedfromtheHTML:

<scriptsrc="http://ajax.googleapis.com/ajax/libs/jquery/1.6.0/jquery.min.js"></script>
<scriptsrc="latlon.js">/*Latitude/Longitudeformulae*/</script>
<scriptsrc="geo.js">/*Geodesyrepresentationconversions*/</script>
<script>
$(document).ready(function(){
$('#calc-dist').click(function(){
varp1=newLatLon(Geo.parseDMS($('#lat1').val()),Geo.parseDMS($('#lon1').val()));
varp2=newLatLon(Geo.parseDMS($('#lat2').val()),Geo.parseDMS($('#lon2').val()));
$('#result-distance').html(p1.distanceTo(p2)+'km');
});
});
</script>
...
<form>
Lat1:<inputtype="text"name="lat1"id="lat1">Lon1:<inputtype="text"name="lon1"id="lon1">
Lat2:<inputtype="text"name="lat2"id="lat2">Lon2:<inputtype="text"name="lon2"id="lon2">
<buttonid="calc-dist">Calculatedistance</button>
<outputid="result-distance"></output>
</form>
--------------------------------------------------------------------------------


/*-----------------------------------------------*/
/*Latitude/longitudesphericalgeodesyformulae&scripts(c)ChrisVeness2002-2011*/
/*-www.movable-type.co.uk/scripts/latlong.html*/
/**/
/*Sampleusage:*/
/*varp1=newLatLon(51.5136,-0.0983);*/
/*varp2=newLatLon(51.4778,-0.0015);*/
/*vardist=p1.distanceTo(p2);//inkm*/
/*varbrng=p1.bearingTo(p2);//indegreesclockwisefromnorth*/
/*...etc*/
/*-----------------------------------------------*/

/*-----------------------------------------------*/
/*Notethatminimalerrorcheckingisperformedinthisexamplecode!*/
/*-----------------------------------------------*/

/**
*Createsapointontheearth'ssurfaceatthesuppliedlatitude/longitude
*
*@constructor
*@param{Number}lat:latitudeinnumericdegrees
*@param{Number}lon:longitudeinnumericdegrees
*@param{Number}[rad=6371]:radiusofearthifdifferentvalueisrequiredfromstandard6,371km
*/
functionLatLon(lat,lon,rad){
if(typeof(rad)=='undefined')rad=6371;//earth'smeanradiusinkm
//onlyacceptnumbersorvalidnumericstrings
this._lat=typeof(lat)=='number'?lat:typeof(lat)=='string'&&lat.trim()!=''?+lat:NaN;
this._lon=typeof(lon)=='number'?lon:typeof(lon)=='string'&&lon.trim()!=''?+lon:NaN;
this._radius=typeof(rad)=='number'?rad:typeof(rad)=='string'&&trim(lon)!=''?+rad:NaN;
}

/**
*Returnsthedistancefromthispointtothesuppliedpoint,inkm
*(usingHaversineformula)
*
*from:Haversineformula-R.W.Sinnott,"VirtuesoftheHaversine",
*SkyandTelescope,vol68,no2,1984
*
*@param{LatLon}point:Latitude/longitudeofdestinationpoint
*@param{Number}[precision=4]:noofsignificantdigitstouseforreturnedvalue
*@returns{Number}Distanceinkmbetweenthispointanddestinationpoint
*/
LatLon.prototype.distanceTo=function(point,precision){
//default4sigfigsreflectstypical0.3%accuracyofsphericalmodel
if(typeofprecision=='undefined')precision=4;

varR=this._radius;
varlat1=this._lat.toRad(),lon1=this._lon.toRad();
varlat2=point._lat.toRad(),lon2=point._lon.toRad();
vardLat=lat2-lat1;
vardLon=lon2-lon1;

vara=Math.sin(dLat/2)*Math.sin(dLat/2)+
Math.cos(lat1)*Math.cos(lat2)*
Math.sin(dLon/2)*Math.sin(dLon/2);
varc=2*Math.atan2(Math.sqrt(a),Math.sqrt(1-a));
vard=R*c;
returnd.toPrecisionFixed(precision);
}

/**
*Returnsthe(initial)bearingfromthispointtothesuppliedpoint,indegrees
*seehttp://williams.best.vwh.net/avform.htm#Crs*
*@param{LatLon}point:Latitude/longitudeofdestinationpoint
*@returns{Number}InitialbearingindegreesfromNorth
*/
LatLon.prototype.bearingTo=function(point){
varlat1=this._lat.toRad(),lat2=point._lat.toRad();
vardLon=(point._lon-this._lon).toRad();

vary=Math.sin(dLon)*Math.cos(lat2);
varx=Math.cos(lat1)*Math.sin(lat2)-
Math.sin(lat1)*Math.cos(lat2)*Math.cos(dLon);
varbrng=Math.atan2(y,x);

return(brng.toDeg()+360)%360;
}

/**
*Returnsfinalbearingarrivingatsupplieddestinationpointfromthispoint;thefinalbearing
*willdifferfromtheinitialbearingbyvaryingdegreesaccordingtodistanceandlatitude
*
*@param{LatLon}point:Latitude/longitudeofdestinationpoint
*@returns{Number}FinalbearingindegreesfromNorth
*/
LatLon.prototype.finalBearingTo=function(point){
//getinitialbearingfromsuppliedpointbacktothispoint...
varlat1=point._lat.toRad(),lat2=this._lat.toRad();
vardLon=(this._lon-point._lon).toRad();

vary=Math.sin(dLon)*Math.cos(lat2);
varx=Math.cos(lat1)*Math.sin(lat2)-
Math.sin(lat1)*Math.cos(lat2)*Math.cos(dLon);
varbrng=Math.atan2(y,x);

//...&reverseitbyadding180°
return(brng.toDeg()+180)%360;
}

/**
*Returnsthemidpointbetweenthispointandthesuppliedpoint.
*seehttp://mathforum.org/library/drmath/view/51822.htmlforderivation
*
*@param{LatLon}point:Latitude/longitudeofdestinationpoint
*@returns{LatLon}Midpointbetweenthispointandthesuppliedpoint
*/
LatLon.prototype.midpointTo=function(point){
lat1=this._lat.toRad(),lon1=this._lon.toRad();
lat2=point._lat.toRad();
vardLon=(point._lon-this._lon).toRad();

varBx=Math.cos(lat2)*Math.cos(dLon);
varBy=Math.cos(lat2)*Math.sin(dLon);

lat3=Math.atan2(Math.sin(lat1)+Math.sin(lat2),
Math.sqrt((Math.cos(lat1)+Bx)*(Math.cos(lat1)+Bx)+By*By));
lon3=lon1+Math.atan2(By,Math.cos(lat1)+Bx);
lon3=(lon3+3*Math.PI)%(2*Math.PI)-Math.PI;//normaliseto-180..+180º

returnnewLatLon(lat3.toDeg(),lon3.toDeg());
}

/**
*Returnsthedestinationpointfromthispointhavingtravelledthegivendistance(inkm)onthe
*giveninitialbearing(bearingmayvarybeforedestinationisreached)
*
*seehttp://williams.best.vwh.net/avform.htm#LL*
*@param{Number}brng:Initialbearingindegrees
*@param{Number}dist:Distanceinkm
*@returns{LatLon}Destinationpoint
*/
LatLon.prototype.destinationPoint=function(brng,dist){
dist=typeof(dist)=='number'?dist:typeof(dist)=='string'&&dist.trim()!=''?+dist:NaN;
dist=dist/this._radius;//convertdisttoangulardistanceinradians
brng=brng.toRad();//
varlat1=this._lat.toRad(),lon1=this._lon.toRad();

varlat2=Math.asin(Math.sin(lat1)*Math.cos(dist)+
Math.cos(lat1)*Math.sin(dist)*Math.cos(brng));
varlon2=lon1+Math.atan2(Math.sin(brng)*Math.sin(dist)*Math.cos(lat1),
Math.cos(dist)-Math.sin(lat1)*Math.sin(lat2));
lon2=(lon2+3*Math.PI)%(2*Math.PI)-Math.PI;//normaliseto-180..+180º

returnnewLatLon(lat2.toDeg(),lon2.toDeg());
}

/**
*Returnsthepointofintersectionoftwopathsdefinedbypointandbearing
*
*seehttp://williams.best.vwh.net/avform.htm#Intersection*
*@param{LatLon}p1:Firstpoint
*@param{Number}brng1:Initialbearingfromfirstpoint
*@param{LatLon}p2:Secondpoint
*@param{Number}brng2:Initialbearingfromsecondpoint
*@returns{LatLon}Destinationpoint(nullifnouniqueintersectiondefined)
*/
LatLon.intersection=function(p1,brng1,p2,brng2){
brng1=typeofbrng1=='number'?brng1:typeofbrng1=='string'&&trim(brng1)!=''?+brng1:NaN;
brng2=typeofbrng2=='number'?brng2:typeofbrng2=='string'&&trim(brng2)!=''?+brng2:NaN;
lat1=p1._lat.toRad(),lon1=p1._lon.toRad();
lat2=p2._lat.toRad(),lon2=p2._lon.toRad();
brng13=brng1.toRad(),brng23=brng2.toRad();
dLat=lat2-lat1,dLon=lon2-lon1;

dist12=2*Math.asin(Math.sqrt(Math.sin(dLat/2)*Math.sin(dLat/2)+
Math.cos(lat1)*Math.cos(lat2)*Math.sin(dLon/2)*Math.sin(dLon/2)));
if(dist12==0)returnnull;

//initial/finalbearingsbetweenpoints
brngA=Math.acos((Math.sin(lat2)-Math.sin(lat1)*Math.cos(dist12))/
(Math.sin(dist12)*Math.cos(lat1)));
if(isNaN(brngA))brngA=0;//protectagainstrounding
brngB=Math.acos((Math.sin(lat1)-Math.sin(lat2)*Math.cos(dist12))/
(Math.sin(dist12)*Math.cos(lat2)));

if(Math.sin(lon2-lon1)>0){
brng12=brngA;
brng21=2*Math.PI-brngB;
}else{
brng12=2*Math.PI-brngA;
brng21=brngB;
}

alpha1=(brng13-brng12+Math.PI)%(2*Math.PI)-Math.PI;//angle2-1-3
alpha2=(brng21-brng23+Math.PI)%(2*Math.PI)-Math.PI;//angle1-2-3

if(Math.sin(alpha1)==0&&Math.sin(alpha2)==0)returnnull;//infiniteintersections
if(Math.sin(alpha1)*Math.sin(alpha2)<0)returnnull;//ambiguousintersection

//alpha1=Math.abs(alpha1);
//alpha2=Math.abs(alpha2);
//...EdWilliamstakesabsofalpha1/alpha2,butseemstobreakcalculation?

alpha3=Math.acos(-Math.cos(alpha1)*Math.cos(alpha2)+
Math.sin(alpha1)*Math.sin(alpha2)*Math.cos(dist12));
dist13=Math.atan2(Math.sin(dist12)*Math.sin(alpha1)*Math.sin(alpha2),
Math.cos(alpha2)+Math.cos(alpha1)*Math.cos(alpha3))
lat3=Math.asin(Math.sin(lat1)*Math.cos(dist13)+
Math.cos(lat1)*Math.sin(dist13)*Math.cos(brng13));
dLon13=Math.atan2(Math.sin(brng13)*Math.sin(dist13)*Math.cos(lat1),
Math.cos(dist13)-Math.sin(lat1)*Math.sin(lat3));
lon3=lon1+dLon13;
lon3=(lon3+3*Math.PI)%(2*Math.PI)-Math.PI;//normaliseto-180..+180º

returnnewLatLon(lat3.toDeg(),lon3.toDeg());
}

/*-----------------------------------------------*/

/**
*Returnsthedistancefromthispointtothesuppliedpoint,inkm,travellingalongarhumbline
*
*seehttp://williams.best.vwh.net/avform.htm#Rhumb*
*@param{LatLon}point:Latitude/longitudeofdestinationpoint
*@returns{Number}Distanceinkmbetweenthispointanddestinationpoint
*/
LatLon.prototype.rhumbDistanceTo=function(point){
varR=this._radius;
varlat1=this._lat.toRad(),lat2=point._lat.toRad();
vardLat=(point._lat-this._lat).toRad();
vardLon=Math.abs(point._lon-this._lon).toRad();

vardPhi=Math.log(Math.tan(lat2/2+Math.PI/4)/Math.tan(lat1/2+Math.PI/4));
varq=(!isNaN(dLat/dPhi))?dLat/dPhi:Math.cos(lat1);//E-WlinegivesdPhi=0
//ifdLonover180°takeshorterrhumbacross180°meridian:
if(dLon>Math.PI)dLon=2*Math.PI-dLon;
vardist=Math.sqrt(dLat*dLat+q*q*dLon*dLon)*R;

returndist.toPrecisionFixed(4);//4sigfigsreflectstypical0.3%accuracyofsphericalmodel
}

/**
*Returnsthebearingfromthispointtothesuppliedpointalongarhumbline,indegrees
*
*@param{LatLon}point:Latitude/longitudeofdestinationpoint
*@returns{Number}BearingindegreesfromNorth
*/
LatLon.prototype.rhumbBearingTo=function(point){
varlat1=this._lat.toRad(),lat2=point._lat.toRad();
vardLon=(point._lon-this._lon).toRad();

vardPhi=Math.log(Math.tan(lat2/2+Math.PI/4)/Math.tan(lat1/2+Math.PI/4));
if(Math.abs(dLon)>Math.PI)dLon=dLon>0?-(2*Math.PI-dLon):(2*Math.PI+dLon);
varbrng=Math.atan2(dLon,dPhi);

return(brng.toDeg()+360)%360;
}

/**
*Returnsthedestinationpointfromthispointhavingtravelledthegivendistance(inkm)onthe
*givenbearingalongarhumbline
*
*@param{Number}brng:BearingindegreesfromNorth
*@param{Number}dist:Distanceinkm
*@returns{LatLon}Destinationpoint
*/
LatLon.prototype.rhumbDestinationPoint=function(brng,dist){
varR=this._radius;
vard=parseFloat(dist)/R;//d=angulardistancecoveredonearth'ssurface
varlat1=this._lat.toRad(),lon1=this._lon.toRad();
brng=brng.toRad();

varlat2=lat1+d*Math.cos(brng);
vardLat=lat2-lat1;
vardPhi=Math.log(Math.tan(lat2/2+Math.PI/4)/Math.tan(lat1/2+Math.PI/4));
varq=(!isNaN(dLat/dPhi))?dLat/dPhi:Math.cos(lat1);//E-WlinegivesdPhi=0
vardLon=d*Math.sin(brng)/q;
//checkforsomedaftbuggergoingpastthepole
if(Math.abs(lat2)>Math.PI/2)lat2=lat2>0?Math.PI-lat2:-(Math.PI-lat2);
lon2=(lon1+dLon+3*Math.PI)%(2*Math.PI)-Math.PI;

returnnewLatLon(lat2.toDeg(),lon2.toDeg());
}

/*-----------------------------------------------*/

/**
*Returnsthelatitudeofthispoint;signednumericdegreesifnoformat,otherwiseformat&dp
*asperGeo.toLat()
*
*@param{String}[format]:Returnvalueas'd','dm','dms'
*@param{Number}[dp=0|2|4]:Noofdecimalplacestodisplay
*@returns{Number|String}Numericdegreesifnoformatspecified,otherwisedeg/min/sec
*
*@requiresGeo
*/
LatLon.prototype.lat=function(format,dp){
if(typeofformat=='undefined')returnthis._lat;

returnGeo.toLat(this._lat,format,dp);
}

/**
*Returnsthelongitudeofthispoint;signednumericdegreesifnoformat,otherwiseformat&dp
*asperGeo.toLon()
*
*@param{String}[format]:Returnvalueas'd','dm','dms'
*@param{Number}[dp=0|2|4]:Noofdecimalplacestodisplay
*@returns{Number|String}Numericdegreesifnoformatspecified,otherwisedeg/min/sec
*
*@requiresGeo
*/
LatLon.prototype.lon=function(format,dp){
if(typeofformat=='undefined')returnthis._lon;

returnGeo.toLon(this._lon,format,dp);
}

/**
*Returnsastringrepresentationofthispoint;formatanddpasperlat()/lon()
*
*@param{String}[format]:Returnvalueas'd','dm','dms'
*@param{Number}[dp=0|2|4]:Noofdecimalplacestodisplay
*@returns{String}Comma-separatedlatitude/longitude
*
*@requiresGeo
*/
LatLon.prototype.toString=function(format,dp){
if(typeofformat=='undefined')format='dms';

if(isNaN(this._lat)||isNaN(this._lon))return'-,-';

returnGeo.toLat(this._lat,format,dp)+','+Geo.toLon(this._lon,format,dp);
}

/*-----------------------------------------------*/

//----extendNumberobjectwithmethodsforconvertingdegrees/radians

/**Convertsnumericdegreestoradians*/
if(typeof(Number.prototype.toRad)==="undefined"){
Number.prototype.toRad=function(){
returnthis*Math.PI/180;
}
}

/**Convertsradianstonumeric(signed)degrees*/
if(typeof(Number.prototype.toDeg)==="undefined"){
Number.prototype.toDeg=function(){
returnthis*180/Math.PI;
}
}

/**
*Formatsthesignificantdigitsofanumber,usingonlyfixed-pointnotation(noexponential)
*
*@param{Number}precision:Numberofsignificantdigitstoappearinthereturnedstring
*@returns{String}Astringrepresentationofnumberwhichcontainsprecisionsignificantdigits
*/
if(typeof(Number.prototype.toPrecisionFixed)==="undefined"){
Number.prototype.toPrecisionFixed=function(precision){
if(isNaN(this))return'NaN';
varnumb=this<0?-this:this;//can'ttakelogof-venumber...
varsign=this<0?'-':'';

if(numb==0){//can'ttakelogofzero,justformatwithprecisionzeros
varn='0.';
while(precision--)n+='0';
returnn
}

varscale=Math.ceil(Math.log(numb)*Math.LOG10E);//noofdigitsbeforedecimal
varn=String(Math.round(numb*Math.pow(10,precision-scale)));
if(scale>0){//addtrailingzeros&insertdecimalasrequired
l=scale-n.length;
while(l-->0)n=n+'0';
if(scale<n.length)n=n.slice(0,scale)+'.'+n.slice(scale);
}else{//prefixdecimalandleadingzerosifrequired
while(scale++<0)n='0'+n;
n='0.'+n;
}
returnsign+n;
}
}

/**Trimswhitespacefromstring(q.v.blog.stevenlevithan.com/archives/faster-trim-javascript)*/
if(typeof(String.prototype.trim)==="undefined"){
String.prototype.trim=function(){
returnString(this).replace(/^\s\s*/,'').replace(/\s\s*$/,'');
}
}

/*-----------------------------------------------*/


/*-----------------------------------------------*/
/*Geodesyrepresentationconversionfunctions(c)ChrisVeness2002-2011*/
/*-www.movable-type.co.uk/scripts/latlong.html*/
/**/
/*Sampleusage:*/
/*varlat=Geo.parseDMS('51°28′40.12″N');*/
/*varlon=Geo.parseDMS('000°00′05.31″W');*/
/*varp1=newLatLon(lat,lon);*/
/*-----------------------------------------------*/

varGeo={};//Geonamespace,representingstaticclass

/**
*Parsesstringrepresentingdegrees/minutes/secondsintonumericdegrees
*
*Thisisveryflexibleonformats,allowingsigneddecimaldegrees,ordeg-min-secoptionally
*suffixedbycompassdirection(NSEW).Avarietyofseparatorsareaccepted(eg3º37'09"W)
*orfixed-widthformatwithoutseparators(eg0033709W).Secondsandminutesmaybeomitted.
*(Noteminimalvalidationisdone).
*
*@param{String|Number}dmsStr:Degreesordeg/min/secinvarietyofformats
*@returns{Number}Degreesasdecimalnumber
*@throws{TypeError}dmsStrisanobject,perhapsDOMobjectwithout.value?
*/
Geo.parseDMS=function(dmsStr){
if(typeofdeg=='object')thrownewTypeError('Geo.parseDMS-dmsStris[DOM?]object');

//checkforsigneddecimaldegreeswithoutNSEW,ifsoreturnitdirectly
if(typeofdmsStr==='number'&&isFinite(dmsStr))returnNumber(dmsStr);

//stripoffanysignorcompassdir'n&splitoutseparated/m/s
vardms=String(dmsStr).trim().replace(/^-/,'').replace(/[NSEW]$/i,'').split(/[^0-9.,]+/);
if(dms[dms.length-1]=='')dms.splice(dms.length-1);//fromtrailingsymbol

if(dms=='')returnNaN;

//andconverttodecimaldegrees...
switch(dms.length){
case3://interpret3-partresultasd/m/s
vardeg=dms[0]/1+dms[1]/60+dms[2]/3600;
break;
case2://interpret2-partresultasd/m
vardeg=dms[0]/1+dms[1]/60;
break;
case1://justd(possiblydecimal)ornon-separateddddmmss
vardeg=dms[0];
//checkforfixed-widthunseparatedformateg0033709W
//if(/[NS]/i.test(dmsStr))deg='0'+deg;//-normaliseN/Sto3-digitdegrees
//if(/[0-9]{7}/.test(deg))deg=deg.slice(0,3)/1+deg.slice(3,5)/60+deg.slice(5)/3600;
break;
default:
returnNaN;
}
if(/^-|[WS]$/i.test(dmsStr.trim()))deg=-deg;//take'-',westandsouthas-ve
returnNumber(deg);
}

/**
*Convertdecimaldegreestodeg/min/secformat
*-degree,prime,double-primesymbolsareadded,butsignisdiscarded,thoughnocompass
*directionisadded
*
*@private
*@param{Number}deg:Degrees
*@param{String}[format=dms]:Returnvalueas'd','dm','dms'
*@param{Number}[dp=0|2|4]:Noofdecimalplacestouse-default0fordms,2fordm,4ford
*@returns{String}degformattedasdeg/min/secsaccordingtospecifiedformat
*@throws{TypeError}degisanobject,perhapsDOMobjectwithout.value?
*/
Geo.toDMS=function(deg,format,dp){
if(typeofdeg=='object')thrownewTypeError('Geo.toDMS-degis[DOM?]object');
if(isNaN(deg))return'NaN';//giveuphereifwecan'tmakeanumberfromdeg

//defaultvalues
if(typeofformat=='undefined')format='dms';
if(typeofdp=='undefined'){
switch(format){
case'd':dp=4;break;
case'dm':dp=2;break;
case'dms':dp=0;break;
default:format='dms';dp=0;//beforgivingoninvalidformat
}
}

deg=Math.abs(deg);//(unsignedresultreadyforappendingcompassdir'n)

switch(format){
case'd':
d=deg.toFixed(dp);//rounddegrees
if(d<100)d='0'+d;//padwithleadingzeros
if(d<10)d='0'+d;
dms=d+'\u00B0';//addºsymbol
break;
case'dm':
varmin=(deg*60).toFixed(dp);//convertdegreestominutes&round
vard=Math.floor(min/60);//getcomponentdeg/min
varm=(min%60).toFixed(dp);//padwithtrailingzeros
if(d<100)d='0'+d;//padwithleadingzeros
if(d<10)d='0'+d;
if(m<10)m='0'+m;
dms=d+'\u00B0'+m+'\u2032';//addº,'symbols
break;
case'dms':
varsec=(deg*3600).toFixed(dp);//convertdegreestoseconds&round
vard=Math.floor(sec/3600);//getcomponentdeg/min/sec
varm=Math.floor(sec/60)%60;
vars=(sec%60).toFixed(dp);//padwithtrailingzeros
if(d<100)d='0'+d;//padwithleadingzeros
if(d<10)d='0'+d;
if(m<10)m='0'+m;
if(s<10)s='0'+s;
dms=d+'\u00B0'+m+'\u2032'+s+'\u2033';//addº,',"symbols
break;
}

returndms;
}

/**
*Convertnumericdegreestodeg/min/seclatitude(suffixedwithN/S)
*
*@param{Number}deg:Degrees
*@param{String}[format=dms]:Returnvalueas'd','dm','dms'
*@param{Number}[dp=0|2|4]:Noofdecimalplacestouse-default0fordms,2fordm,4ford
*@returns{String}Deg/min/seconds
*/
Geo.toLat=function(deg,format,dp){
varlat=Geo.toDMS(deg,format,dp);
returnlat==''?'':lat.slice(1)+(deg<0?'S':'N');//knockoffinitial'0'forlat!
}

/**
*Convertnumericdegreestodeg/min/seclongitude(suffixedwithE/W)
*
*@param{Number}deg:Degrees
*@param{String}[format=dms]:Returnvalueas'd','dm','dms'
*@param{Number}[dp=0|2|4]:Noofdecimalplacestouse-default0fordms,2fordm,4ford
*@returns{String}Deg/min/seconds
*/
Geo.toLon=function(deg,format,dp){
varlon=Geo.toDMS(deg,format,dp);
returnlon==''?'':lon+(deg<0?'W':'E');
}

/**
*Convertnumericdegreestodeg/min/secasabearing(0º..360º)
*
*@param{Number}deg:Degrees
*@param{String}[format=dms]:Returnvalueas'd','dm','dms'
*@param{Number}[dp=0|2|4]:Noofdecimalplacestouse-default0fordms,2fordm,4ford
*@returns{String}Deg/min/seconds
*/
Geo.toBrng=function(deg,format,dp){
deg=(Number(deg)+360)%360;//normalise-vevaluesto180º..360º
varbrng=Geo.toDMS(deg,format,dp);
returnbrng.replace('360','0');//justincaseroundingtookusupto360º!
}

/*-----------------------------------------------*/


综合可参考资料：

http://www.catalina-capri-25s.org/tech/latlongcalc.asp

http://www.movable-type.co.uk/scripts/latlong.html

http://www.cpearson.com/excel/LatLong.aspx

http://www.anycalculator.com/longitude.htm

http://www.movable-type.co.uk/scripts/latlong-vincenty.html

http://www.movable-type.co.uk/scripts/latlong-convert-coords.html
                                            

                                        

                        
                        内容来自用户分享和网络整理，不保证内容的准确性，如有侵权内容，可联系管理员处理 点击这里给我发消息 
                    

                    

                    

                        

                         标签： 
                                                

                        

                    


                

            


            

                
相关文章推荐

                

                
            

        

        

            

            
            

                

                    
新的分享

                    

                        
                    

                

            


            

                

                    
章节导航