distancespherical-geometry
I'm trying to study how to find the distance between two points if their latitudes and longitudes are given.
I have a question about the proof of part of the haversine formula given at Math Forum. It says that the length of chord AD, two points at the same latitude, lat1, on a sphere of radius 1, is
2*sin(dlon/2)*cos(lat1)
but I couldn't get how they obtained it. Could you help me?
Using the Pythagorean formula on positions given in latitude and longitude makes as little sense as, say, computing the area of a circle using the formula for a square: although it produces a number, there is no reason to suppose it ought to work.
Although at small scales any smooth surface looks like a plane, the accuracy of the Pythagorean formula depends on the coordinates used. When those coordinates are latitude and longitude on a sphere (or ellipsoid), we can expect that
Distances along lines of longitude will be reasonably accurate.
Distances along the Equator will be reasonably accurate.
All other distances will be erroneous, in rough proportion to the differences in latitude and longitude.
The error depends on the start and end point of the distance calculations. However, because both the sphere and ellipsoid have a circular symmetry around the axis, the error depends only on the difference of the longitudes, so to study this error we might as well take the point of origin to be at the Prime Meridian. Because both the sphere and ellipsoid are symmetric under a north-south reflection, we only need to study points of origin in the southern hemisphere. For any such point we may draw a contour map of the relative error, equal to [Pythagorean calculation] / [True distance].
The Pythagorean formula, using the mean radius of the earth, is
Pythagorean distance = 6371000. * Sqrt[dx^2 + dy^2]] * pi / 180 meters
where dx is the difference in longitudes and dy is the difference in latitudes, both in degrees. (The difference in longitude values is reduced modulo 360 to give the correct value of dx when crossing the antimeridian; not doing so would introduce artificially large errors that tell us nothing about the Pythagorean formula itself.)
The following plots show the relative error compared to the correct distance on the WGS 84 ellipsoid for latitudes from -70 to 0 in increments of 10 degrees. The horizontal coordinate is the difference in longitudes and the vertical coordinate is the latitude of the destination. Light regions have relatively small error: the contour lines are at 1, 1.01, 1.02, 1.05, 1.1, 1.2, 1.5, 2, etc. (The pure white areas in the corners are places where the error goes beyond the range of these contours.) The red dots show the point of origin.
The vertical white bands testify to the correctness of expectation (1): Pythagorean distances are accurate when there is a small difference in longitudes. The horizontal white bands at low latitudes confirm expectation (2): near the Equator, horizontal distances are reasonably accurate. Otherwise, as witnessed by the extensive darker regions, at all other distances the Pythagorean formula is bad.
We can make quantitative estimates of the maximum error attained for pairs of nearby points (within, say, a few hundred kilometers of each other). Scale--using an appropriate value for the radius--is true along the meridian but along a circle of latitude it errs approximately by the secant of the latitude. For example, at a latitude of 40 degrees the secant is 1.31, implying the Pythagorean formula will give distances about 31% too large in the east-west direction. (This is evident in the upper right contour plot, for a point of origin at -40 degrees latitude, where the region immediately east-west of the red dot lies between the 1.2 and 1.5 contours.) Short distances in all other directions will be too large by some amount between 0% and 31%; longer distances may err by even more (as the contour plots show).
With points as close as these we can ignore what type of spheroid the earth is and use the eqirectangular projection Pythagoras formula is all that is required.
The following two functions are in PHP.
function Equirectangular($lat1,$lng1,$lat2,$lng2){
$lng1 = abs($lng1);
$lng2 = abs($lng2);
$alpha = $lng2-$lng1;
$x = deg2rad($alpha) * cos(deg2rad($lat1+$lat2)/2);
$y = deg2rad($lat1-$lat2);
$R = 6372.8*1000; // gives d in km
$distance = sqrt($x*$x + $y*$y) * $R;
return $distance;
}
function Haversine($lat1,$lng1,$lat2,$lng2) {
$deltaLat = $lat2 - $lat1 ;
$deltaLng = $lng2 - $lng1 ;
$earthRadius = 6372.8 *1000; // 3959 in miles 6371 in meters.
$alpha = $deltaLat/2;
$beta = $deltaLng/2;
$a = sin(deg2rad($alpha)) * sin(deg2rad($alpha)) + cos(deg2rad($lat1)) * cos(deg2rad($lat2)) * sin(deg2rad($beta)) * sin(deg2rad($beta)) ;
$c = 2 * atan2(sqrt($a), sqrt(1-$a));
$distance = $earthRadius * $c;
return $distance;
}
As you can Equirectangular() is simpler than Haversine()and can be used for small distances. I have multiplied earth radius by 100 to convert to metres.
Result. Equirectangular 0.211513 m Haversine 0.211513 m The following code can be modified to suit your database
PDOPyth.php
Use PDO with parameters to query database using Equirectangular approximation formula for range.
Save exceptions to file.
02/03/2016
*/
require("dbinfo2.php");
//Create coordinates
$lat = 52;
$lng =-2;
$radius = 25;
//Connect to database
$dbh = new PDO("mysql:host=$host;dbname=$database", $username);
$dbh->setAttribute(PDO::ATTR_ERRMODE, PDO::ERRMODE_EXCEPTION);
try {
// Prepare statement
//SELECT (SQRT( POW(RADIANS(lng1 -lng2), 2) + POW(RADIANS(lat1 - lat2)* COS(RADIANS(lat1 - lat2)/2),2 ))*6372
$stmt1 = $dbh->prepare("SELECT name, lat, lng, ( SQRT( POW(RADIANS(? - lng), 2) + POW(RADIANS(? - lat)* COS(RADIANS(? - lat)/2),2 ))*6372 ) AS distance FROM gbstn HAVING distance < ? ORDER BY distance ");
// Assign parameters
$stmt1->bindParam(1,$lng);
$stmt1->bindParam(2,$lat);
$stmt1->bindParam(3,$lat);
$stmt1->bindParam(4,$radius);
//Execute query
$stmt1->setFetchMode(PDO::FETCH_ASSOC);
$stmt1->execute();
$data = array();
// Iterate through the rows, adding nodes for each
while($row = $stmt1->fetch()) {
$data[] = $row;
}
echo $lat ." ".$lng;
echo "<BR>" ;
echo json_encode($data);
}
catch(PDOException $e) {
echo "I'm sorry I'm afraid you can't do that.". $e->getMessage() ;// Remove or modify after testing
file_put_contents('PDOErrors.txt',date('[Y-m-d H:i:s]').", mapSelect.php, ". $e->getMessage()."\r\n", FILE_APPEND);
}
//Close the connection
$dbh = null;
?>
Best Answer
The radius, r, of the small circle joining all points at latitude, φ is
r = R cos φ
where R is the radius of the sphere. That simplifies to
r = cos φ
if we assume a "unit sphere" (R = 1) for convenience.
The chord length of a straight line, AD, joining two points on the same latitude is
AD = 2 r sin dλ/2
where dλ is the difference in longitude of A and D. Thus
AD = 2 R cos φ sin dλ/2
or
AD = 2 cos φ sin dλ/2
if R = 1