Calculate the azimuth between two points (from Aviation Formulary V1.45):
DELIMITER //
DROP FUNCTION IF EXISTS azimuth//
CREATE FUNCTION azimuth(lng1 DOUBLE, lat1 DOUBLE, lng2 DOUBLE, lat2 DOUBLE)
RETURNS DOUBLE
DETERMINISTIC
BEGIN
DECLARE az DOUBLE;
DECLARE lat_rad1 DOUBLE;
DECLARE lat_rad2 DOUBLE;
DECLARE dLon DOUBLE;
DECLARE var DOUBLE;
SET var = lng1;
SET lat_rad1 := RADIANS(lat1);
SET lat_rad2 := RADIANS(lat2);
SET dlon := RADIANS(lng2-lng1);
SET az := atan2(
sin(dLon)*cos(lat_rad2),
(cos(lat_rad1)*sin(lat_rad2)) - (sin(lat_rad1)*cos(lat_rad2)*cos(dLon))
);
SET az := (DEGREES(az) + 360) % 360;
RETURN az;
END//
DELIMITER ;
Calculate perpendicular bearing of an azimuth (note: MySQL was not correctly calculating values for 90° and 270°, so I hard coded the correct values):
DELIMITER //
DROP FUNCTION IF EXISTS perpendicularBearing//
CREATE FUNCTION perpendicularBearing(azimuth DOUBLE, side VARCHAR(5))
RETURNS DOUBLE
BEGIN
DECLARE perp_az DOUBLE;
SET perp_az = DEGREES(atan(1/tan(RADIANS(azimuth))*-1));
IF side = "left" THEN
SET perp_az = 180 + (180 -(180 - perp_az));
END IF;
IF azimuth = 90 THEN
SET perp_az = 0;
END IF;
IF azimuth = 270 THEN
SET perp_az = 0;
END IF;
IF azimuth = -90 THEN
SET perp_az = 0;
END IF;
IF azimuth = -270 THEN
SET perp_az = 0;
END IF;
RETURN perp_az;
END//
DELIMITER ;
Calculate the coordinates of a point given a start point, an azimuth, and a distance in feet (from Aviation Formulary V1.45):
DELIMITER //
DROP FUNCTION IF EXISTS destinationPoint//
CREATE FUNCTIOn destinationPoint(azimuth DOUBLE, distance DOUBLE, lat DOUBLE, lng DOUBLE)
RETURNS Point
BEGIN
DECLARE az_rad DOUBLE;
DECLARE lat_rad DOUBLE;
DECLARE lng_rad DOUBLE;
DECLARE dist_km DOUBLE;
DECLARE dist_rad DOUBLE;
DECLARE lat_dest DOUBLE;
DECLARE lng_dest DOUBLE;
DECLARE dest_point Point;
SET dist_km = distance * 0.0003048;
SET dist_rad = dist_km/6371;
SET lat_rad = RADIANS(lat);
SET lng_rad = RADIANS(lng);
SET az_rad = RADIANS(azimuth);
SET lat_dest = asin(sin(lat_rad)*cos(dist_rad) +
cos(lat_rad)*sin(dist_rad)*cos(az_rad) );
SET lng_dest = lng_rad + atan2(sin(az_rad)*sin(dist_rad)*cos(lat_rad),
cos(dist_rad)-sin(lat_rad)*sin(lat_rad));
SET lng_dest = (lng_dest+3*pi())%(2*pi()) - pi();
SET dest_point = Point(DEGREES(lng_dest),DEGREES(lat_dest));
RETURN dest_point;
END//
DELIMITER ;
Finally, Vincenty Distance calculations from bramsi at forge.mysql.com:
DELIMITER ;;
DROP FUNCTION IF EXISTS `vd`;;
CREATE FUNCTION `vd`(lng1 DOUBLE, lat1 DOUBLE, lng2 DOUBLE, lat2 DOUBLE, metric VARCHAR(2)) RETURNS DOUBLE
DETERMINISTIC
COMMENT 'Vincenty Distance WGS-84 http://code.google.com/p/geopy/'
BEGIN
DECLARE gcdx DOUBLE;
DECLARE lng_rad1 DOUBLE;
DECLARE lat_rad1 DOUBLE;
DECLARE lng_rad2 DOUBLE;
DECLARE lat_rad2 DOUBLE;
DECLARE wgs84_major DOUBLE;
DECLARE wgs84_minor DOUBLE;
DECLARE wgs84_flattening DOUBLE;
DECLARE delta_lng DOUBLE;
DECLARE reduced_lat1 DOUBLE;
DECLARE reduced_lat2 DOUBLE;
DECLARE sin_reduced1 DOUBLE;
DECLARE cos_reduced1 DOUBLE;
DECLARE sin_reduced2 DOUBLE;
DECLARE cos_reduced2 DOUBLE;
DECLARE lambda_lng DOUBLE;
DECLARE lambda_prime DOUBLE;
DECLARE iter_limit INT;
DECLARE sin_lambda_lng DOUBLE;
DECLARE cos_lambda_lng DOUBLE;
DECLARE sin_sigma DOUBLE;
DECLARE cos_sigma DOUBLE;
DECLARE sigma DOUBLE;
DECLARE sin_alpha DOUBLE;
DECLARE cos_sq_alpha DOUBLE;
DECLARE cos2_sigma_m DOUBLE;
DECLARE C DOUBLE;
DECLARE u_sq DOUBLE;
DECLARE A DOUBLE;
DECLARE B DOUBLE;
DECLARE delta_sigma DOUBLE;
SET lng_rad1 := RADIANS(lng1);
SET lat_rad1 := RADIANS(lat1);
SET lng_rad2 := RADIANS(lng2);
SET lat_rad2 := RADIANS(lat2);
SET wgs84_major := 6378.137;
SET wgs84_minor := 6356.7523142;
SET wgs84_flattening := 1 / 298.257223563;
SET delta_lng := lng_rad2 - lng_rad1;
SET reduced_lat1 := atan((1 - wgs84_flattening) * tan(lat_rad1));
SET reduced_lat2 := atan((1 - wgs84_flattening) * tan(lat_rad2));
SET sin_reduced1 := sin(reduced_lat1);
SET cos_reduced1 := cos(reduced_lat1);
SET sin_reduced2 := sin(reduced_lat2);
SET cos_reduced2 := cos(reduced_lat2);
SET lambda_lng := delta_lng;
SET lambda_prime := 2 * pi();
SET iter_limit = 20;
WHILE abs(lambda_lng - lambda_prime) > pow(10, -11) and iter_limit > 0 DO
SET sin_lambda_lng := sin(lambda_lng);
SET cos_lambda_lng := cos(lambda_lng);
SET sin_sigma := sqrt(pow((cos_reduced2 * sin_lambda_lng), 2) +
pow((cos_reduced1 * sin_reduced2 - sin_reduced1 *
cos_reduced2 * cos_lambda_lng), 2));
IF sin_sigma = 0 THEN
RETURN 0;
END IF;
SET cos_sigma := (sin_reduced1 * sin_reduced2 +
cos_reduced1 * cos_reduced2 * cos_lambda_lng);
SET sigma := atan2(sin_sigma, cos_sigma);
SET sin_alpha := cos_reduced1 * cos_reduced2 * sin_lambda_lng / sin_sigma;
SET cos_sq_alpha := 1 - pow(sin_alpha, 2);
IF cos_sq_alpha != 0 THEN
SET cos2_sigma_m := cos_sigma - 2 * (sin_reduced1 * sin_reduced2 /
cos_sq_alpha);
ELSE
SET cos2_sigma_m := 0.0;
END IF;
SET C := wgs84_flattening / 16.0 * cos_sq_alpha * (4 + wgs84_flattening * (4 - 3 * cos_sq_alpha));
SET lambda_prime := lambda_lng;
SET lambda_lng := (delta_lng + (1 - C) * wgs84_flattening * sin_alpha *
(sigma + C * sin_sigma *
(cos2_sigma_m + C * cos_sigma *
(-1 + 2 * pow(cos2_sigma_m, 2)))));
SET iter_limit := iter_limit - 1;
END WHILE;
IF iter_limit = 0 THEN
RETURN NULL;
END IF;
SET u_sq := cos_sq_alpha * (pow(wgs84_major, 2) - pow(wgs84_minor, 2)) / pow(wgs84_minor, 2);
SET A := 1 + u_sq / 16384.0 * (4096 + u_sq * (-768 + u_sq *
(320 - 175 * u_sq)));
SET B := u_sq / 1024.0 * (256 + u_sq * (-128 + u_sq * (74 - 47 * u_sq)));
SET delta_sigma := (B * sin_sigma *
(cos2_sigma_m + B / 4. *
(cos_sigma * (-1 + 2 * pow(cos2_sigma_m, 2)) -
B / 6. * cos2_sigma_m * (-3 + 4 * pow(sin_sigma, 2)) *
(-3 + 4 * pow(cos2_sigma_m, 2)))));
SET gcdx := wgs84_minor * A * (sigma - delta_sigma);
IF metric = 'km' THEN
RETURN gcdx;
ELSEIF metric = 'mi' THEN
RETURN gcdx * 0.621371192;
ELSEIF metric = 'nm' THEN
RETURN gcdx / 1.852;
ELSE
RETURN gcdx;
END IF;
END;;
DELIMITER ;
