Last active
August 29, 2015 14:13
-
-
Save oliverheilig/3808a17393053054b015 to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| // Calcluate airline distance for two mercator points. | |
| // This approximation formula is sufficiently accurate for | |
| // our needs for distances of up to 600 km and | |
| // 80° latitude (the error is never more than 5% even for extreme values). | |
| using System; | |
| public class Program | |
| { | |
| public static void Main() | |
| { | |
| // our latitude/longitude values for test | |
| double lat1 = 49.0; | |
| double lon1 = 8.0; | |
| double lat2 = 50.0; | |
| double lon2 = 9.0; | |
| // convert to PTV Mercator, so assume we have Mercator as input | |
| double x1, y1, x2, y2; | |
| LatLonToSphereMercator(lat1, lon1, out x1, out y1); | |
| LatLonToSphereMercator(lat2, lon2, out x2, out y2); | |
| // calclate mercator distance with Pythagoras | |
| double mercDist = Math.Sqrt(Math.Pow(x2 - x1, 2) + Math.Pow(y2 - y1, 2)); | |
| // the center of the y values to apply the correction factor | |
| double centerY = (y1 + y2) / 2; | |
| // calculate the corresponding latitude | |
| double centerLat = LatFromMercator(centerY); | |
| // the latitude determines our correction factor | |
| double correctionFactor = Math.Cos(centerLat * Math.PI / 180.0); | |
| // apply the correction factor to real distance | |
| double realDist = mercDist * correctionFactor; | |
| Console.WriteLine("Distance Between {0}/{1} and {2}/{3}", x1, y1, x2, y2); | |
| Console.WriteLine("Pythagoras: {0}", mercDist); | |
| Console.WriteLine("Correction Factor: {0}", correctionFactor); | |
| Console.WriteLine("Real: {0}", realDist); | |
| } | |
| // Project a geographic coordinate on the (spherical) Mercator map. | |
| // Using the mean between major and minor axis here ("PTV standard"). | |
| // You could also use the major axis 6378137 ("Google standard"). | |
| public static void LatLonToSphereMercator( | |
| double latitude, double longitude, out double x, out double y) | |
| { | |
| x = 6371000.0 * longitude * Math.PI / 180; | |
| y = 6371000.0 * Math.Log(Math.Tan(Math.PI / 4 + latitude * Math.PI / 360)); | |
| } | |
| // Caclculate the latitude for a mercator Y | |
| public static double LatFromMercator(double y) | |
| { | |
| return (360 / Math.PI) * (Math.Atan(Math.Exp(y / 6371000.0)) - (Math.PI / 4)); | |
| } | |
| } |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment