Recently I needed to calculate distance from one point to another set of points in order to find the nearest point and its distance to the origin point. I opted to project all points to the
azimuthal equidistant projection. The most important property of this projection is that all distances from the center of the projection to all other points represents the true distance (assuming a spheric earth) to all other projected points. To calculate the distance you just have to calculate the euclidean distance and multiply it with the average earth radius (6371000.0 meters).
The formulas for the azimuthal equidistant projection can be found at
mathworld.wolfram.com. This formula can then be directly translated into F# like this:
open System
module AzimuthalEquidistantProjection =
let inline degToRad d = 0.0174532925199433 * d; // (1.0/180.0 * Math.PI) * d
let project centerlon centerlat lon lat =
// http://mathworld.wolfram.com/AzimuthalEquidistantProjection.html
// http://www.radicalcartography.net/?projectionref
let t:float = degToRad lat
let l:float = degToRad lon
let t1 = degToRad centerlat // latitude center of projection
let l0 = degToRad centerlon // longitude center of projection
let c = acos ((sin t1) * (sin t) + (cos t1) * (cos t) * (cos (l-l0)))
let k = c / (sin c)
let x = k * (cos t) * (sin (l-l0))
let y = k * ((cos t1) * (sin t) - (sin t1) * (cos t) * (cos (l-l0)))
(x, y)
I also implemented the azimuthal equidistant projection by using
units of measure. I create one set of measures to distinguish between degrees and radians and another one to distinguish between x and y.
open System
module AzimuthalEquidistantProjectionWithMeasures =
[<Measure>] type deg // degrees
[<Measure>] type rad // radians
[<Measure>] type x
[<Measure>] type y
type lon = float<x deg>
type lat = float<y deg>
let inline degToRad (d:float<'u deg>) = d*0.0174532925199433<rad/deg>
let cos (d:float<'u rad>) = Math.Cos(float d)
let sin (d:float<'u rad>) = Math.Sin(float d)
let project (centerlon:lon) (centerlat:lat) (lon:lon) (lat:lat) =
// http://mathworld.wolfram.com/AzimuthalEquidistantProjection.html
// http://www.radicalcartography.net/?projectionref
let t = degToRad lat
let l = degToRad lon
let t1 = degToRad centerlat // latitude center of projection
let l0 = degToRad centerlon // longitude center of projection
let c = acos ((sin t1) * (sin t) + (cos t1) * (cos t) * (cos (l-l0)))
let k = c / (sin c)
let x:float<x> = k * (cos t) * (sin (l-l0))
|> LanguagePrimitives.FloatWithMeasure
let y:float<y> = k * ((cos t1) * (sin t) - (sin t1) * (cos t) * (cos (l-l0)))
|> LanguagePrimitives.FloatWithMeasure
(x, y)
In the past I've used the
python OGR bindings and the
ESRI Projection Engine for my map projection needs but this time I needed a pure python implementation so I translated the above code and optimized it a bit by precalculating some values that we'll need when we project the points during the initialization of my projection class. A similarly optimized version in F# is on
fssnip.net.
from math import cos, sin, acos, radians
class Point(object):
def __init__(self,x,y):
self.x = x
self.y = y
class AzimuthalEquidistantProjection(object):
"""
http://mathworld.wolfram.com/AzimuthalEquidistantProjection.html
http://www.radicalcartography.net/?projectionref
"""
def __init__(self, center):
self.center = center
self.t1 = radians(center.y) ## latitude center of projection
self.l0 = radians(center.x) ## longitude center of projection
self.cost1 = cos(self.t1)
self.sint1 = sin(self.t1)
def project(self, point):
t = radians(point.y)
l = radians(point.x)
costcosll0 = cos(t) * cos(l-self.l0)
sint = sin(t)
c = acos ((self.sint1) * (sint) + (self.cost1) * costcosll0)
k = c / sin(c)
x = k * cos(t) * sin(l-self.l0)
y = k * (self.cost1 * sint - self.sint1 * costcosll0)
return Point(x, y)
import unittest
class Test_AzimuthalEquidistantProjection(unittest.TestCase):
def test_project(self):
p = AzimuthalEquidistantProjection(Point(0.0,0.0))
r = p.project(Point(1.0,1.0))
self.assertAlmostEqual(0.01745152022, r.x)
self.assertAlmostEqual(0.01745417858, r.y)
p = AzimuthalEquidistantProjection(Point(1.0,2.0))
r = p.project(Point(3.0,4.0))
self.assertAlmostEqual(0.03482860733, r.x)
self.assertAlmostEqual(0.03494898734, r.y)
p = AzimuthalEquidistantProjection(Point(-10.0001, 80.0001))
r = p.project(Point(7.935, 63.302))
self.assertAlmostEqual(0.1405127567, r.x)
self.assertAlmostEqual(-0.263406547, r.y)
if __name__ == '__main__':
unittest.main()
I've also implemented the straightforward projection code in
Julia. Julia is a promising language built for technical computing. It's still only on version 0.2 but a lot of
packages have already been created for the language.
function degToRad(d)
d*0.0174532925199433
end
function project(centerlon, centerlat, lon, lat)
t = degToRad(lat)
l = degToRad(lon)
t1 = degToRad(centerlat)
l0 = degToRad(centerlon)
c = acos(sin(t1) * sin(t) + cos(t1) * cos(t) * cos(l-l0))
k = c / sin(c)
x = k * cos(t) * sin(l-l0)
y = k * (cos(t1) * sin(t) - sin(t1) * cos(t) * cos(l-l0))
x, y
end
More Python, F#, Julia and geospatial code will be posted here so make sure to subscribe to my e-mail updates on the left. I also occasionally tweet as
@gissolved.
