geography class: great circle distance differs when using degrees vs radians

[sneakyimp]

I'm working on a geography class to help locate nearby businesses and calculate distances. To assist with generating SQL queries to locate businesses in a database, I've written functions get_bounding_box_rads and get_bounding_box_degrees. The idea is that you supply the lat and lon of a point and supply unitless values (you could specify miles or kilometers or furlongs, etc) to specify the radius of the sphere and the range and this function will return a bounding box of minimum and maximum latitudes and longitudes so that you can write an efficient query to retrieve records from your db. Efficient in the sense that the query has simple greater than / less than comparisons so that indexes can be used rather than some insane great circle distance calculation that results in a nasty table scan.

My problem is that the bounding box returns different longitudes (curiously, the latitudes seem identical so far) for the calculation performed in degrees vs the calculation performed in radians. I think this difference arises in the circle distance calculation which I borrowed from this wikipedia article, which says:

A formula that is accurate for all distances is the following special case of the Vincenty formula for an ellipsoid with equal major and minor axes:

https://wikimedia.org/api/rest_v1/media/math/render/svg/87cea288a5b6e80757bc81375c3b6a38a30a5184

Can anyone tell me the cause of the difference and how I might avoid it? If not, which is more accurate, degrees or radians? It's only about a 1 in 1000 difference, but that seems non trivial to me. I wonder if this is due to a 32-bit or 64-bit rounding error, or if perhaps I've made a mistake somewhere.

You should be able to simply run the following script and see what I mean

<?php
/**
 * A class to encapsulate and standardize geography-related functions.
 *
 * @author sneakyimp
 */
class OGH_geography {

	/**
	 * Radius of the earth in miles
	 * @var float
	 */
	const EARTH_RADIUS_MI = 3958.8;

	/**
	 * outer limit of latitude in radians
	 * @var float
	 */
	const LAT_LIMIT_RAD = M_PI/2;

	/**
	 * outer limit of latitude in degrees
	 * @var float
	 */
	const LAT_LIMIT_DEGREE = 90;

	/**
	 * outer limit of longitude in radians
	 * @var float
	 */
	const LON_LIMIT_RAD = M_PI;

	/**
	 * outer limit of longitude in degrees
	 * @var float
	 */
	const LON_LIMIT_DEGREE = 180;

	/**
	 * Calculates great circle distance between two points on a sphere
	 * @see https://en.wikipedia.org/wiki/Great-circle_distance
	 * @param float $lat1 Point 1 latitude in rads
	 * @param float $lon1 Point 1 longitude in rads
	 * @param float $lat2 Point 2 latitude in rads
	 * @param float $lon2 Point 2 longitude in rads
	 * @param float $radius The unitless radius of the sphere (could be km or miles)
	 * @return float unitless distance (depends on units in which radius is specified)
	 */
	public static function great_circle_distance_rads($lat1, $lon1, $lat2, $lon2, $radius) : float {
		// broadly, distance = radius * central_angle

		// using the 'special case of the vincenty formula we get
		// central_angle = arctan(sqrt((cos(lat2) * sin(lon1)) ** 2 + ((cos(lat1) * sin(lat2)) - (sin(lat1) * cos(lat2) * cos(abs(lon2 - lon1)))) ** 2) / (sin(lat1) * sin(lat2)) + (cos(lat1) * cos(lat2) * cos(abs(lon2 - lon1))))
		$coslat1 = cos($lat1);
		$coslat2 = cos($lat2);
		$sinlat1 = sin($lat1);
		$sinlat2 = sin($lat2);
		$delta_lon = abs($lon2 - $lon1);
		$sin_delta_lon = sin($delta_lon);
		$cos_delta_lon = cos($delta_lon);

		$is1 = ($coslat2 * $sin_delta_lon) ** 2;
		$is2 = (($coslat1 * $sinlat2) - ($sinlat1 * $coslat2 * $cos_delta_lon)) ** 2;
		$sq_root = sqrt($is1 + $is2);
		$denom = ($sinlat1 * $sinlat2) + ($coslat1 * $coslat2 * $cos_delta_lon);
		//$central_angle_radians = atan($sq_root / $denom);
		$central_angle_radians = atan2($sq_root, $denom);

		return $central_angle_radians * $radius;
	}

	/**
	 * Returns distance units per degree or latitude
	 * latitude, unlike longitude has constant distance per unit, so doesn't depend on location
	 * @param float $radius The unitless radius of the sphere
	 * @return float distance/length units per degree
	 */
	public static function distance_per_latitude_degree($radius) :float {
		// more subtle approach, use great circle distance calculation for points
		// of identical longitude, 1 deg of latitude apart
		$rads_per_deg = deg2rad(1);
		return self::great_circle_distance_rads(0, 0, $rads_per_deg, 0, $radius);
	}

	/**
	 * Returns distance units per radian of latitude
	 * latitude, unlike longitude has constant distance per unit, so doesn't depend on location
	 * @param float $radius The unitless radius of the sphere
	 * @return float distance/length units per radian
	 */
	public static function distance_per_latitude_rad($radius) :float {
		// more subtle approach, use great circle distance calculation for points
		// of identical longitude, 1 rad of latitude apart
		return self::great_circle_distance_rads(0, 0, 1, 0, $radius);
	}

	/**
	 * Returns distance units per radian of longitude
	 * distance longitudinally varies with the latitude, so we need additional param of $lat
	 * @param float $lat latitude, in radians, at which we want the longitudinal distance
	 * @param float $radius The unitless radius of the sphere
	 * @return float distance/length units per radian
	 */
	public static function distance_per_longitude_rad($lat, $radius) :float {
		// use great circle distance calculation for points
		// of identical latitude, 1 rad of longitude apart
		return self::great_circle_distance_rads($lat, 0, $lat, 1, $radius);
	}

	/**
	 * Returns distance units per degree of longitude
	 * distance longitudinally varies with the latitude, so we need additional param of $lat
	 * @param float $lat latitude, in degrees, at which we want the longitudinal distance
	 * @param float $radius The unitless radius of the sphere
	 * @return float distance/length units per degree
	 */
	public static function distance_per_longitude_degree($lat, $radius) :float {
		// use great circle distance calculation for points
		// of identical latitude, 1 degree of longitude apart
		$rads_per_deg = deg2rad(1);
		$lat_rads = deg2rad($lat);
		return self::great_circle_distance_rads($lat_rads, 0, $lat_rads, $rads_per_deg, $radius);
	}

	/**
	 * Given a point on the globe and a unitless range, will return upper and lower lat/lon in radians for
	 * a roughly square bounding box, useful for constructing a query to find nearby points
	 * @param float $radius unitless radius of globe
	 * @param float $lat latitude of center in rads
	 * @param float $lon longitude of center in rads
	 * @param float $range unitless range distance from center
	 * @return array associative array containing lat/lon floats in radians with keys min_lat, max_lat, min_lon, max_lon
	 */
	public static function get_bounding_box_rads($lat, $lon, $radius, $range) {
		if (!is_numeric($lat) || $lat > self::LAT_LIMIT_RAD || $lat < -self::LAT_LIMIT_RAD) {
			throw new \Exception("$lat is not a valid latitude");
		}
		if (!is_numeric($lon) || $lon > self::LON_LIMIT_RAD || $lon < -self::LON_LIMIT_RAD) {
			throw new \Exception("$lon is not a valid longitude");
		}
		if (!is_numeric($radius)) {
			throw new \Exception("$radius is not a valid radius");
		}
		if (!is_numeric($range)) {
			throw new \Exception("$radius is not a valid range");
		}
		// sanity check...this prob shouldn't be used for more than 1/20 of the globe
		$max_range = $radius/20;
		if ($range < 0 || $range > $max_range) {
			throw new \Exception("$range is out of bounds. Must be between 0 and $max_range");
		}


		$lat_dist_per_rad = self::distance_per_latitude_rad($radius);
		$lon_dist_per_rad = self::distance_per_longitude_rad($lat, $radius);

		$lat_adjust = $range/$lat_dist_per_rad;
		$lon_adjust = $range/$lon_dist_per_rad;

		$min_lat = $lat - $lat_adjust;
		// watch out for min_lat out of bounds
		$min_lat = ($min_lat < -self::LAT_LIMIT_RAD) ? -self::LAT_LIMIT_RAD : $min_lat;
		$max_lat = $lat + $lat_adjust;
		// watch out for max_lat out of bounds
		$max_lat = ($max_lat > self::LAT_LIMIT_RAD) ? self::LAT_LIMIT_RAD : $max_lat;

		$min_lon = $lon - $lon_adjust;
		// watch out for min_lon out of bounds
		$min_lon = ($min_lon <= -self::LON_LIMIT_RAD) ? ($min_lon + (2 * M_PI)) : $min_lon;
		$max_lon = $lon + $lon_adjust;
		// watch out for max_lon out of bounds
		$max_lon = ($max_lon > self::LON_LIMIT_RAD) ? ($max_lon - (2 * M_PI)) : $max_lon;

		return [
			'min_lat' => $min_lat,
			'max_lat' => $max_lat,
			'min_lon' => $min_lon,
			'max_lon' => $max_lon
		];
	}


	/**
	 * Given a point on the globe and a unitless range, will return upper and lower lat/lon in degrees for
	 * a roughly square bounding box, useful for constructing a query to find nearby points
	 * @param float $radius unitless radius of globe
	 * @param float $lat latitude of center in degrees
	 * @param float $lon longitude of center in degrees
	 * @param float $range unitless range distance from center
	 * @return array associative array containing lat/lon floats in degrees with keys min_lat, max_lat, min_lon, max_lon
	 */
	public static function get_bounding_box_degrees($lat, $lon, $radius, $range) {
		if (!is_numeric($lat) || $lat > self::LAT_LIMIT_DEGREE || $lat < -self::LAT_LIMIT_DEGREE) {
			throw new \Exception("$lat is not a valid latitude");
		}
		if (!is_numeric($lon) || $lon > self::LON_LIMIT_DEGREE || $lon < -self::LON_LIMIT_DEGREE) {
			throw new \Exception("$lon is not a valid longitude");
		}
		if (!is_numeric($radius)) {
			throw new \Exception("$radius is not a valid radius");
		}
		if (!is_numeric($range)) {
			throw new \Exception("$radius is not a valid range");
		}
		// sanity check...this prob shouldn't be used for more than 1/20 of the globe
		$max_range = $radius/20;
		if ($range < 0 || $range > $max_range) {
			throw new \Exception("$range is out of bounds. Must be between 0 and $max_range");
		}


		$lat_dist_per_degree = self::distance_per_latitude_degree($radius);
		$lon_dist_per_degree = self::distance_per_longitude_degree($lat, $radius);


		$lat_adjust = $range/$lat_dist_per_degree;
		$lon_adjust = $range/$lon_dist_per_degree;


		$min_lat = $lat - $lat_adjust;
		// watch out for min_lat out of bounds
		$min_lat = ($min_lat < -self::LAT_LIMIT_DEGREE) ? -self::LAT_LIMIT_DEGREE : $min_lat;
		$max_lat = $lat + $lat_adjust;
		// watch out for max_lat out of bounds
		$max_lat = ($max_lat > self::LAT_LIMIT_DEGREE) ? self::LAT_LIMIT_DEGREE : $max_lat;

		$min_lon = $lon - $lon_adjust;
		// watch out for min_lon out of bounds
		$min_lon = ($min_lon <= -self::LON_LIMIT_DEGREE) ? ($min_lon + 360) : $min_lon;
		$max_lon = $lon + $lon_adjust;
		// watch out for max_lon out of bounds
		$max_lon = ($max_lon > self::LON_LIMIT_DEGREE) ? ($max_lon - 360) : $max_lon;

		return [
				'min_lat' => $min_lat,
				'max_lat' => $max_lat,
				'min_lon' => $min_lon,
				'max_lon' => $max_lon
		];
	}
} // class OGH_geography


$lat = 34.078991;
$lon = -118.263804;

$lat_rad = deg2rad($lat);
$lon_rad = deg2rad($lon);

$box = OGH_geography::get_bounding_box_rads($lat_rad, $lon_rad, 3958.8, 5);
var_dump($box);
foreach($box as $key => $val) {
	$val = rad2deg($val);
	echo "$key $val\n";
}

$box = OGH_geography::get_bounding_box_degrees($lat, $lon, 3958.8, 5);
var_dump($box);
foreach($box as $key => $val) {
	echo "$key $val\n";
}

[Weedpacket]

Note that being off by 0.001° in longitude corresponds to about 111m measured at the equator, and less closer toward the poles. If you're worried about that difference you'd need to be taking into account the fact that the Earth isn't spherical.

The only real difference between the two code paths are the extra deg2rad conversions in distance_per_longitude_degree() and the conversion back to degrees at the end. The size of the error (0.00123783299685°, give or take) is the same across all longitudes at the given latitude. But it varies with latitude (from 0 at the equator to practically infinite at the poles) and the northern edge of the bounding box is at a different latitude from the southern.

[sneakyimp]

EDIT: I fat-fingered response before. Here is complete response.

Thanks, @Weedpacket for your response & insight here. It did occur to me to add some conditionals to try and avoid certain weirdness near the poles, but it now occurs to me that as your search approaches the poles very closely, there's not enough longitude available to extend to the specified range. E.g., if you specify a range of 100 miles and you get close enough to the pole, the longitude line circling the globe at that high latitude might only be a fraction of that 100-mile range, and my logic leads to longitude ranges that exceed 180 degrees.

E.g., these values:

$lat = 89.9;
$lon = -118.263804;
$box = OGH_geography::get_bounding_box_degrees($lat, $lon, 3958.8, 100);

yield a globe-circling longitude range:

min_lat 88.452698304711
max_lat 90
min_lon -587.51753836652
max_lon 350.98993036652

I suppose I'll stick in a sanity check

if (abs($max_lat - $min_lat) > 360) {
    throw new Exception('bounding box longitude out of range'); // or something
}

Weedpacket The only real difference between the two code paths are the extra deg2rad conversions

This is presumably done by multiplying/dividing by 2xPI. I'm guessing this results in some tiny difference in values, either caused by an inexact M_PI or by the limitations of 32 bit (or 64 bit) architecture and that tiny difference gets amplified by the great circle distance calculation?

Weedpacket the northern edge of the bounding box is at a different latitude from the southern.

Can you clarify what you mean here? I'm aware that distance-per-longitude will differ a tiny amount between one's selected search point, the northern longitude, and the southern latitude, unless your selected search point happens to be at precisely latitude zero.

I will note that the maximum latitude in my data set is zip code 99723:

99723,71.253861,-156.798246

Which yields this output when i specify a range of 100 miles:

// rad
min_lat 69.806559304711
max_lat 72.701162695289
min_lon -161.47625951744
max_lon -152.12023248256

// degrees
min_lat 69.806559304711
max_lat 72.701162695289
min_lon -161.30175506031
max_lon -152.29473693969

min_lon and max_lon differ by -0.17450445713001 degrees.

My function here:

$dist_per_deg =  OGH_geography::distance_per_longitude_degree(71.253861, 3958.8);
echo "dist per degree: $dist_per_deg\n";
$error = $dist_per_deg * 0.17450445713001;
echo "error: $error\n";

says that at latitude 71.253861, there are 22.204907031574 miles per degree of longitude and this means a difference between degrees vs rads of 3.9 miles.

dist per degree: 22.204907031574
error: 3.8748552471672

I think this is tolerable for a range of 100 miles, but I'm surprised at the size of the error.

[Weedpacket]

sneakyimp This is presumably done by multiplying/dividing by 2xPI

π/180

I'm going to have to trace through what's happening again, because there are a couple of things bothering me.

1. Lines of latitude are not great circles: measuring by $radius along a line of latitude generally won't get you a point that is $radius GC distance from your starting point and vice versa. The circumference of a latitude circle is just 2πrsin(lat), where r is the Earth's radius. It's pretty straightforward then to convert between angular measure and distance because the circumference of the circle is the conversion factor.
2. Contrariwise, all longitude circles are great circles — i.e., they all have the same radius r; again, converting between distances and angles is just a matter of using the fact that the circumference of such a circle is 2πr.
3. Couldn't roughly half your code be replaced by rad2deg(calculate_in_radians(deg2rad($degrees))) type constructions?

[sneakyimp]

Thanks, again, for your generous insights, @Weedpacket.

Weedpacket π/180

DOH yes I bricked on that one. My trig muscles are rusty.

Weedpacket 1 - Lines of latitude are not great circles...

I had not bothered to look up what the term 'great circle' meant, but see what you mean now. My first impulse in calculating distance per degree of longitude was to use some simple sin/cos calculation, but I doubted myself, thinking that was too easy. I thought 'but what about the curvature of the earth...' and rather thoughtlessly just used the great circle approach out of laziness. I've modified my distance_per_longitude_rad and distance_per_longitude_degree functions to use a simple formula as you suggest (but circumference_at_lat is 2πrcos(lat), yes?)

    public static function distance_per_longitude_rad($lat, $radius) :float {
       return $radius * cos($lat);
    }
    public static function distance_per_longitude_degree($lat, $radius) :float {
        return self::distance_per_longitude_rad(deg2rad($lat), $radius) * (M_PI / 180);
    }    

This adjustment appears to eliminate any difference when I run this code:

$lat = 34.078991;
$lon = -118.263804;

$lat_rad = deg2rad($lat);
$lon_rad = deg2rad($lon);

$box = OGH_geography::get_bounding_box_rads($lat_rad, $lon_rad, 3958.8, 5);
foreach($box as $key => $val) {
    $val = rad2deg($val);
    echo "$key $val\n";
}

$box = OGH_geography::get_bounding_box_degrees($lat, $lon, 3958.8, 5);
foreach($box as $key => $val) {
    echo "$key $val\n";
}

the output:

$ php foo.php
min_lat 34.006625915236
max_lat 34.151356084764
min_lon -118.35117339305
max_lon -118.17643460695

min_lat 34.006625915236
max_lat 34.151356084764
min_lon -118.35117339305
max_lon -118.17643460695

Weedpacket 2 - Contrariwise, all longitude circles are great circles...

It sounds like you're suggesting that distance_per_latitude_x functions should just divide the globe's circumference by 2PI or 360 and avoid the elaborate great circle calculation. I'll make this adjustment.

Weedpacket 3 - Couldn't roughly half your code be replaced by rad2deg...

This is probably true, and I'll probably tidy up here momentarily. My first run at this class was intended to provide a class that would expose convenient functions for calculations in degrees or rads at the top level of one's code, without the need to loop and clutter one's controllers with a bunch of rad2deg / deg2rad stuff. Nearly all geographic APIs (e.g. google maps) seem to express lat/lon in degrees but all the PHP functions use rads. 🤷I also wrote some static methods where it assumes you are talking about earth and you want miles so you don't have to supply any radius parameters. I'm sure it could use some cleanup, in any case.

In the end, I still wonder where the discrepancy arose. I suspect some tiny rounding error or lack of precision introduced by rad2deg/deg2rad was getting magnified by the great circle distance calculation somehow. The wikipedia page I linked shows various formulae, but says the one I used "is accurate for all distances". 🤔

[Weedpacket]

sneakyimp (but circumference_at_lat is 2πrcos(lat), yes?)

Uh, yeah: should be maximum at the equator and zero at the poles. I should have drawn myself something to look at while I was writing.

[NogDog]

C'mon...the Earth is flat, why are you wasting all this effort on spherical geometry? 😆

Seriously, though, if it's feasible for the project you're working on, I'd look into using the preferred geo-spatial extension, plug-in, library, whatever for the DBMS you're using. They'll provide some sort of "spatial" column type, and functions you can populate them with something like some_function([latitude], [longitude]), along with functions to get the distance between any two such values.

I did this about 6 years ago or so, in that case using a Neo4j graph database, setting up a node for each Zip Code in the US. (The data we wanted to compare just had zip codes, not lat/long.) It was a proof of concept thing that worked a treat, but we never ended up using it, and I've probably forgotten 90% of what I did now. 🙂

[sneakyimp]

NogDog C'mon...the Earth is flat, why are you wasting all this effort on spherical geometry? 😆

This is funny, and yet profoundly sad. I almost feel you should put a disclaimer on your joke, or we are likely to get some more QAnon weirdos dying in rocket crashes, trying to prove the world is flat.

I have been curious about the geospatial stuff, and it appears that both MySQL and MariaDB offer geometry type functionality, but it looks like a pretty big chunk of documentation -- and I'd point out that the census data i found for zip code lat/lon doesn't appear to offer any SQL data.

If you have any quick-start suggestions about getting acquainted with such stuff, I'd appreciate it.

[Weedpacket]

sneakyimp weirdos dying in rocket crashes, trying to prove the world is flat.

Hey, they're making the effort to test their convictions; give them credit for that.

I mean, I don't see how their rockets made of welded-together garbage cans are going to get them up to the airliner cruising altitudes they'd need to see the curvature for themselves, but....

[sneakyimp]

Weedpacket I mean, I don't see how their rockets made of welded-together garbage cans are going to get them up to the airliner cruising altitudes they'd need to see the curvature for themselves, but....

I don't expect these flat earthers are going to learn much, if anything, about the shape of earth beyond what Eratosthenes did 22 centuries ago. They may earn a Darwin award, though.

In other news, my business search (search for a record store near you!) is up and running on my site. The adjusted distance_per_longitude_x functions appear to have resolved any weirdness.